Free Books    Introduction to Digital Filters  

Pole-Zero Analysis

This chapter discusses pole-zero analysis of digital filters. Every digital filter can be specified by its poles and zeros (together with a gain factor). Poles and zeros give useful insights into a filter's response, and can be used as the basis for digital filter design. This chapter additionally presents the Durbin step-down recursion for checking filter stability by finding the reflection coefficients, including matlab code.

Going back to Eq.$ \,$(6.5), we can write the general transfer function for the recursive LTI digital filter as

$\displaystyle H(z) = g\frac{1 + \beta_1 z^{-1}+ \cdots + \beta_M z^{-M}}{1 + a_1 z^{-1}+ \cdots + a_N z^{-N}} \protect$ (9.1)

which is the same as Eq.$ \,$(6.5) except that we have factored out the leading coefficient $ b_0$ in the numerator (assumed to be nonzero) and called it g. (Here $ \beta_i \isdeftext b_i/b_0$.) In the same way that $ z^2 + 3z + 2$ can be factored into $ (z + 1)(z + 2)$, we can factor the numerator and denominator to obtain

$\displaystyle H(z) = g\frac{(1-q_1z^{-1})(1-q_2z^{-1})\cdots(1-q_Mz^{-1})}{(1-p_1z^{-1})(1-p_2z^{-1})\cdots(1-p_Nz^{-1})}. \protect$ (9.2)

Assume, for simplicity, that none of the factors cancel out. The (possibly complex) numbers $ \{q_1,\ldots,q_M\}$ are the roots, or zeros, of the numerator polynomial. When $ z$ is set to any of these values, the transfer function evaluates to 0. For this reason, the numerator roots $ q_i$ are called the zeros of the filter. In other words, the zeros of the numerator of an irreducible transfer-function are called the zeros of the transfer-function. Similarly, when $ z$ approaches any root of the denominator polynomial, the magnitude of the transfer function approaches infinity. Consequently, the denominator roots $ \{p_1, \ldots, p_N\}$ are called the poles of the filter.

The term ``pole'' makes sense when one plots the magnitude of $ H(z)$ as a function of z. Since $ z$ is complex, it may be taken to lie in a plane (the $ z$ plane). The magnitude of $ H(z)$ is real and therefore can be represented by distance above the $ z$ plane. The plot appears as an infinitely thin surface spanning in all directions over the $ z$ plane. The zeros are the points where the surface dips down to touch the $ z$ plane. At high altitude, the poles look like thin, well, ``poles'' that go straight up forever, getting thinner the higher they go.

Notice that the $ M+1$ feedforward coefficients from the general difference quation, Eq.$ \,$(5.1), give rise to $ M$ zeros. Similarly, the $ N$ feedback coefficients in Eq.$ \,$(5.1) give rise to $ N$ poles. Recall that we defined the filter order as the maximum of $ N$ and $ M$ in Eq.$ \,$(6.5). Therefore, the filter order equals the number of poles or zeros, whichever is greater.

Filter Order = Transfer Function Order

Recall that the order of a polynomial is defined as the highest power of the polynomial variable. For example, the order of the polynomial $ p(x)=1+2x+3x^2$ is 2. From Eq.$ \,$(8.1), we see that $ M$ is the order of the transfer-function numerator polynomial in $ z^{-1}$. Similarly, $ N$ is the order of the denominator polynomial in $ z^{-1}$.

A rational function is any ratio of polynomials. That is, $ R(z)$ is a rational function if it can be written as

$\displaystyle R(z)\eqsp \frac{P(z)}{Q(z)} $

for finite-order polynomials $ P(z)$ and $ Q(z)$. The order of a rational function is defined as the maximum of its numerator and denominator polynomial orders. As a result, we have the following simple rule:
$\textstyle \parbox{0.9\textwidth}{\emph{The order of an LTI filter is the order of its transfer function.}}$

It turns out the transfer function can be viewed as a rational function of either $ z^{-1}$ or $ z$ without affecting order. Let $ K=\max\{M,N\}$ denote the order of a general LTI filter with transfer function $ H(z)$ expressible as in Eq.$ \,$(8.1). Then multiplying $ H(z)$ by $ z^K/z^K$ gives a rational function of $ z$ (as opposed to $ z^{-1}$) that is also order $ K$ when viewed as a ratio of polynomials in $ z$. Another way to reach this conclusion is to consider that replacing $ z$ by $ z^{-1}$ is a conformal map [57] that inverts the $ z$-plane with respect to the unit circle. Such a transformation clearly preserves the number of poles and zeros, provided poles and zeros at $ z=\infty$ and $ z=0$ are either both counted or both not counted.

Graphical Computation of
Amplitude Response from
Poles and Zeros

Now consider what happens when we take the factored form of the general transfer function, Eq.$ \,$(8.2), and set $ z$ to $ e^{j\omega T}$ to get the frequency response in factored form:

$\displaystyle H(e^{j\omega T}) = g\frac{(1-q_1e^{-j\omega T})(1-q_2e^{-j\omega ... ...ga T})}{(1-p_1e^{-j\omega T})(1-p_2e^{-j\omega T})\cdots(1-p_Ne^{-j\omega T})} $

As usual for the frequency response, we prefer the polar form for this expression. Consider first the amplitude response $ G(\omega) \isdeftext \left\vert H(e^{j\omega T})\right\vert$.
$\displaystyle G(\omega)$ $\displaystyle =$ $\displaystyle \left\vert g\right\vert\frac{\left\vert 1-q_1e^{-j\omega T}\right... ... 1-p_2e^{-j\omega T}\right\vert\cdots\left\vert 1-p_Ne^{-j\omega T}\right\vert}$  
  $\displaystyle =$ $\displaystyle \left\vert g\right\vert \frac{\left\vert e^{-jM\omega T}\right\ve... ...vert e^{j\omega T}-p_2\right\vert\cdots\left\vert e^{j\omega T}-p_N\right\vert}$  
  $\displaystyle =$ $\displaystyle \left\vert g\right\vert \frac{\left\vert e^{j\omega T}-q_1\right\... ...\omega T}-p_2\right\vert\cdots\left\vert e^{j\omega T}-p_N\right\vert} \protect$ (9.3)

In the complex plane, the number $ z = x + jy$ is plotted at the coordinates $ (x, y)$ [84]. The difference of two vectors $ u = x_1 + jy_1$ and $ v = x_2 + jy_2$ is $ u - v = (x_1 - x_2) + j(y_1 - y_2)$, as shown in Fig.8.1. Translating the origin of the vector $ u-v$ to the tip of $ v$ shows that $ u-v$ is an arrow drawn from the tip of $ v$ to the tip of $ u$. The length of a vector is unaffected by translation away from the origin. However, the angle of a translated vector must be measured relative to a translated copy of the real axis. Thus the term $ e^{j\omega T} - q_i$ may be drawn as an arrow from the $ i$th zero to the point $ e^{j\omega T}$ on the unit circle, and $ e^{j\omega T} - p_i$ is an arrow from the $ i$th pole. Therefore, each term in Eq.$ \,$(8.3) is the length of a vector drawn from a pole or zero to a single point on the unit circle, as shown in Fig.8.2 for two poles and two zeros. In summary:

$\textstyle \parbox{0.8\textwidth}{% The frequency response magnitude (amplitude... ... by the product of lengths of vectors drawn from the poles to $e^{j\omega T}$.}$

Figure 8.1: Treatment of complex numbers as vectors in a plane.
\begin{figure}\input fig/kfig2p12.pstex_t \end{figure}

Figure 8.2: Measurement of amplitude response from a pole-zero diagram. A pole is represented in the complex plane by `X'; a zero, by `O'.
\begin{figure}\input fig/kfig2p13.pstex_t \end{figure}

For example, the dc gain is obtained by multiplying the lengths of the lines drawn from all poles and zeros to the point $ z = 1$. The filter gain at half the sampling rate is the product of the lengths of these lines when drawn to the point $ z = -1$. For an arbitrary frequency $ f$ Hz, we draw arrows from the poles and zeros to the point $ z = e^{j2\pi fT}$. Thus, at the frequency where the arrows in Fig.8.2 join, (which is slightly less than one-eighth the sampling rate) the gain of this two-pole two-zero filter is $ G(\omega) = (d_1d_2)/(d_3d_4)$. Figure 8.3 gives the complete amplitude response for the poles and zeros shown in Fig.8.2. Before looking at that, it is a good exercise to try sketching it by inspection of the pole-zero diagram. It is usually easy to sketch a qualitatively accurate amplitude-response directly from the poles and zeros (to within a scale factor).

Figure: Amplitude response obtained by traversing the entire upper semicircle in Fig.8.2 with the point $ e^{j\omega T}$. The point of the amplitude obtained in that figure is marked by a heavy dot. For real filters, precisely the same curve is obtained if the lower half of the unit circle is traversed, since $ G(\omega ) = G( - \omega )$. Thus, plotting the response over positive frequencies only is sufficient for real filters.
\begin{figure}\input fig/kfig2p14b.pstex_t \end{figure}

Graphical Phase Response Calculation

The phase response is almost as easy to evaluate graphically as is the amplitude response:

\begin{eqnarray*} \Theta(\omega) &\isdef & \angle \left\{g \frac{(1-q_1e^{-j\ome... ...gle(e^{j\omega T}-p_2)-\cdots-\angle(e^{j\omega T}-p_N) \protect \end{eqnarray*}

If $ g$ is real, then $ \angle g$ is either 0 or $ \pi $. Terms of the form $ e^{j\omega T}-z$ can be interpreted as a vector drawn from the point $ z$ to the point $ e^{j\omega T}$ in the complex plane. The angle of $ e^{j\omega T}-z$ is the angle of the constructed vector (where a vector pointing horizontally to the right has an angle of 0). Therefore, the phase response at frequency $ f$ Hz is again obtained by drawing lines from all the poles and zeros to the point $ e^{j2\pi f T}$, as shown in Fig.8.4. The angles of the lines from the zeros are added, and the angles of the lines from the poles are subtracted. Thus, at the frequency $ \omega$ the phase response of the two-pole two-zero filter in the figure is $ \Theta(\omega) = \theta_1+\theta_2-\theta_3-\theta_4$.

Figure 8.4: Measurement of phase response from a pole-zero diagram.
\begin{figure}\input fig/kfig2p15.pstex_t \end{figure}

Note that an additional phase of $ (N - M)2\pi fT$ radians appears when the number of poles is not equal to the number of zeros. This factor comes from writing the transfer function as

$\displaystyle H(z) = gz^{(N-M)}\frac{(z-q_1)(z-q_2)\cdots(z-q_M)}{(z-p_1)(z-p_2)\cdots(z-p_N)} $

and may be thought of as arising from $ N - M$ additional zeros at $ z=0$ when $ N > M$, or $ M - N$ poles at $ z=0$ when $ M>N$. Strictly speaking, every digital filter has an equal number of poles and zeros when those at $ z=0$ and $ z=\infty$ are counted. It is customary, however, when discussing the number of poles and zeros a filter has, to neglect these, since they correspond to pure delay and do not affect the amplitude response. Figure 8.5 gives the phase response for this two-pole two-zero example.

Figure 8.5: Phase response obtained from Fig.8.4 for positive frequencies. The point of the phase response corresponding to the arrows in that figure is marked by a heavy dot. For real filters, the phase response is odd ( $ \Theta (-\omega ) = -\Theta (\omega )$), so the curve shown here may be reflected through 0 and negated to obtain the plot for negative frequencies.
\begin{figure}\input fig/kfig2p16.pstex_t \end{figure}

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.

Bandwidth of One Pole

A typical formula relating 3-dB bandwidth $ B$ (in Hz) to the pole radius $ R\in[0,1)$ is

$\displaystyle \zbox {R = e^{- \pi B T}} \protect$ (9.7)

where $ T$ denotes the sampling interval as usual. Solving for $ B$ gives

$\displaystyle B = - \frac{\ln(R)}{\pi T}. $

In §E.6 (see also §B.1.3), it is shown that $ B$ is the 3-dB bandwidth or half-power bandwidth, in the limit, as the sampling rate goes to infinity.

Time Constant of One Pole

A useful approximate formula giving the decay time-constant9.4 $ \tau$ (in seconds) in terms of a pole radius $ R\in[0,1)$ is

$\displaystyle \zbox {\tau\approx \frac{T}{1-R}} \protect$ (9.8)

where $ T$ denotes the sampling interval in seconds, and we assume $ T\ll\tau$.

The exact relation between $ \tau$ and $ R$ is obtained by sampling an exponential decay:

$\displaystyle e^{-t/\tau} \;\rightarrow\; e^{-nT/\tau} \;\isdef \; R^n $

Thus, setting $ n=1$ yields

$\displaystyle R = e^{-T/\tau}. $

Expanding the right-hand side in a Taylor series and neglecting terms higher than first order gives

$\displaystyle e^{-\frac{T}{\tau}} = 1 - \frac{T}{\tau} + \frac{1}{2!}\left(\fr... ...{1}{3!}\left(\frac{T}{\tau}\right)^3 + \cdots \;\approx\; 1 - \frac{T}{\tau}, $

which derives $ R\approx 1-T/\tau$. Solving for $ \tau$ then gives Eq.$ \,$(8.8). From its derivation, we see that the approximation is valid for $ T\ll\tau$. Thus, as long as the impulse response of a pole $ p$ ``rings'' for many samples, the formula $ \tau\approx T/(1-\vert p\vert)$ should well estimate the time-constant of decay in seconds. The time-constant estimate in samples is of course $ 1/(1-\vert p\vert)$. For higher-order systems, the approximate decay time is $ 1/(1-R_{\mbox{max}})$, where $ R_{\mbox{max}}$ is the largest pole magnitude (closest to the unit circle) in the (stable) system.

Unstable Poles--Unit Circle Viewpoint

We saw in §8.4 that an LTI filter is stable if and only if all of its poles are strictly inside the unit circle ($ \vert z\vert=1$) in the complex $ z$ plane. In particular, a pole $ p$ outside the unit circle ($ \vert p\vert>1$) gives rise to an impulse-response component proportional to $ p^n$ which grows exponentially over time $ n$. We also saw in §6.2 that the z transform of a growing exponential does not not converge on the unit circle in the $ z$ plane. However, this was the case for a causal exponential $ u(n)p^n$, where $ u(n)$ is the unit-step function (which switches from 0 to 1 at time 0). If the same exponential is instead anticausal, i.e., of the form $ u(-n)p^n$, then, as we'll see in this section, its z transform does exist on the unit circle, and the pole is in exactly the same place as in the causal case. Therefore,to unambiguously invert a z transform, we must know its region of convergence. The critical question is whether the region of convergence includes the unit circle: If it does, then each pole outside the unit circle corresponds to an anticausal, finite energy, exponential, while each pole inside corresponds to the usual causal decaying exponential.

Geometric Series

The essence of the situation can be illustrated using a simple geometric series. Let $ R$ be any real (or complex) number. Then we have

$\displaystyle \frac{1}{1-R} \eqsp 1 + R + R^2 + R^3 + \cdots \quad < \infty$   when$\displaystyle \quad\vert R\vert<1. $

In other words, the geometric series $ 1 + R + R^2 + R^3 + \cdots$ is guaranteed to be summable when $ \vert R\vert<1$, and in that case, the sum is given by $ 1/(1-R)$. On the other hand, if $ \vert R\vert>1$, we can rewrite $ 1/(1-R)$ as $ -R^{-1}/(1-R^{-1})$ to obtain

$\displaystyle \frac{1}{1-R} \eqsp \frac{-R^{-1}}{1-R^{-1}} \eqsp -R^{-1}\left[1 + R^{-1} + R^{-2} + R^{-3} + \cdots \right] $

which is summable when $ \vert R\vert>1$. Thus, $ 1/(1-R)$ is a valid closed-form sum whether or not $ \vert R\vert$ is less than or greater than 1. When $ \vert R\vert<1$, it is the sum of the causal geometric series in powers of $ R$. When $ \vert R\vert>1$, it is the sum of the causal geometric series in powers of $ R^{-1}$, or, an anticausal geometric series in (negative) powers of $ R$.

One-Pole Transfer Functions

We can apply the same analysis to a one-pole transfer function. Let $ p\in{\bf C}$ denote any real or complex number:

$\displaystyle H(z) \eqsp \frac{1}{1-pz^{-1}} \eqsp 1 + pz^{-1}+ pz^{-2}+ pz^{-3} + \cdots $

The convergence criterion is now $ \vert pz^{-1}\vert<1$, or $ \vert z\vert>\vert p\vert$. For the region of convergence to include the unit circle (our frequency axis), we must have $ \vert p\vert<1$, which is our usual stability criterion for a pole at $ z=p$. The inverse z transform is then the causal decaying sampled exponential

$\displaystyle H(z) \;\longleftrightarrow\; h(n) = u(n)p^n $

Now consider the rewritten case:

\begin{eqnarray*} \frac{1}{1-pz^{-1}} &=& \frac{-p^{-1}z}{1-p^{-1}z} \\ &=& -p^... ...cdots\right]\\ &\leftrightarrow& - u(-n-1)p^n,\quad n\in{\bf Z} \end{eqnarray*}

where the inverse z transform is the inverse bilateral z transform. In this case, the convergence criterion is $ \vert p^{-1}z\vert<1$, or $ \vert z\vert<\vert p\vert$, and this region includes the unit circle when $ \vert p\vert>1$.

In summary, when the region-of-convergence of the z transform is assumed to include the unit circle of the $ z$ plane, poles inside the unit circle correspond to stable, causal, decaying exponentials, while poles outside the unit circle correspond to anticausal exponentials that decay toward time $ -\infty$, and stop before time zero.

Figure 8.8 illustrates the two types of exponentials (causal and anticausal) that correspond to poles (inside and outside the unit circle) when the z transform region of convergence is defined to include the unit circle.

myFourFiguresToWidthpolesout11polesout21polesout12polesout220.52Left column: Causal exponential decay, pole at $ p=0.9$. Right column: Anticausal exponential decay, pole at $ p=1/0.9$. Top: Pole-zero diagram. Bottom: Corresponding impulse response, assuming the region of convergence includes the unit circle in the $ z$ plane.

Poles and Zeros of the Cepstrum

The complex cepstrum of a sequence $ h(n)$ is typically defined as the inverse Fourier transform of its log spectrum [60]

$\displaystyle {\tilde h}(n)\isdef \frac{1}{2\pi}\int_{-\pi}^\pi \ln[H(e^{j\omega})] e^{j\omega n}d\omega, $

where $ \ln(x)$ denotes the natural logarithm (base $ e$) of $ x$, and $ H(e^{j\omega})$ denotes the Fourier transform (DTFT) of $ h(n)$,

$\displaystyle H(e^{j\omega})\isdef \sum_{n=-\infty}^\infty h(n)e^{-j\omega n}. $

(The real cepstrum, in contrast, is the inverse Fourier transform of the log-magnitude spectrum.) An equivalent definition (when the DTFT exists) is to define the complex cepstrum of $ h$ as the inverse z transform of $ \ln H(z)$. The cepstrum has numerous applications in digital signal processing including speech modeling [60] and pitch detection [34]. In §11.7, we use the cepstrum to compute minimum-phase spectra corresponding to a given spectral magnitude--an important tool in digital filter design.

From Eq.$ \,$(8.2), the log z transform can be written in terms of the factored form as

$\displaystyle \ln H(z)$ $\displaystyle =$ $\displaystyle \ln(g)$  
    $\displaystyle + \ln(1-q_1z^{-1}) + \ln(1-q_2z^{-1}) + \cdots + \ln (1-q_Mz^{-1})$  
    $\displaystyle - \ln(1-p_1z^{-1}) - \ln(1-p_2z^{-1}) - \cdots - \ln(1-p_Nz^{-1})$  
  $\displaystyle =$ $\displaystyle \ln(g) + \sum_{m=1}^M\ln(1-q_mz^{-1}) - \sum_{n=1}^N\ln(1-p_nz^{-1})$  
  $\displaystyle =$ $\displaystyle \ln(g) - \sum_{m=1}^M\ln\left(\frac{1}{1-q_mz^{-1}}\right) + \sum... ...eft(\frac{1}{1-p_nz^{-1}}\right),\qquad{} % force eqn no. to next line \protect$ (9.9)

where $ q_m$ denotes the $ m$th zero and $ p_n$ denotes the $ n$th pole of the z transform $ H(z)$. Applying the Maclaurin series expansion

$\displaystyle \ln\left(\frac{1}{1-x}\right) = x + \frac{x^2}{2} + \frac{x^3}{3} + \cdots + \frac{x^k}{k} + \cdots,\qquad \left\vert x\right\vert<1, $

we have that each pole and zero term in Eq.$ \,$(8.9) may be expanded as

\begin{eqnarray*} \ln\left(\frac{1}{1-q_iz^{-1}}\right) &=& \sum_{n=1}^\infty \f... ...^{-n}, \qquad \left\vert z\right\vert>\left\vert p_i\right\vert. \end{eqnarray*}

Since the region of convergence of the z transform must include the unit circle (where the spectrum (DTFT) is defined), we see that the Maclaurin expansion gives us the inverse z transform of all terms of Eq.$ \,$(8.9) corresponding to poles and zeros inside the unit circle of the $ z$ plane. Since the poles must be inside the unit circle anyway for stability, this restriction is normally not binding for the poles. However, zeros outside the unit circle--so-called ``non-minimum-phase zeros''--are used quite often in practice.

For a zero (or pole) outside the unit circle, we may rewrite the corresponding term of Eq.$ \,$(8.9) as

\begin{eqnarray*} \ln\left(\frac{1}{1-q_iz^{-1}}\right) &=& \ln\left(\frac{q_i^{... ...n,\qquad \left\vert z\right\vert<\left\vert q_i^{-1}\right\vert, \end{eqnarray*}

where we used the Maclaurin series expansion for $ \ln(1/(1-x))$ once again with the region of convergence including the unit circle. The infinite sum in this expansion is now the bilateral z transform of an anticausal sequence, as discussed in §8.7. That is, the time-domain sequence is zero for nonnegative times ($ n\geq 0$) and the sequence decays in the direction of time minus-infinity. The factored-out terms $ -z/q_i$ and $ -z/p_i$, for all poles and zeros outside the unit circle, can be collected together and associated with the overall gain factor $ g$ in Eq.$ \,$(8.9), resulting in a modified scaling and time-shift for the original sequence $ h(n)$ which can be dealt with separately [60].

When all poles and zeros are inside the unit circle, the complex cepstrum is causal and can be expressed simply in terms of the filter poles and zeros as

$\displaystyle {\tilde h}(n) = \left\{\begin{array}{ll} \ln(g), & n=0 \\ [5pt] \... ...ystyle\sum_{k=1}^M \frac{q_k^n}{n}, & n=1,2,3,\ldots\,, \\ \end{array}\right. $

where $ N$ is the number of poles and $ M$ is the number of zeros. Note that when $ N > M$, there are really $ N - M$ additional zeros at $ z=0$, but these contribute zero to the complex cepstrum (since $ q_i=0$). Similarly, when $ M>N$, there are $ M - N$ additional poles at $ z=0$ which also contribute zero (since $ p_i=0$).

In summary, each stable pole contributes a positive decaying exponential (weighted by $ 1/n$) to the complex cepstrum, while each zero inside the unit circle contributes a negative weighted-exponential of the same type. The decaying exponentials start at time 1 and have unit amplitude (ignoring the $ 1/n$ weighting) in the sense that extrapolating them to time 0 (without the $ 1/n$ weighting) would use the values $ p_i^0 = 1$ and $ -q_i^0 = -1$. The decay rates are faster when the poles and zeros are well inside the unit circle, but cannot decay slower than $ 1/n$.

On the other hand, poles and zeros outside the unit circle contribute anticausal exponentials to the complex cepstrum, negative for the poles and positive for the zeros.

Conversion to Minimum Phase

As discussed in §11.7, any spectrum can be converted to minimum-phase form (without affecting the spectral magnitude) by computing its cepstrum and replacing any anticausal components with corresponding causal components. In other words, the anticausal part of the cepstrum, if any, is ``flipped'' about time zero so that it adds to the causal part. Doing this corresponds to reflecting non-minimum phase zeros (and any unstable poles) inside the unit circle in a manner that preserves spectral magnitude. The original spectral phase is then replaced by the unique minimum phase corresponding to the given spectral magnitude.

A matlab listing for computing a minimum-phase spectrum from the magnitude spectrum is given in §J.11.

Hilbert Transform Relations

Closely related to the cepstrum are the so-called Hilbert transform relations that relate the real and imaginary parts of the spectra of causal signals. In particular, for minimum-phase spectra, the cepstrum is causal, and this means that the log-magnitude and phase form a Hilbert-transform pair. Methods for designing allpass filters have been based on this relationship (see §B.2.2). For more about cepstra and Hilbert transform relations, see [60].

Pole-Zero Analysis Problems

See http://ccrma.stanford.edu/~jos/filtersp/Pole_Zero_Analysis_Problems.html.

Next Section:
Implementation Structures for Recursive Digital Filters
Previous Section:
Frequency Response Analysis