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.1Now 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.3Here, 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


In addition to the denominator polynomial , we need its
flip:
The recursion begins by setting the




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