Free Books
   Free Books    Introduction to Digital Filters  

Transfer Function Analysis

This chapter discusses filter transfer functions and associated analysis. The transfer function provides an algebraic representation of a linear, time-invariant (LTI) filter in the frequency domain:
$\textstyle \parbox{0.8\textwidth}{The \emph{transfer function}\index{transfer f... ..., and $X(z)$\ denotes the {\it z} transform of the filter input signal $x(n)$.}$

The transfer function is also called the system function [60]. Let $ h(n)$ denote the impulse response of the filter. It turns out (as we will show) that the transfer function is equal to the z transform of the impulse response $ h(n)$:

$\displaystyle \zbox {H(z) = \frac{Y(z)}{X(z)}} $

Since multiplying the input transform $ X(z)$ by the transfer function $ H(z)$ gives the output transform $ Y(z)$, we see that $ H(z)$ embodies the transfer characteristics of the filter--hence the name. It remains to define ``z transform'', and to prove that the z transform of the impulse response always gives the transfer function, which we will do by proving the convolution theorem for z transforms.

The Z Transform

The bilateral z transform of the discrete-time signal $ x(n)$ is defined to be

$\displaystyle X(z) \isdefs \sum_{n=-\infty}^\infty x(n) z^{-n} \qquad\hbox{(bilateral {\it z} transform)} \protect$ (7.1)

where $ z$ is a complex variable. Since signals are typically defined to begin (become nonzero) at time $ n = 0$, and since filters are often assumed to be causal,7.1 the lower summation limit given above may be written as 0 rather than $ -\infty$ to yield the unilateral z transform:

$\displaystyle X(z) \isdefs \sum_{n=0}^\infty x(n) z^{-n} \qquad\hbox{(unilateral {\it z} transform)}$ (7.2)

The unilateral z transform is most commonly used. For inverting z transforms, see §6.8. Recall (§4.1) that the mathematical representation of a discrete-time signal $ x(n)$ maps each integer $ n\in{\bf Z}$ to a complex number ( $ x(n)\in{\bf C}$) or real number ( $ x(n)\in{\bf R}$). The z transform of $ x$, on the other hand, $ X(z)$, maps every complex number $ z\in{\bf C}$ to a new complex number $ X(z)\in{\bf C}$. On a higher level, the z transform, viewed as a linear operator, maps an entire signal $ x$ to its z transform $ X$. We think of this as a ``function to function'' mapping. We may say $ X$ is the z transform of $ x$ by writing

$\displaystyle \zbox {X \leftrightarrow x} $

or, using operator notation,

$\displaystyle X(z) = {\cal Z}_z\{x(\cdot)\} $

which can be abbreviated as

$\displaystyle X = {\cal Z}\{x\}. $

One also sees the convenient but possibly misleading notation $ X(z) \leftrightarrow x(n)$, in which $ n$ and $ z$ must be understood as standing for the entire domains $ n\in{\bf Z}$ and $ z\in{\bf C}$, as opposed to denoting particular fixed values. The z transform of a signal $ x$ can be regarded as a polynomial in $ z^{-1}$, with coefficients given by the signal samples. For example, the signal

$\displaystyle x(n) = \left\{\begin{array}{ll} n+1, & 0\leq n \leq 2 \\ [5pt] 0, & \mbox{otherwise} \\ \end{array}\right. $

has the z transform $ X(z) = 1 + 2z^{-1}+ 3z^{-2} = 1 + 2z^{-1}+ 3(z^{-1})^2$.

Existence of the Z Transform

The z transform of a finite-amplitude signal $ x$ will always exist provided (1) the signal starts at a finite time and (2) it is asymptotically exponentially bounded, i.e., there exists a finite integer $ n_f$, and finite real numbers $ A\geq 0$ and $ \sigma$, such that $ \left\vert x(n)\right\vert<A\exp(\sigma n)$ for all $ n\geq n_f$. The bounding exponential may even be growing with $ n$ ($ \sigma>0$). These are not the most general conditions for existence of the z transform, but they suffice for most practical purposes. For a signal $ x(n)$ growing as $ \exp(\sigma n)$, for $ \sigma>0$, one would naturally expect the z transform $ X(z)$ to be defined only in the region $ \left\vert z\right\vert>\exp(\sigma)$ of the complex plane. This is expected because the infinite series

$\displaystyle \sum_{n=0}^\infty e^{\sigma n} z^{-n} = \sum_{n=0}^\infty \left(\frac{e^{\sigma}}{z}\right)^n $

requires $ \left\vert z\right\vert>\exp(\sigma)$ to ensure convergence. Since $ \sigma<0\,\Leftrightarrow\,\exp(\sigma)<1$ for a decaying exponential, we see that the region of convergence of the $ z$ transform of a decaying exponential always includes the unit circle of the $ z$ plane. More generally, it turns out that, in all cases of practical interest, the domain of $ X(z)$ can be extended to include the entire complex plane, except at isolated ``singular'' points7.2 at which $ \vert X(z)\vert$ approaches infinity (such as at $ z=\exp(\sigma)$ when $ x(n)=\exp(\sigma n)$). The mathematical technique for doing this is called analytic continuation, and it is described in §D.1 as applied to the Laplace transform (the continuous-time counterpart of the z transform). A point to note, however, is that in the extension region (all points $ z$ such that $ \left\vert z\right\vert<\exp(\sigma)$ in the above example), the signal component corresponding to each singularity inside the extension region is ``flipped'' in the time domain. That is, ``causal'' exponentials become ``anticausal'' exponentials, as discussed in §8.7. The z transform is discussed more fully elsewhere [52,60], and we will derive below only what we will need.


Shift and Convolution Theorems

In this section, we prove the highly useful shift theorem and convolution theorem for unilateral z transforms. We consider the space of infinitely long, causal, complex sequences $ x(n)\in{\bf C}$, $ n\in{\bf Z}$, with $ x(n)=0$ for $ n < 0$.

Shift Theorem

The shift theorem says that a delay of $ \Delta$ samples in the time domain corresponds to a multiplication by $ z^{-\Delta}$ in the frequency domain:

$\displaystyle {\cal Z}_z\{$SHIFT$\displaystyle _\Delta\{x\}\} \;=\; z^{-\Delta} X(z), \; \Delta\ge 0, $

or, using more common notation,

$\displaystyle \zbox {x(n-\Delta) \;\leftrightarrow\; z^{-\Delta} X(z), \; \Delta\ge 0.} $

Thus, $ x(\cdot - \Delta)$, which is the waveform $ x(\cdot)$ delayed by $ \Delta$ samples, has the z transform $ z^{-\Delta}X(z)$.
Proof:
\begin{eqnarray*} {\cal Z}_z\{\mbox{{\sc Shift}}_\Delta\{x\}\} &\isdef & \sum_{n... ...ty}x(m) z^{-m} \\ &\isdef & z^{-\Delta} X(z), % \quad\pfendmath \end{eqnarray*}
where we used the causality assumption $ x(m)=0$ for $ m<0$.

Convolution Theorem

The convolution theorem for z transforms states that for any (real or) complex causal signals $ x$ and $ y$, convolution in the time domain is multiplication in the $ z$ domain, i.e.,

$\displaystyle \zbox {x\ast y \;\leftrightarrow\; X\cdot Y} $

or, using operator notation,

$\displaystyle {\cal Z}_z\{x \ast y\} \;=\; X(z)Y(z), $

where $ X(z)\isdef {\cal Z}_z(x)$, and $ Y(z)\isdef {\cal Z}_z(y)$. (See [84] for a development of the convolution theorem for discrete Fourier transforms.)
Proof:
\begin{eqnarray*} {\cal Z}_z(x\ast y) &\isdef & \sum_{n=0}^{\infty}(x\ast y)_n z... ...(by the Shift Theorem)}\\ &\isdef & X(z)Y(z) % \quad\pfendmath \end{eqnarray*}
The convolution theorem provides a major cornerstone of linear systems theory. It implies, for example, that any stable causal LTI filter (recursive or nonrecursive) can be implemented by convolving the input signal with the impulse response of the filter, as shown in the next section.

Z Transform of Convolution

From Eq.$ \,$(5.5), we have that the output $ y$ from a linear time-invariant filter with input $ x$ and impulse response $ h$ is given by the convolution of $ h$ and $ x$, i.e.,

$\displaystyle y(n) \eqsp (h \ast x)(n) \protect$ (7.3)

where ``$ \ast $'' means convolution as before. Taking the z transform of both sides of Eq.$ \,$(6.3) and applying the convolution theorem from the preceding section gives

$\displaystyle Y(z) \eqsp H(z)X(z) \protect$ (7.4)

where H(z) is the z transform of the filter impulse response. We may divide Eq.$ \,$(6.4) by $ X(z)$ to obtain

$\displaystyle H(z) \eqsp \frac{Y(z)}{X(z)} \;\isdef \; \hbox{transfer function}. $

This shows that, as a direct result of the convolution theorem, the z transform of an impulse response $ h(n)$ is equal to the transfer function $ H(z)=Y(z)/X(z)$ of the filter, provided the filter is linear and time invariant.

Z Transform of Difference Equations

Since z transforming the convolution representation for digital filters was so fruitful, let's apply it now to the general difference equation, Eq.$ \,$(5.1). To do this requires two properties of the z transform, linearity (easy to show) and the shift theorem (derived in §6.3 above). Using these two properties, we can write down the z transform of any difference equation by inspection, as we now show. In §6.8.2, we'll show how to invert by inspection as well. Repeating the general difference equation for LTI filters, we have (from Eq.$ \,$(5.1)) \begin{eqnarrayda} y(n) &=& b_0 x(n) &+& b_1 x(n - 1) + \cdots + b_M x(n - M)\\ & & &-& a_1 y(n - 1) - \cdots - a_N y(n - N).\\ \end{eqnarrayda} Let's take the z transform of both sides, denoting the transform by $ {\cal Z}\{\}$. Because $ {\cal Z}\{\}$ is a linear operator, it may be distributed through the terms on the right-hand side as follows:7.3\begin{eqnarrayda} {\cal Z}_z\{y(\cdot)\} &=& {\cal Z}\{ b_0 x(n) &+& b_1 x(n ... ...-M} X(z)\\ & & &-& a_1 z^{-1}Y(z) - \cdots - a_N z^{-N} Y(z), \end{eqnarrayda} where we used the superposition and scaling properties of linearity given on page [*], followed by use of the shift theorem, in that order. The terms in $ Y(z)$ may be grouped together on the left-hand side to get
\begin{eqnarray*} \lefteqn{Y(z) + a_1 z^{-1}Y(z) + \cdots + a_N z^{-N} Y(z) = }\... ... \\ & & b_0 X(z) + b_1 z^{-1}X(z) + \cdots + b_M z^{-M} X(z). \end{eqnarray*}
Factoring out the common terms $ Y(z)$ and $ X(z)$ gives

$\displaystyle Y(z)\left[1 + a_1 z^{-1}+ \cdots + a_N z^{-N}\right] = X(z)\left[b_0 + b_1 z^{-1}+ \cdots + b_M z^{-M}\right]. $

Defining the polynomials
\begin{eqnarray*} A(z) &\isdef & 1 + a_1\,z^{-1} + \,\cdots\, + a_N\,z^{-N}\\ B(z) &\isdef & b_0 + b_1\,z^{-1}+\,\cdots\,+b_M\,z^{-M}, \end{eqnarray*}
the z transform of the difference equation yields

$\displaystyle A(z)\,Y(z) = B(z)\,X(z). $

Finally, solving for $ Y(z)/X(z)$, which is by definition the transfer function $ H(z)$, gives

$\displaystyle H(z) \isdefs \frac{Y(z)}{X(z)} \eqsp \frac{b_0 + b_1 z^{-1}+ \cdo... ...^{-M}}{1 + a_1 z^{-1}+ \cdots + a_N z^{-N}} \isdefs \frac{B(z)}{A(z)}. \protect$ (7.5)

Thus, taking the z transform of the general difference equation led to a new formula for the transfer function in terms of the difference equation coefficients. (Now the minus signs for the feedback coefficients in the difference equation Eq.$ \,$(5.1) are explained.)

Factored Form

By the fundamental theorem of algebra, every $ N$th order polynomial can be factored into a product of $ N$ first-order polynomials. Therefore, Eq.$ \,$(6.5) above can be written in factored form as

$\displaystyle H(z) = b_0\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$ (7.6)

The numerator roots $ \{q_1,\ldots,q_M\}$ are called the zeros of the transfer function, and the denominator roots $ \{p_1, \ldots, p_N\}$ are called the poles of the filter. Poles and zeros are discussed further in Chapter 8.

Series and Parallel Transfer Functions

The transfer function conveniently captures the algebraic structure of a filtering operation with respect to series or parallel combination. Specifically, we have the following cases:
  1. Transfer functions of filters in series multiply together.
  2. Transfer functions of filters in parallel sum together.

Series Case

Figure 6.1: Series combination of transfer functions $ H_1(z)$ and $ H_2(z)$ to produce the combined transfer function $ H(z)=H_1(z)H_2(z)$.
\begin{figure}\input fig/series.pstex_t \end{figure}
Figure 6.1 illustrates the series connection of two filters $ H_1(z)=V(z)/X(z)$ and $ H_2(z)=Y(z)/V(z)$. The output $ v(n)$ from filter 1 is used as the input to filter 2. Therefore, the overall transfer function is

$\displaystyle H(z) \isdefs \frac{Y(z)}{X(z)} \eqsp \frac{H_2(z)V(z)}{X(z)} \eqsp H_2(z)H_1(z). $

In summary, if the output of filter $ H_1(z)$ is given as input to filter $ H_2(z)$ (a series combination), as shown in Fig.6.1, the overall transfer function is $ H(z)=H_1(z)H_2(z)$--transfer functions of filters connected in series multiply together.

Parallel Case

Figure 6.2: Parallel combination of transfer functions $ H_1(z)$ and $ H_2(z)$, yielding $ H(z)=H_1(z)+H_2(z)$.
\begin{figure}\input fig/parallel.pstex_t \end{figure}
Figure 6.2 illustrates the parallel combination of two filters. The filters $ H_1(z)$ and $ H_2(z)$ are driven by the same input signal $ x(n)$, and their respective outputs $ y_1(n)$ and $ y_2(n)$ are summed. The transfer function of the parallel combination is therefore

$\displaystyle H(z) \isdefs \frac{Y(z)}{X(z)} \eqsp \frac{Y_1(z) + Y_2(z)}{X(z)} \eqsp \frac{Y_1(z)}{X(z)} + \frac{Y_2(z)}{X(z)} \isdefs H_1(z)+H_2(z). $

where we needed only linearity of the z transform to have that $ {\cal Z}\{y_1+y_2\} = {\cal Z}\{y_1\}+{\cal Z}\{y_2\}$.

Series Combination is Commutative

Since multiplication of complex numbers is commutative, we have

$\displaystyle H_1(z)H_2(z)=H_2(z)H_1(z), $

which implies that any ordering of filters in series results in the same overall transfer function. Note, however, that the numerical performance of the overall filter is usually affected by the ordering of filter stages in a series combination [103]. Chapter 9 further considers numerical performance of filter implementation structures. By the convolution theorem for z transforms, commutativity of a product of transfer functions implies that convolution is commutative:

$\displaystyle h_1 \ast h_2 \;\leftrightarrow\; H_1\cdot H_2 \;=\; H_2\cdot H_1 \;\leftrightarrow\; h_2 \ast h_1 $

Partial Fraction Expansion

An important tool for inverting the z transform and converting among digital filter implementation structures is the partial fraction expansion (PFE). The term ``partial fraction expansion'' refers to the expansion of a rational transfer function into a sum of first and/or second-order terms. The case of first-order terms is the simplest and most fundamental:

$\displaystyle H(z) \isdefs \frac{B(z)}{A(z)} \eqsp \sum_{i=1}^{N} \frac{r_i}{1-p_iz^{-1}} \protect$ (7.7)

where
\begin{eqnarray*} B(z) &=& b_0 + b_1 z^{-1}+ b_2z^{-2}+ \cdots + b_M z^{-M}\\ A(z) &=& 1 + a_1 z^{-1}+ a_2z^{-2}+ \cdots + a_N z^{-N} \end{eqnarray*}
and $ M<N$. (The case $ M\geq N$ is addressed in the next section.) The denominator coefficients $ p_i$ are called the poles of the transfer function, and each numerator $ r_i$ is called the residue of pole $ p_i$. Equation (6.7) is general only if the poles $ p_i$ are distinct. (Repeated poles are addressed in §6.8.5 below.) Both the poles and their residues may be complex. The poles may be found by factoring the polynomial $ A(z)$ into first-order terms,7.4e.g., using the roots function in matlab. The residue $ r_i$ corresponding to pole $ p_i$ may be found analytically as

$\displaystyle r_i \eqsp \left.(1-p_iz^{-1})H(z)\right\vert _{z=p_i} \protect$ (7.8)

when the poles $ p_i$ are distinct. The matlab function residuez7.5 will find poles and residues computationally, given the difference-equation (transfer-function) coefficients. Note that in Eq.$ \,$(6.8), there is always a pole-zero cancellation at $ z=p_i$. That is, the term $ (1-p_iz^{-1})$ is always cancelled by an identical term in the denominator of $ H(z)$, which must exist because $ H(z)$ has a pole at $ z=p_i$. The residue $ r_i$ is simply the coefficient of the one-pole term $ 1/(1-p_i z^{-1})$ in the partial fraction expansion of $ H(z)$ at $ z=p_i$. The transfer function is $ r_i/(1-p_i z^{-1})$, in the limit, as $ z\to p_i$.

Example

Consider the two-pole filter

$\displaystyle H(z) \eqsp \frac{1}{(1-z^{-1})(1-0.5z^{-1})}. $

The poles are $ p_1=1$ and $ p_2=0.5$. The corresponding residues are then
\begin{eqnarray*} r_1 &=& \left.(1-z^{-1})H(z)\right\vert _{z=1} \eqsp \left.\f... ...\ &=& \left.\frac{1}{1-z^{-1}}\right\vert _{z=0.5} \eqsp -1\,. \end{eqnarray*}
We thus conclude that

$\displaystyle H(z) \eqsp \frac{2}{1-z^{-1}} - \frac{1}{1-0.5z^{-1}}. $

As a check, we can add the two one-pole terms above to get

$\displaystyle \frac{2}{1-z^{-1}} - \frac{1}{1-0.5z^{-1}} \eqsp \frac{2-z^{-1}- 1 + z^{-1}}{(1-z^{-1})(1-0.5z^{-1})} \eqsp \frac{1}{(1-z^{-1})(1-0.5z^{-1})} $

as expected.

Complex Example

To illustrate an example involving complex poles, consider the filter

$\displaystyle H(z) \eqsp \frac{g}{1+z^{-2}}, $

where $ g$ can be any real or complex value. (When $ g$ is real, the filter as a whole is real also.) The poles are then $ p_1=j$ and $ p_2=-j$ (or vice versa), and the factored form can be written as

$\displaystyle H(z) \eqsp \frac{g}{(1-jz^{-1})(1+jz^{-1})}. $

Using Eq.$ \,$(6.8), the residues are found to be
\begin{eqnarray*} r_1 &=& \left.(1-jz^{-1})H(z)\right\vert _{z=j} \eqsp \left.\... ...eft.\frac{g}{1-jz^{-1}}\right\vert _{z=-j} \eqsp \frac{g}{2}\,. \end{eqnarray*}
Thus,

$\displaystyle H(z) \eqsp \frac{g/2}{1-jz^{-1}} + \frac{g/2}{1+jz^{-1}}. $

A more elaborate example of a partial fraction expansion into complex one-pole sections is given in §3.12.1.

PFE to Real, Second-Order Sections

When all coefficients of $ A(z)$ and $ B(z)$ are real (implying that $ H(z)=B(z)/A(z)$ is the transfer function of a real filter), it will always happen that the complex one-pole filters will occur in complex conjugate pairs. Let $ (r,p)$ denote any one-pole section in the PFE of Eq.$ \,$(6.7). Then if $ p$ is complex and $ H(z)$ describes a real filter, we will also find $ (\overline{r},\pc)$ somewhere among the terms in the one-pole expansion. These two terms can be paired to form a real second-order section as follows:
\begin{eqnarray*} H(z) &=& \frac{r}{1-pz^{-1}} + \frac{\overline{r}}{1-\pc z^{-1... ...box{re}\left\{p\right\}z^{-1}+ \left\vert p\right\vert^2 z^{-2}} \end{eqnarray*}
Expressing the pole $ p$ in polar form as $ p=Re^{j\theta}$, and the residue as $ r=Ge^{j\phi}$, the last expression above can be rewritten as

$\displaystyle H(z) \eqsp 2G\frac{\cos(\phi)-\cos(\phi-\theta)z^{-1}}{1-2R\,\cos(\theta)z^{-1}+ R^2 z^{-2}}. $

The use of polar-form coefficients is discussed further in the section on two-pole filters (§B.1.3). Expanding a transfer function into a sum of second-order terms with real coefficients gives us the filter coefficients for a parallel bank of real second-order filter sections. (Of course, each real pole can be implemented in its own real one-pole section in parallel with the other sections.) In view of the foregoing, we may conclude that every real filter with $ M<N$ can be implemented as a parallel bank of biquads.7.6 However, the full generality of a biquad section (two poles and two zeros) is not needed because the PFE requires only one zero per second-order term. To see why we must stipulate $ M<N$ in Eq.$ \,$(6.7), consider the sum of two first-order terms by direct calculation:

$\displaystyle H_2(z) \eqsp \frac{r_1}{1-p_1z^{-1}} + \frac{r_2}{1-p_2z^{-1}} \eqsp \frac{(r_1 + r_2) - (r_1 p_2 + r_2 p_1) z^{-1}}{(1-p_1z^{-1})(1-p_2z^{-1})}$ (7.9)

Notice that the numerator order, viewed as a polynomial in $ z^{-1}$, is one less than the denominator order. In the same way, it is easily shown by mathematical induction that the sum of $ N$ one-pole terms $ r_i/(1-p_i z^{-1})$ can produce a numerator order of at most $ M=N-1$ (while the denominator order is $ N$ if there are no pole-zero cancellations). Following terminology used for analog filters, we call the case $ M<N$ a strictly proper transfer function.7.7 Thus, every strictly proper transfer function (with distinct poles) can be implemented using a parallel bank of two-pole, one-zero filter sections.

Inverting the Z Transform

The partial fraction expansion (PFE) provides a simple means for inverting the z transform of rational transfer functions. The PFE provides a sum of first-order terms of the form

$\displaystyle H_i(z) \eqsp \frac{r_i}{1-p_iz^{-1}}. $

It is easily verified that such a term is the z transform of

$\displaystyle h_i(n) \eqsp r_i p_i^n, \quad n=0,1,2,\ldots\,. $

Thus, the inverse z transform of $ H(z)$ is simply

$\displaystyle h(n) \eqsp \sum_{i=1}^N h_i(n) \eqsp \sum_{i=1}^N r_i p_i^n, \quad n=0,1,2,\ldots\,. $

Thus, the impulse response of every strictly proper LTI filter (with distinct poles) can be interpreted as a linear combination of sampled complex exponentials. Recall that a uniformly sampled exponential is the same thing as a geometric sequence. Thus, $ h$ is a linear combination of $ N$ geometric sequences. The term ratio of the $ i$th geometric sequence is the $ i$th pole, $ p_i$, and the coefficient of the $ i$th sequence is the $ i$th residue, $ r_i$. In the improper case, discussed in the next section, we additionally obtain an FIR part in the z transform to be inverted:

$\displaystyle F(z) \eqsp f_0 + f_1z^{-1}+ f_2z^{-2}+ \cdots + f_K z^{-K} \;\longleftrightarrow\; [f_0,f_1,\ldots,f_K,0,0,\ldots]. $

The FIR part (a finite-order polynomial in $ z^{-1}$) is also easily inverted by inspection. The case of repeated poles is addressed in §6.8.5 below.

FIR Part of a PFE

When $ M\geq N$ in Eq.$ \,$(6.7), we may perform a step of long division of $ B(z)/A(z)$ to produce an FIR part in parallel with a strictly proper IIR part:

$\displaystyle H(z) \isdefs \frac{B(z)}{A(z)} \eqsp F(z) + \sum_{i=1}^{N} \frac{r_i}{1-p_iz^{-1}} \protect$ (7.10)

where
\begin{eqnarray*} B(z) &=& b_0 + b_1 z^{-1}+ b_2z^{-2}+ \cdots + b_M z^{-M}\\ A... ...=& f_0 + f_1z^{-1}+ f_2z^{-2}+ \cdots + f_K z^{-K}, \quad K=M-N. \end{eqnarray*}
When $ M<N$, we define $ F(z)=0$. This type of decomposition is computed by the residuez function (a matlab function for computing a complete partial fraction expansion, as illustrated in §6.8.8 below). An alternate FIR part is obtained by performing long division on the reversed polynomial coefficients to obtain

$\displaystyle H(z) \eqsp F(z) + z^{-(K+1)}\sum_{i=1}^{N} \frac{r_i}{1-p_iz^{-1}}, \protect$ (7.11)

where $ K=M-N\geq 0$ is again the order of the FIR part. This type of decomposition is computed (as part of the PFE) by residued, described in §J.6 and illustrated numerically in §6.8.8 below. We may compare these two PFE alternatives as follows: Let $ A_N$ denote $ A(z)$, $ F_K\isdeftext F(z)$, and $ B_M\isdeftext B(z)$. (I.e., we use a subscript to indicate polynomial order, and `$ (z)$' is omitted for notational simplicity.) Then for $ K=M-N\geq 0$ we have two cases:
\begin{eqnarray*} (1) && H(z) \eqsp F_K + \frac{B^\prime_{N-1}}{A_N} \eqsp \frac... ..._N} \eqsp \frac{F_K A_N + z^{-(K+1)}B^{\prime\prime}_{N-1}}{A_N} \end{eqnarray*}
In the first form, the $ B^\prime_{N-1}$ coefficients are ``left justified'' in the reconstructed numerator, while in the second form they are ``right justified''. The second form is generally more efficient for modeling purposes, since the numerator of the IIR part ( $ B^{\prime\prime}_{N-1}(z)$) can be used to match additional terms in the impulse response after the FIR part $ F_K(z)$ has ``died out''. In summary, an arbitrary digital filter transfer function $ H(z)$ with $ N$ distinct poles can always be expressed as a parallel combination of complex one-pole filters, together with a parallel FIR part when $ M\geq N$. When there is an FIR part, the strictly proper IIR part may be delayed such that its impulse response begins where that of the FIR part leaves off. In artificial reverberation applications, the FIR part may correspond to the early reflections, while the IIR part provides the late reverb, which is typically dense, smooth, and exponentially decaying [86]. The predelay (``pre-delay'') control in some commercial reverberators is the amount of pure delay at the beginning of the reverberator's impulse response. Thus, neglecting the early reflections, the order of the FIR part can be viewed as the amount of predelay for the IIR part.

Example: The General Biquad PFE

The general second-order case with $ M=N=2$ (the so-called biquad section) can be written when $ b_0\ne 0$ as

$\displaystyle H(z) \eqsp g\frac{1 + b_1 z^{-1}+ b_2 z^{-2}}{1 + a_1 z^{-1}+ a_2 z^{-2}}. $

To perform a partial fraction expansion, we need to extract an order 0 (length 1) FIR part via long division. Let $ d=z^{-1}$ and rewrite $ H(z)$ as a ratio of polynomials in $ d$:

$\displaystyle H(d^{-1}) \eqsp g\frac{b_2 d^2 + b_1 d + 1 }{a_2 d^2 + a_1 d + 1} $

Then long division gives % For typesetting long division --- NEEDED WITHIN THE MAKEIMAGE ENV? % (raw TeX,... ...} & & b_1-\frac{b_2}{a_2}a_1 & 1-\frac{b_2}{a_2} & \end{array}\end{displaymath}
yielding

$\displaystyle H(d^{-1}) \eqsp g\frac{b_2}{a_2} + g\frac{\left(b_1-\frac{b_2}{a_2}a_1\right)d+ \left(1-\frac{b_2}{a_2}\right)}{a_2d^2 + a_1d + 1} $

or

$\displaystyle H(z) \eqsp g\frac{b_2}{a_2} + g\frac{\left(1-\frac{b_2}{a_2}\right) +\left(b_1-\frac{b_2}{a_2}a_1\right)z^{-1}}{1 + a_1z^{-1}+ a_2z^{-2}}. $

The delayed form of the partial fraction expansion is obtained by leaving the coefficients in their original order. This corresponds to writing $ H(z)$ as a ratio of polynomials in $ z$:

$\displaystyle H(z) \eqsp g\frac{z^2 + b_1 z + b_2 }{z^2 + a_1 z + a_2} $

Long division now looks like % For typesetting long division --- NEEDED WITHIN THE MAKEIMAGE ENV?\begin{dis... ...rule width 22\digitwidth}} & & b_1-a_1 & b_2-a_2 & \end{array}\end{displaymath}
giving

$\displaystyle H(z) \eqsp g + z^{-1}g\frac{(b_1-a_1) + (b_2-a_2)z^{-1}}{1 + a_1 z^{-1}+ a_2 z^{-2}}. $

Numerical examples of partial fraction expansions are given in §6.8.8 below. Another worked example, in which the filter $ y(n) = x(n) + 0.5^3 x(n-3) - 0.9^5 y(n-5)$ is converted to a set of parallel, second-order sections is given in §3.12. See also §9.2 regarding conversion to second-order sections in general, and §G.9.1 (especially Eq.$ \,$(G.22)) regarding a state-space approach to partial fraction expansion.

Alternate PFE Methods

Another method for finding the pole residues is to write down the general form of the PFE, obtain a common denominator, expand the numerator terms to obtain a single polynomial, and equate like powers of $ z^{-1}$. This gives a linear system of $ N$ equations in $ N$ unknowns $ r_i$, $ i=1,\ldots,N$. Yet another method for finding residues is by means of Taylor series expansions of the numerator $ B(z)$ and denominator $ A(z)$ about each pole $ p_i$, using l'Hôpital's rule.. Finally, one can alternatively construct a state space realization of a strictly proper transfer function (using, e.g., tf2ss in matlab) and then diagonalize it via a similarity transformation. (See Appendix G for an introduction to state-space models and diagonalizing them via similarity transformations.) The transfer function of the diagonalized state-space model is trivially obtained as a sum of one-pole terms--i.e., the PFE. In other words, diagonalizing a state-space filter realization implicitly performs a partial fraction expansion of the filter's transfer function. When the poles are distinct, the state-space model can be diagonalized; when there are repeated poles, it can be block-diagonalized instead, as discussed further in §G.10.

Repeated Poles

When poles are repeated, an interesting new phenomenon emerges. To see what's going on, let's consider two identical poles arranged in parallel and in series. In the parallel case, we have

$\displaystyle H_1(z) \eqsp \frac{r_1}{1-pz^{-1}} + \frac{r_2}{1-pz^{-1}} \eqsp \frac{r_1+r_2}{1-pz^{-1}} \isdefs \frac{r_3}{1-pz^{-1}}. $

In the series case, we get

$\displaystyle H_2(z) \eqsp \frac{r_1}{1-pz^{-1}} \cdot \frac{r_2}{1-pz^{-1}} \eqsp \frac{r_1r_2}{(1-pz^{-1})^2} \isdefs \frac{r_3}{(1-pz^{-1})^2}. $

Thus, two one-pole filters in parallel are equivalent to a new one-pole filter7.8 (when the poles are identical), while the same two filters in series give a two-pole filter with a repeated pole. To accommodate both possibilities, the general partial fraction expansion must include the terms

$\displaystyle \frac{r_{1,1}}{(1-pz^{-1})^2} + \frac{r_{1,2}}{(1-pz^{-1})} $

for a pole $ p$ having multiplicity 2.

Dealing with Repeated Poles Analytically

A pole of multiplicity $ m_i$ has $ m_i$ residues associated with it. For example,
$\displaystyle H(z)$ $\displaystyle \isdef$ $\displaystyle \frac{7 - 5z^{-1}+ z^{-2}}{\left(1-\frac{1}{2}z^{-1}\right)^3}$  
  $\displaystyle =$ $\displaystyle \frac{1}{\left(1-\frac{1}{2}z^{-1}\right)^3} + \frac{2}{\left(1-\frac{1}{2}z^{-1}\right)^2} + \frac{4}{\left(1-\frac{1}{2}z^{-1}\right)} \protect$ (7.12)

and the three residues associated with the pole $ z=1/2$ are 1, 2, and 4. Let $ r_{ij}$ denote the $ j$th residue associated with the pole $ p_i$, $ j=1,\ldots,m_i$. Successively differentiating $ (1-p_iz^{-1})^{m_i}H(z)$ $ k-1$ times with respect to $ z^{-1}$ and setting $ z=p_i$ isolates the residue $ r_{ik}$:
\begin{eqnarray*} r_{i1} &=& \left.(1-p_iz^{-1})^{m_i}H(z)\right\vert _{z=p_i}\\... ...ac{d^3}{d(z^{-1})^3} (1-p_iz^{-1})^{m_i}H(z)\right\vert _{z=p_i} \end{eqnarray*}
or

$\displaystyle \zbox {r_{ik} = \left.\frac{1}{(k-1)!(-p_i)^{k-1}}\frac{d^{k-1}}{d(z^{-1})^{k-1}} (1-p_iz^{-1})^{m_i}H(z)\right\vert _{z=p_i}} $

Example

For the example of Eq.$ \,$(6.12), we obtain
\begin{eqnarray*} r_{11} &=& \left.\left(1-\frac{1}{2}z^{-1}\right)^3H(z)\right\... ...}{dz^{-1}} (-5 + 2z^{-1})\right\vert _{z^{-1}=2} = 2\cdot 2 = 4. \end{eqnarray*}

Impulse Response of Repeated Poles

In the time domain, repeated poles give rise to polynomial amplitude envelopes on the decaying exponentials corresponding to the (stable) poles. For example, in the case of a single pole repeated twice, we have

$\displaystyle \zbox {\frac{1}{\left(1-pz^{-1}\right)^2} \;\longleftrightarrow\; (n+1) p^n, \quad n=0,1,2,\ldots\,.} $


Proof: First note that

$\displaystyle \frac{d}{dz^{-1}}\left(\frac{1}{1-pz^{-1}}\right) = (-1)(1-pz^{-1})^{-2}(-p) = \frac{p}{\left(1-pz^{-1}\right)^2}\;. $

Therefore,
$\displaystyle \frac{1}{\left(1-pz^{-1}\right)^2}$ $\displaystyle =$ $\displaystyle \frac{1}{p}\, \frac{d}{dz^{-1}}\left(\frac{1}{1-pz^{-1}}\right)$  
  $\displaystyle =$ $\displaystyle \frac{1}{p}\, \frac{d}{dz^{-1}} \left(1 + pz^{-1}+ p^2z^{-2}+ p^3 z^{-3} + \cdots \right)$  
  $\displaystyle =$ $\displaystyle \frac{1}{p} \left(0 + p + 2p^2z^{-1}+ 3p^3z^{-2}+ \cdots \right)$  
  $\displaystyle =$ $\displaystyle 1 + 2pz^{-1}+ 3p^2z^{-2}+ \cdots$  
  $\displaystyle =$ $\displaystyle \sum_{n=0}^{\infty}(n+1)p^n z^{-n}$  
  $\displaystyle \isdef$ $\displaystyle {\cal Z}\left\{(n+1)p^n\right\} \;\longleftrightarrow\; (n+1)p^n.$ (7.13)

Note that $ n+1$ is a first-order polynomial in $ n$. Similarly, a pole repeated three times corresponds to an impulse-response component that is an exponential decay multiplied by a quadratic polynomial in $ n$, and so on. As long as $ \vert p\vert<1$, the impulse response will eventually decay to zero, because exponential decay always overtakes polynomial growth in the limit as $ n$ goes to infinity.

So What's Up with Repeated Poles?

In the previous section, we found that repeated poles give rise to polynomial amplitude-envelopes multiplying the exponential decay due to the pole. On the other hand, two different poles can only yield a convolution (or sum) of two different exponential decays, with no polynomial envelope allowed. This is true no matter how closely the poles come together; the polynomial envelope can occur only when the poles merge exactly. This might violate one's intuitive expectation of a continuous change when passing from two closely spaced poles to a repeated pole. To study this phenomenon further, consider the convolution of two one-pole impulse-responses $ h_1(n) = p_1^n$ and $ h_2(n) = p_2^n$:

$\displaystyle h(n) \isdef (h_1\ast h_2)(n) = \sum_{m=0}^n h_1(m)h_2(n-m) = \sum... ...^n p_1^{m}p_2^{n-m} = p_2^n\sum_{m=0}^n \left(\frac{p_1}{p_2}\right)^m \protect$ (7.14)

The finite limits on the summation result from the fact that both $ h_1$ and $ h_2$ are causal. Recall the closed-form sum of a truncated geometric series:

$\displaystyle \sum_{m=0}^n r^m = \frac{1-r^{n+1}}{1-r} $

Applying this to Eq.$ \,$(6.14) yields

$\displaystyle h(n) = p_2^n \frac{1-(p_1/p_2)^{n+1}}{1-(p_1/p_2)} = \frac{p_2^{n+1}-p_1^{n+1}}{p_2-p_1} = \frac{p_1^{n+1}-p_2^{n+1}}{p_1-p_2}. $

Note that the result is symmetric in $ p_1$ and $ p_2$. If $ \left\vert p_1\right\vert>\left\vert p_2\right\vert$, then $ h(n)$ becomes proportional to $ p_1^n$ for large $ n$, while if $ \left\vert p_2\right\vert>\left\vert p_1\right\vert$, it becomes instead proportional to $ p_2^n$. Going back to Eq.$ \,$(6.14), we have

$\displaystyle h(n) = p_2^n\sum_{m=0}^n \left(\frac{p_1}{p_2}\right)^m = p_1^n\sum_{m=0}^n \left(\frac{p_2}{p_1}\right)^m.$ (7.15)

Setting $ p_1=p_2=p$ yields

$\displaystyle h(n) = (n+1)p^n$ (7.16)

which is the first-order polynomial amplitude-envelope case for a repeated pole. We can see that the transition from ``two convolved exponentials'' to ``single exponential with a polynomial amplitude envelope'' is perfectly continuous, as we would expect. We also see that the polynomial amplitude-envelopes fundamentally arise from iterated convolutions. This corresponds to the repeated poles being arranged in series, rather than in parallel. The simplest case is when the repeated pole is at $ p=1$, in which case its impulse response is a constant:

$\displaystyle \frac{1}{1-z^{-1}} \eqsp 1 + z^{-1}+ z^{-2}+ \cdots \;\longleftrightarrow\; [1,1,1,\ldots] $

The convolution of a constant with itself is a ramp:

$\displaystyle h_1(n)\eqsp \sum_{m=0}^n 1\cdot 1 \eqsp n+1 $

The convolution of a constant and a ramp is a quadratic, and so on:7.9
\begin{eqnarray*} h_2(n)&=&\sum_{m=0}^n (m+1)\cdot 1 \eqsp \frac{(n+1)(n+2)}{2}\... ...+1)(m+2)}{2}\cdot 1\eqsp \frac{(n+1)(n+2)(n+3)}{3!}\\ &\cdots& \end{eqnarray*}

Alternate Stability Criterion

In §5.6 (page [*]), a filter was defined to be stable if its impulse response $ h(n)$ decays to 0 in magnitude as time $ n$ goes to infinity. In §6.8.5, we saw that the impulse response of every finite-order LTI filter can be expressed as a possible FIR part (which is always stable) plus a linear combination of terms of the form $ a_i(n)p_i^n$, where $ a_i(n)$ is some finite-order polynomial in $ n$, and $ p_i$ is the $ i$th pole of the filter. In this form, it is clear that the impulse response always decays to zero when each pole is strictly inside the unit circle of the $ z$ plane, i.e., when $ \vert p_i\vert<1$. Thus, having all poles strictly inside the unit circle is a sufficient criterion for filter stability. If the filter is observable (meaning that there are no pole-zero cancellations in the transfer function from input to output), then this is also a necessary criterion. A transfer function with no pole-zero cancellations is said to be irreducible. For example, $ H(z) = (1+z^{-1})/(1-z^{-1})$ is irreducible, while $ H(z) = (1-z^{-2})/(1-2z^{-2}+z^{-2})$ is reducible, since there is the common factor of $ (1-z^{-1})$ in the numerator and denominator. Using this terminology, we may state the following stability criterion:
$\textstyle \parbox{0.8\textwidth}{\emph{An irreducible transfer function $H(z)$\ is stable if and only if its poles have magnitude less than one.}}$
This characterization of stability is pursued further in §8.4, and yet another stability test (most often used in practice) is given in §8.4.1.

Summary of the Partial Fraction Expansion

In summary, the partial fraction expansion can be used to expand any rational z transform

$\displaystyle H(z) \eqsp \frac{B(z)}{A(z)} \eqsp \frac{b_0 + b_1 z^{-1}+ \cdots + b_M z^{-M}}{1 + a_1 z^{-1}+ \cdots + a_N z^{-N}} $

as a sum of first-order terms

$\displaystyle H(z) \isdefs \frac{B(z)}{A(z)} \eqsp \sum_{i=1}^{N} \frac{r_i}{1-p_iz^{-1}} \protect$ (7.17)

for $ M<N$, and

$\displaystyle H(z) \eqsp F(z) + z^{-(K+1)}\sum_{i=1}^{N}\frac{r_i}{1-p_iz^{-1}} \protect$ (7.18)

for $ M\geq N$, where the term $ z^{-(K+1)}$ is optional, but often preferred. For real filters, the complex one-pole terms may be paired up to obtain second-order terms with real coefficients. The PFE procedure occurs in two or three steps:
  1. When $ M\geq N$, perform a step of long division to obtain an FIR part $ F(z)$ and a strictly proper IIR part $ B^\prime(z)/A(z)$.
  2. Find the $ N$ poles $ p_i$, $ i=1,\ldots,N$ (roots of $ A(z)$).
  3. If the poles are distinct, find the $ N$ residues $ r_i$, $ i=1,\ldots,N$ from

    $\displaystyle r_i = \left.(1-p_iz^{-1})\frac{B(z)}{A(z)}\right\vert _{z=p_i} $

  4. If there are repeated poles, find the additional residues via the method of §6.8.5, and the general form of the PFE is

    $\displaystyle H(z) \eqsp F(z) + z^{-(K+1)}\sum_{i=1}^{N_p}\sum_{k=1}^{m_i}\frac{r_{i,k}}{(1-p_iz^{-1})^k} \protect$ (7.19)

    where $ N_p$ denotes the number of distinct poles, and $ m_i\ge 1$ denotes the multiplicity of the $ i$th pole.
In step 2, the poles are typically found by factoring the denominator polynomial $ A(z)$. This is a dangerous step numerically which may fail when there are many poles, especially when many poles are clustered close together in the $ z$ plane. The following matlab code illustrates factoring $ A(z) = 1 - z^{-3}$ to obtain the three roots, $ p_k=e^{jk2\pi/3}$, $ k=0,1,2$:
A = [1 0 0 -1];  % Filter denominator polynomial
poles = roots(A) % Filter poles
See Chapter 9 for additional discussion regarding digital filters implemented as parallel sections (especially §9.2.2).

Software for Partial Fraction Expansion

Figure 6.3 illustrates the use of residuez (§J.5) for performing a partial fraction expansion on the transfer function

$\displaystyle H(z) \eqsp \frac{1 + 0.5^3 z^{-3}}{1 + 0.9^5z^{-5}} $

The complex-conjugate terms can be combined to obtain two real second-order sections, giving a total of one real first-order section in parallel with two real second-order sections, as discussed and depicted in §3.12.
Figure 6.3: Use of residuez to perform a partial fraction expansion of an IIR filter transfer function $ H(z)=B(z)/A(z)$.
  
B = [1 0 0 0.125];
A = [1 0 0 0 0 0.9^5];
[r,p,f] = residuez(B,A)
% r =
%   0.16571
%   0.22774 - 0.02016i
%   0.22774 + 0.02016i
%   0.18940 + 0.03262i
%   0.18940 - 0.03262i
%
% p =
%   -0.90000
%   -0.27812 - 0.85595i
%   -0.27812 + 0.85595i
%    0.72812 - 0.52901i
%    0.72812 + 0.52901i
%
% f = [](0x0)

Example 2

For the filter
$\displaystyle H(z)$ $\displaystyle \isdef$ $\displaystyle \frac{2+6z^{-1}+6z^{-2}+2z^{-3}}{1-2z^{-1}+z^{-2}}$ (7.20)
  $\displaystyle =$ $\displaystyle (2+10z^{-1}) + z^{-2}\left[\frac{8}{1-z^{-1}} + \frac{16}{(1-z^{-1})^2}\right] \protect$ (7.21)

we obtain the output of residued (§J.6) shown in Fig.6.4. In contrast to residuez, residued delays the IIR part until after the FIR part. In contrast to this result, residuez returns r=[-24;16] and f=[10;2], corresponding to the PFE

$\displaystyle H(z) = 10+2z^{-1}-\frac{24}{1-z^{-1}} + \frac{16}{(1-z^{-1})^2},$ (7.22)

in which the FIR and IIR parts have overlapping impulse responses. See Sections J.5 and J.6 starting on page [*] for listings of residuez, residued and related discussion.
Figure 6.4: Use of residued to perform a partial fraction expansion of an IIR filter transfer function $ H(z)=B(z)/A(z)$.
  
B=[2 6 6 2]; A=[1 -2 1];
[r,p,f,m] = residued(B,A)
% r =
%    8
%   16
%
% p =
%   1
%   1
%
% f =
%    2  10
%
% m =
%   1
%   2

Polynomial Multiplication in Matlab

The matlab function conv (convolution) can be used to perform polynomial multiplication. For example:
B1 = [1 1];   % 1st row of Pascal's triangle
B2 = [1 2 1]; % 2nd row of Pascal's triangle
B3 = conv(B1,B2) % 3rd row
% B3 = 1  3  3  1
B4 = conv(B1,B3) % 4th row
% B4 = 1  4  6  4  1
% ...
The matlab conv(B1,B2) is identical to filter(B1,1,B2), except that conv returns the complete convolution of its two input vectors, while filter truncates the result to the length of the ``input signal'' B2.7.10 Thus, if B2 is zero-padded with length(B1)-1 zeros, it will return the complete convolution:
B1 = [1 2 3];
B2 = [4 5 6 7];
conv(B1,B2)
% ans = 4  13  28  34  32  21
filter(B1,1,B2)
% ans = 4  13  28  34
filter(B1,1,[B2,zeros(1,length(B1)-1)])
% ans = 4  13  28  34  32  21

Polynomial Division in Matlab

The matlab function deconv (deconvolution) can be used to perform polynomial long division in order to split an improper transfer function into its FIR and strictly proper parts:
B = [ 2 6 6 2]; % 2*(1+1/z)^3
A = [ 1 -2 1];  % (1-1/z)^2
[firpart,remainder] = deconv(B,A)
% firpart =
%   2  10
% remainder =
%    0    0   24   -8
Thus, this example finds that $ H(z)$ is as written in Eq.$ \,$(6.21). This result can be checked by obtaining a common denominator in order to recalculate the direct-form numerator:
Bh = remainder + conv(firpart,A)
%  = 2 6 6 2
The operation deconv(B,A) can be implemented using filter in a manner analogous to the polynomial multiplication case (see §6.8.8 above):
firpart = filter(B,A,[1,zeros(1,length(B)-length(A))])
%       = 2 10
remainder = B - conv(firpart,A)
%         =  0 0 24 -8
That this must work can be seen by looking at Eq.$ \,$(6.21) and noting that the impulse-response of the remainder (the strictly proper part) does not begin until time $ n=2$, so that the first two samples of the impulse-response come only from the FIR part. In summary, we may conveniently use convolution and deconvolution to perform polynomial multiplication and division, respectively, such as when converting transfer functions to various alternate forms. When carrying out a partial fraction expansion on a transfer function having a numerator order which equals or exceeds the denominator order, a necessary preliminary step is to perform long division to obtain an FIR filter in parallel with a strictly proper transfer function. This section describes how an FIR part of any length can be extracted from an IIR filter, and this can be used for PFEs as well as for more advanced applications [].

Problems

See http://ccrma.stanford.edu/~jos/filtersp/Transfer_Function_Analysis_Problems.html.
Next Section:
Frequency Response Analysis
Previous Section:
Time Domain Digital Filter Representations