## 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

The signal
has the *z* transform

where the last step holds for
, which is
true whenever
. Thus, the transfer function consists of a single pole at
, and it exists for .^{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)); |

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