## Stability Revisited

As defined earlier in §5.6 (page ), a filter is said to be*stable*if its impulse response decays to 0 as 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 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

*z*transform

^{9.1}Now consider what happens when we let become greater than 1. The pole of moves outside the unit circle, and the impulse response has an exponentially

*increasing*amplitude. (Note .) Thus, the definition of stability is violated. Since the

*z*transform exists only for , we see that 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:

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 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.3}Here, we will settle for a simple recipe for computing the reflection coefficients from the transfer-function denominator polynomial . 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 denote the th-order denominator polynomial of the recursive filter transfer function :We have introduced the new subscript because the step-down procedure is defined recursively in polynomial order. We will need to keep track of polynomials orders between 1 and . In addition to the denominator polynomial , we need its

*flip*:

The recursion begins by setting the th reflection coefficient to . If , the recursion halts prematurely, and the filter is declared unstable. (Equivalently, the polynomial is declared non-minimum phase, as defined in Chapter 11.) Otherwise, if , the polynomial order is decremented by 1 to yield as follows (recall that is

*monic*):

Next is set to , and the recursion continues until is reached, or is found for some . Whenever , 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.

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 |

% 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)); |

**Next Section:**

Bandwidth of One Pole

**Previous Section:**

Graphical Phase Response Calculation