Testing Filter Stability in Matlab

Figure 8.6 gives a listing of a matlab function stabilitycheck for testing the stability of a digital filter using the Durbin step-down recursion. Figure 8.7 lists a main program for testing stabilitycheck against the more prosaic method of factoring the transfer-function denominator and measuring the pole radii. The Durbin recursion is far faster than the method based on root-finding.

Figure 8.6: Matlab function for testing digital filter stability by computing its reflection coefficients.

 

function [stable] = stabilitycheck(A);

N = length(A)-1; % Order of A(z)
stable = 1;      % stable unless shown otherwise
A = A(:);        % make sure it's a column vector
for i=N:-1:1
  rci=A(i+1);
  if abs(rci) >= 1
    stable=0;
    return;
  end
  A = (A(1:i) - rci * A(i+1:-1:2))/(1-rci^2);
% disp(sprintf('A[%d]=',i)); A(1:i)'
end

Figure 8.7: Test program (matlab) for function stabilitycheck.

 

% TSC - test function stabilitycheck, comparing against
%       pole radius computation

N = 200; % polynomial order
M = 20; % number of random polynomials to generate
disp('Random polynomial test');
nunstp = 0; % count of unstable A polynomials
sctime = 0; % total time in stabilitycheck()
rftime = 0; % total time computing pole radii
for pol=1:M
  A = [1; rand(N,1)]'; % random polynomial
  tic; 
  stable = stabilitycheck(A); 
  et=toc;  % Typ. 0.02 sec Octave/Linux, 2.8GHz Pentium
  sctime = sctime + et;
  % Now do it the old fashioned way
  tic; 
  poles = roots(A); % system poles
  pr = abs(poles);  % pole radii
  unst = (pr >= 1); % bit vector
  nunst = sum(unst);% number of unstable poles
  et=toc; % Typ. 2.9 sec Octave/Linux, 2.8GHz Pentium
  rftime = rftime + et;
  if stable, nunstp = nunstp + 1; end
  if (stable & nunst>0) | (~stable & nunst==0)
    error('*** stabilitycheck() and poles DISAGREE ***'); 
  end
end
disp(sprintf(...
     ['Out of %d random polynomials of order %d,',...
      ' %d were unstable'], M,N,nunstp));

When run in Octave over Linux 2.4 on a 2.8 GHz Pentium PC, the Durbin recursion is approximately 140 times faster than brute-force root-finding, as measured by the program listed in Fig.8.7.