Free Books
   Free Books    Introduction to Digital Filters  

Stability Revisited

As defined earlier in §5.6 (page [*]), a filter is said to be stable if its impulse response $ h(n)$ decays to 0 as $ n$ goes to infinity. In terms of poles and zeros, an irreducible filter transfer function is stable if and only if all its poles are inside the unit circle in the $ z$ plane (as first discussed in §6.8.6). This is because the transfer function is the z transform of the impulse response, and if there is an observable (non-canceled) pole outside the unit circle, then there is an exponentially increasing component of the impulse response. To see this, consider a causal impulse response of the form

$\displaystyle h(n) = R^n e^{j\omega nT}, \qquad n=0,1,2,\ldots\,. $

This signal is a damped complex sinusoid when $ 0 < R < 1$. It oscillates with zero-crossing rate $ 2\omega/2\pi=\omega/\pi$ zeros per second, and it has an exponentially decaying amplitude envelope. If $ R>1$, then the amplitude envelope increases exponentially as $ R^n$.

The signal $ h(n)=R^n e^{j\omega n T}$ has the z transform

\begin{eqnarray*} H(z) &=& \sum_{n=0}^\infty R^n e^{j\omega nT} z^{-n}\\ &=& \... ...ga T}z^{-1}\right)^{n}\\ &=& \frac{1}{1-Re^{j\omega T}z^{-1}}, \end{eqnarray*}

where the last step holds for $ \left\vert R e^{j\omega T}z^{-1}\right\vert<1$, which is true whenever $ R < \left\vert z\right\vert$. Thus, the transfer function consists of a single pole at $ z = Re^{j\omega T}$, and it exists for $ \vert z\vert>R$.9.1Now consider what happens when we let $ R$ become greater than 1. The pole of $ H(z)$ moves outside the unit circle, and the impulse response has an exponentially increasing amplitude. (Note $ \left\vert h(n)\right\vert = R^n$.) Thus, the definition of stability is violated. Since the z transform exists only for $ \left\vert z\right\vert > R$, we see that $ R\geq 1$ implies that the z transform no longer exists on the unit circle, so that the frequency response becomes undefined!

The above one-pole analysis shows that a one-pole filter is stable if and only if its pole is inside the unit circle. In the case of an arbitrary transfer function, inspection of its partial fraction expansion (§6.8) shows that the behavior near any pole approaches that of a one-pole filter consisting of only that pole. Therefore, all poles must be inside the unit circle for stability.

In summary, we can state the following:

$\textstyle \parbox{0.8\textwidth}{A necessary and sufficient condition for the ... ...oles of its irreducible transfer function lie strictly inside the unit circle.}$
Isolated poles on the unit circle may be called marginally stable. The impulse response component corresponding to a single pole on the unit circle never decays, but neither does it grow.9.2 In physical modeling applications, marginally stable poles occur often in lossless systems, such as ideal vibrating string models [86].

Computing Reflection Coefficients to
Check Filter Stability

Since we know that a recursive filter is stable if and only if all its poles have magnitude less than 1, an obvious method for checking stability is to find the roots of the denominator polynomial $ A(z)$ in the filter transfer function [Eq.$ \,$(7.4)]. If the moduli of all roots are less than 1, the filter is stable. This test works fine for low-order filters (e.g., on the order of 100 poles or less), but it may fail numerically at higher orders because the roots of a polynomial are very sensitive to round-off error in the polynomial coefficients [62]. It is therefore of interest to use a stability test that is faster and more reliable numerically than polynomial root-finding. Fortunately, such a test exists based on the filter reflection coefficients.

It is a mathematical fact [48] that all poles of a recursive filter are inside the unit circle if and only if all its reflection coefficients (which are always real) are strictly between -1 and 1. The full theory associated with reflection coefficients is beyond the scope of this book, but can be found in most modern treatments of linear prediction [48,47] or speech modeling [92,19,69]. An online derivation appears in [86].9.3Here, we will settle for a simple recipe for computing the reflection coefficients from the transfer-function denominator polynomial $ A(z)$. This recipe is called the step-down procedure, Schur-Cohn stability test, or Durbin recursion [48], and it is essentially the same thing as the Schur recursion (for allpass filters) or Levinson algorithm (for autocorrelation functions of autoregressive stochastic processes) [38].

Step-Down Procedure

Let $ A_N(z)$ denote the $ N$th-order denominator polynomial of the recursive filter transfer function $ H(z)=B(z)/A_N(z)$:

$\displaystyle A_N(z) \isdef 1 + a_{N,1}\,z^{-1}+ a_{N,2}\, z^{-2}+ \cdots + a_{N,N-1}\,z^{-(N-1)} + a_{N,N}\,z^{-N} \protect$ (9.4)

We have introduced the new subscript $ N$ because the step-down procedure is defined recursively in polynomial order. We will need to keep track of polynomials orders between 1 and $ N$.

In addition to the denominator polynomial $ A_N(z)$, we need its flip:

$\displaystyle \tilde{A}_N(z)$ $\displaystyle \isdef$ $\displaystyle z^{-N}A_N(1/z)$  
  $\displaystyle =$ $\displaystyle a_{N,N} + a_{N,N-1}\,z^{-1}+ a_{N,N-2}\, z^{-2}+ \cdots$  
    $\displaystyle \quad + a_{N,2}\,z^{-(N-2)} + a_{N,1}\,z^{-(N-1)} + z^{-N} \protect$ (9.5)

The recursion begins by setting the $ N$th reflection coefficient to $ k_N = a_{N,N}$. If $ \left\vert k_N\right\vert\geq 1$, the recursion halts prematurely, and the filter is declared unstable. (Equivalently, the polynomial $ A_N(z)$ is declared non-minimum phase, as defined in Chapter 11.)

Otherwise, if $ \left\vert k_N\right\vert<1$, the polynomial order is decremented by 1 to yield $ A_{N-1}(z)$ as follows (recall that $ A_N(z)$ is monic):

$\displaystyle A_{N-1}(z) = \frac{A_N(z) - k_N \tilde{A}_N(z)}{1-k_N^2} = \frac{A_N(z) - z^{-N} k_N A_N(1/z)}{1-k_N^2} \protect$ (9.6)

Next $ k_{N-1}$ is set to $ a_{N-1,N-1}$, and the recursion continues until $ k_1=a_{1,1}$ is reached, or $ \left\vert k_i\right\vert\geq 1$ is found for some $ i$.

Whenever $ \left\vert k_m\right\vert=1$, the recursion halts prematurely, and the filter is usually declared unstable (at best it is marginally stable, meaning that it has at least one pole on the unit circle).

Note that the reflection coefficients can also be used to implement the digital filter in what are called lattice or ladder structures [48]. Lattice/ladder filters have superior numerical properties relative to direct-form filter structures based on the difference equation. As a result, they can be very important for fixed-point implementations such as in custom VLSI or low-cost (fixed-point) signal processing chips. Lattice/ladder structures are also a good point of departure for computational physical models of acoustic systems such as vibrating strings, wind instrument bores, and the human vocal tract [81,16,48].

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.

Next Section:
Bandwidth of One Pole
Previous Section:
Graphical Phase Response Calculation