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:
The transfer function is also called the system function [60].
Let 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
:
![$\displaystyle \zbox {H(z) = \frac{Y(z)}{X(z)}}
$](http://www.dsprelated.com/josimages_new/filters/img622.png)
![$ X(z)$](http://www.dsprelated.com/josimages_new/filters/img306.png)
![$ H(z)$](http://www.dsprelated.com/josimages_new/filters/img308.png)
![$ Y(z)$](http://www.dsprelated.com/josimages_new/filters/img307.png)
![$ H(z)$](http://www.dsprelated.com/josimages_new/filters/img308.png)
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 is
defined to be
where
![$ z$](http://www.dsprelated.com/josimages_new/filters/img45.png)
![$ n = 0$](http://www.dsprelated.com/josimages_new/filters/img133.png)
![$ -\infty$](http://www.dsprelated.com/josimages_new/filters/img582.png)
![]() |
(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 maps each integer
to a complex
number (
) or real number (
). The z transform
of
, on the other hand,
, maps every complex number
to a new complex number
. On a higher
level, the z transform, viewed as a linear operator, maps an entire
signal
to its z transform
. We think of this as a ``function to
function'' mapping. We may say
is the z transform of
by writing
![$\displaystyle \zbox {X \leftrightarrow x}
$](http://www.dsprelated.com/josimages_new/filters/img629.png)
![$\displaystyle X(z) = {\cal Z}_z\{x(\cdot)\}
$](http://www.dsprelated.com/josimages_new/filters/img630.png)
![$\displaystyle X = {\cal Z}\{x\}.
$](http://www.dsprelated.com/josimages_new/filters/img631.png)
![$ X(z)
\leftrightarrow x(n)$](http://www.dsprelated.com/josimages_new/filters/img632.png)
![$ n$](http://www.dsprelated.com/josimages_new/filters/img89.png)
![$ z$](http://www.dsprelated.com/josimages_new/filters/img45.png)
![$ n\in{\bf Z}$](http://www.dsprelated.com/josimages_new/filters/img387.png)
![$ z\in{\bf C}$](http://www.dsprelated.com/josimages_new/filters/img626.png)
The z transform of a signal can be regarded as a polynomial in
, 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.
$](http://www.dsprelated.com/josimages_new/filters/img633.png)
![$ X(z) = 1 + 2z^{-1}+ 3z^{-2} = 1 + 2z^{-1}+ 3(z^{-1})^2$](http://www.dsprelated.com/josimages_new/filters/img634.png)
Existence of the Z Transform
The z transform of a finite-amplitude
signal 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
, and finite real numbers
and
,
such that
for all
. The
bounding exponential may even be growing with
(
). These are
not the most general conditions for existence of the z transform, but they
suffice for most practical purposes.
For a signal growing as
, for
, one
would naturally expect the z transform
to be defined only in the
region
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
$](http://www.dsprelated.com/josimages_new/filters/img640.png)
![$ \left\vert z\right\vert>\exp(\sigma)$](http://www.dsprelated.com/josimages_new/filters/img639.png)
![$ \sigma<0\,\Leftrightarrow\,\exp(\sigma)<1$](http://www.dsprelated.com/josimages_new/filters/img641.png)
![$ z$](http://www.dsprelated.com/josimages_new/filters/img45.png)
![$ z$](http://www.dsprelated.com/josimages_new/filters/img45.png)
More generally, it turns out that, in all cases of practical interest,
the domain of can be extended to include the
entire complex plane, except at isolated ``singular''
points7.2 at which
approaches
infinity (such as at
when
).
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
such that
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
,
, with
for
.
Shift Theorem
The shift theorem says that a delay of samples
in the time domain corresponds to a multiplication by
in the frequency domain:
![$\displaystyle {\cal Z}_z\{$](http://www.dsprelated.com/josimages_new/filters/img649.png)
![$\displaystyle _\Delta\{x\}\} \;=\; z^{-\Delta} X(z), \; \Delta\ge 0,
$](http://www.dsprelated.com/josimages_new/filters/img650.png)
![$\displaystyle \zbox {x(n-\Delta) \;\leftrightarrow\; z^{-\Delta} X(z), \; \Delta\ge 0.}
$](http://www.dsprelated.com/josimages_new/filters/img651.png)
![$ x(\cdot - \Delta)$](http://www.dsprelated.com/josimages_new/filters/img652.png)
![$ x(\cdot)$](http://www.dsprelated.com/josimages_new/filters/img107.png)
![$ \Delta$](http://www.dsprelated.com/josimages_new/filters/img647.png)
![$ z^{-\Delta}X(z)$](http://www.dsprelated.com/josimages_new/filters/img653.png)
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*}](http://www.dsprelated.com/josimages_new/filters/img654.png)
where we used the causality assumption for
.
Convolution Theorem
The convolution theorem for z transforms states that for any (real or)
complex causal signals and
,
convolution in the time domain is multiplication in the
domain, i.e.,
![$\displaystyle \zbox {x\ast y \;\leftrightarrow\; X\cdot Y}
$](http://www.dsprelated.com/josimages_new/filters/img657.png)
![$\displaystyle {\cal Z}_z\{x \ast y\} \;=\; X(z)Y(z),
$](http://www.dsprelated.com/josimages_new/filters/img658.png)
![$ X(z)\isdef {\cal Z}_z(x)$](http://www.dsprelated.com/josimages_new/filters/img659.png)
![$ Y(z)\isdef {\cal Z}_z(y)$](http://www.dsprelated.com/josimages_new/filters/img660.png)
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*}](http://www.dsprelated.com/josimages_new/filters/img661.png)
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
from a linear
time-invariant filter with input
and impulse response
is given
by the convolution of
and
, i.e.,
where ``
![$ \ast $](http://www.dsprelated.com/josimages_new/filters/img568.png)
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
where H(z) is the z transform of the filter impulse response. We may divide Eq.
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
![$ X(z)$](http://www.dsprelated.com/josimages_new/filters/img306.png)
![$\displaystyle H(z) \eqsp \frac{Y(z)}{X(z)} \;\isdef \; \hbox{transfer function}.
$](http://www.dsprelated.com/josimages_new/filters/img664.png)
![$ h(n)$](http://www.dsprelated.com/josimages_new/filters/img536.png)
![$ H(z)=Y(z)/X(z)$](http://www.dsprelated.com/josimages_new/filters/img305.png)
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))
Let's take the z transform of both sides, denoting the transform by
. Because
is a linear operator,
it may be distributed through the terms on the right-hand side as
follows:7.3
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
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*}](http://www.dsprelated.com/josimages_new/filters/img668.png)
Factoring out the common terms and
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].
$](http://www.dsprelated.com/josimages_new/filters/img669.png)
![\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*}](http://www.dsprelated.com/josimages_new/filters/img670.png)
the z transform of the difference equation yields
![$\displaystyle A(z)\,Y(z) = B(z)\,X(z).
$](http://www.dsprelated.com/josimages_new/filters/img671.png)
![$ Y(z)/X(z)$](http://www.dsprelated.com/josimages_new/filters/img309.png)
![$ H(z)$](http://www.dsprelated.com/josimages_new/filters/img308.png)
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.
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
Factored Form
By the fundamental theorem of algebra, every th order polynomial
can be factored into a product of
first-order polynomials.
Therefore, Eq.
(6.5) above can be written in
factored form as
The numerator roots
![$ \{q_1,\ldots,q_M\}$](http://www.dsprelated.com/josimages_new/filters/img674.png)
![$ \{p_1, \ldots,
p_N\}$](http://www.dsprelated.com/josimages_new/filters/img675.png)
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:
- Transfer functions of filters in series multiply together.
- Transfer functions of filters in parallel sum together.
Series Case
![]() |
Figure 6.1 illustrates the series connection of two
filters
and
.
The output
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).
$](http://www.dsprelated.com/josimages_new/filters/img680.png)
![$ H_1(z)$](http://www.dsprelated.com/josimages_new/filters/img34.png)
![$ H_2(z)$](http://www.dsprelated.com/josimages_new/filters/img35.png)
![$ H(z)=H_1(z)H_2(z)$](http://www.dsprelated.com/josimages_new/filters/img36.png)
Parallel Case
Figure 6.2 illustrates the parallel combination of two
filters. The filters and
are driven by the
same input signal
, and their respective outputs
and
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).
$](http://www.dsprelated.com/josimages_new/filters/img683.png)
![$ {\cal Z}\{y_1+y_2\} = {\cal Z}\{y_1\}+{\cal Z}\{y_2\}$](http://www.dsprelated.com/josimages_new/filters/img684.png)
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),
$](http://www.dsprelated.com/josimages_new/filters/img685.png)
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
$](http://www.dsprelated.com/josimages_new/filters/img686.png)
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:
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*}](http://www.dsprelated.com/josimages_new/filters/img688.png)
and . (The case
is addressed in the next section.)
The denominator coefficients
are called the poles of the
transfer function, and each numerator
is called the
residue of pole
. Equation (6.7) is general only if the poles
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
into first-order terms,7.4e.g., using the roots function in matlab.
The residue
corresponding to pole
may be found
analytically as
when the poles
![$ p_i$](http://www.dsprelated.com/josimages_new/filters/img691.png)
Note that in Eq.(6.8), there is always a pole-zero cancellation at
. That is, the term
is always cancelled by an
identical term in the denominator of
, which must exist because
has a pole at
. The residue
is simply the
coefficient of the one-pole term
in the partial
fraction expansion of
at
. The transfer function
is
, in the limit, as
.
Example
Consider the two-pole filter
![$\displaystyle H(z) \eqsp \frac{1}{(1-z^{-1})(1-0.5z^{-1})}.
$](http://www.dsprelated.com/josimages_new/filters/img700.png)
![$ p_1=1$](http://www.dsprelated.com/josimages_new/filters/img701.png)
![$ p_2=0.5$](http://www.dsprelated.com/josimages_new/filters/img702.png)
![\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*}](http://www.dsprelated.com/josimages_new/filters/img703.png)
We thus conclude that
![$\displaystyle H(z) \eqsp \frac{2}{1-z^{-1}} - \frac{1}{1-0.5z^{-1}}.
$](http://www.dsprelated.com/josimages_new/filters/img704.png)
![$\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})}
$](http://www.dsprelated.com/josimages_new/filters/img705.png)
Complex Example
To illustrate an example involving complex poles, consider the filter
![$\displaystyle H(z) \eqsp \frac{g}{1+z^{-2}},
$](http://www.dsprelated.com/josimages_new/filters/img706.png)
![$ g$](http://www.dsprelated.com/josimages_new/filters/img435.png)
![$ g$](http://www.dsprelated.com/josimages_new/filters/img435.png)
![$ p_1=j$](http://www.dsprelated.com/josimages_new/filters/img707.png)
![$ p_2=-j$](http://www.dsprelated.com/josimages_new/filters/img708.png)
![$\displaystyle H(z) \eqsp \frac{g}{(1-jz^{-1})(1+jz^{-1})}.
$](http://www.dsprelated.com/josimages_new/filters/img709.png)
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
![\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*}](http://www.dsprelated.com/josimages_new/filters/img710.png)
Thus,
![$\displaystyle H(z) \eqsp \frac{g/2}{1-jz^{-1}} + \frac{g/2}{1+jz^{-1}}.
$](http://www.dsprelated.com/josimages_new/filters/img711.png)
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 and
are real (implying that
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
denote any one-pole
section in the PFE of Eq.
(6.7). Then if
is complex and
describes a real filter, we will also find
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*}](http://www.dsprelated.com/josimages_new/filters/img715.png)
Expressing the pole in polar form as
,
and the residue as
,
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}}.
$](http://www.dsprelated.com/josimages_new/filters/img718.png)
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 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 in Eq.
(6.7), consider the sum of two
first-order terms by direct calculation:
![]() |
(7.9) |
Notice that the numerator order, viewed as a polynomial in
![$ z^{-1}$](http://www.dsprelated.com/josimages_new/filters/img91.png)
![$ N$](http://www.dsprelated.com/josimages_new/filters/img278.png)
![$ r_i/(1-p_i z^{-1})$](http://www.dsprelated.com/josimages_new/filters/img698.png)
![$ M=N-1$](http://www.dsprelated.com/josimages_new/filters/img604.png)
![$ N$](http://www.dsprelated.com/josimages_new/filters/img278.png)
![$ M<N$](http://www.dsprelated.com/josimages_new/filters/img689.png)
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}}.
$](http://www.dsprelated.com/josimages_new/filters/img722.png)
![$\displaystyle h_i(n) \eqsp r_i p_i^n, \quad n=0,1,2,\ldots\,.
$](http://www.dsprelated.com/josimages_new/filters/img723.png)
![$ H(z)$](http://www.dsprelated.com/josimages_new/filters/img308.png)
![$\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\,.
$](http://www.dsprelated.com/josimages_new/filters/img724.png)
![$ h$](http://www.dsprelated.com/josimages_new/filters/img580.png)
![$ N$](http://www.dsprelated.com/josimages_new/filters/img278.png)
![$ i$](http://www.dsprelated.com/josimages_new/filters/img571.png)
![$ i$](http://www.dsprelated.com/josimages_new/filters/img571.png)
![$ p_i$](http://www.dsprelated.com/josimages_new/filters/img691.png)
![$ i$](http://www.dsprelated.com/josimages_new/filters/img571.png)
![$ i$](http://www.dsprelated.com/josimages_new/filters/img571.png)
![$ r_i$](http://www.dsprelated.com/josimages_new/filters/img692.png)
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].
$](http://www.dsprelated.com/josimages_new/filters/img725.png)
![$ z^{-1}$](http://www.dsprelated.com/josimages_new/filters/img91.png)
The case of repeated poles is addressed in §6.8.5 below.
FIR Part of a PFE
When in Eq.
(6.7), we may perform a step of long division
of
to produce an FIR part in parallel with a
strictly proper IIR part:
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*}](http://www.dsprelated.com/josimages_new/filters/img727.png)
When , we define
. 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
where
![$ K=M-N\geq 0$](http://www.dsprelated.com/josimages_new/filters/img730.png)
We may compare these two PFE alternatives as follows:
Let denote
,
, and
.
(I.e., we use a subscript to indicate polynomial order, and `
' is
omitted for notational simplicity.) Then for
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*}](http://www.dsprelated.com/josimages_new/filters/img735.png)
In the first form, the
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 (
) can be used to match additional
terms in the impulse response after the FIR part
has
``died out''.
In summary, an arbitrary digital filter transfer function with
distinct poles can always be expressed as a parallel combination
of complex one-pole filters, together with a parallel FIR part
when
. 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 (the so-called
biquad section) can be written when
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}}.
$](http://www.dsprelated.com/josimages_new/filters/img741.png)
![$ d=z^{-1}$](http://www.dsprelated.com/josimages_new/filters/img742.png)
![$ H(z)$](http://www.dsprelated.com/josimages_new/filters/img308.png)
![$ d$](http://www.dsprelated.com/josimages_new/filters/img743.png)
![$\displaystyle H(d^{-1}) \eqsp g\frac{b_2 d^2 + b_1 d + 1 }{a_2 d^2 + a_1 d + 1}
$](http://www.dsprelated.com/josimages_new/filters/img744.png)
![% 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}](http://www.dsprelated.com/josimages_new/filters/img745.png)
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}
$](http://www.dsprelated.com/josimages_new/filters/img746.png)
![$\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}}.
$](http://www.dsprelated.com/josimages_new/filters/img747.png)
The delayed form of the partial fraction expansion is obtained by
leaving the coefficients in their original order. This corresponds
to writing as a ratio of polynomials in
:
![$\displaystyle H(z) \eqsp g\frac{z^2 + b_1 z + b_2 }{z^2 + a_1 z + a_2}
$](http://www.dsprelated.com/josimages_new/filters/img748.png)
![% 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}](http://www.dsprelated.com/josimages_new/filters/img749.png)
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}}.
$](http://www.dsprelated.com/josimages_new/filters/img750.png)
Numerical examples of partial fraction expansions are given in §6.8.8
below. Another worked example, in which the filter
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 .
This gives a linear system of
equations in
unknowns
,
.
Yet another method for finding residues is by means of Taylor series
expansions of the numerator and denominator
about each
pole
, 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}}.
$](http://www.dsprelated.com/josimages_new/filters/img752.png)
![$\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}.
$](http://www.dsprelated.com/josimages_new/filters/img753.png)
![$\displaystyle \frac{r_{1,1}}{(1-pz^{-1})^2} + \frac{r_{1,2}}{(1-pz^{-1})}
$](http://www.dsprelated.com/josimages_new/filters/img754.png)
![$ p$](http://www.dsprelated.com/josimages_new/filters/img48.png)
Dealing with Repeated Poles Analytically
A pole of multiplicity has
residues associated with it. For example,
and the three residues associated with the pole
![$ z=1/2$](http://www.dsprelated.com/josimages_new/filters/img759.png)
Let denote the
th residue associated with the pole
,
.
Successively differentiating
times with
respect to
and setting
isolates the residue
:
![\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*}](http://www.dsprelated.com/josimages_new/filters/img765.png)
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}}
$](http://www.dsprelated.com/josimages_new/filters/img766.png)
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*}](http://www.dsprelated.com/josimages_new/filters/img767.png)
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\,.}
$](http://www.dsprelated.com/josimages_new/filters/img768.png)
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}\;.
$](http://www.dsprelated.com/josimages_new/filters/img769.png)
![]() |
![]() |
![]() |
|
![]() |
![]() |
||
![]() |
![]() |
||
![]() |
![]() |
||
![]() |
![]() |
||
![]() |
![]() |
(7.13) |
Note that
![$ n+1$](http://www.dsprelated.com/josimages_new/filters/img777.png)
![$ n$](http://www.dsprelated.com/josimages_new/filters/img89.png)
![$ n$](http://www.dsprelated.com/josimages_new/filters/img89.png)
![$ \vert p\vert<1$](http://www.dsprelated.com/josimages_new/filters/img778.png)
![$ n$](http://www.dsprelated.com/josimages_new/filters/img89.png)
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
and
:
The finite limits on the summation result from the fact that both
![$ h_1$](http://www.dsprelated.com/josimages_new/filters/img782.png)
![$ h_2$](http://www.dsprelated.com/josimages_new/filters/img503.png)
![$\displaystyle \sum_{m=0}^n r^m = \frac{1-r^{n+1}}{1-r}
$](http://www.dsprelated.com/josimages_new/filters/img783.png)
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
![$\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}.
$](http://www.dsprelated.com/josimages_new/filters/img784.png)
![$ p_1$](http://www.dsprelated.com/josimages_new/filters/img785.png)
![$ p_2$](http://www.dsprelated.com/josimages_new/filters/img786.png)
![$ \left\vert p_1\right\vert>\left\vert p_2\right\vert$](http://www.dsprelated.com/josimages_new/filters/img787.png)
![$ h(n)$](http://www.dsprelated.com/josimages_new/filters/img536.png)
![$ p_1^n$](http://www.dsprelated.com/josimages_new/filters/img788.png)
![$ n$](http://www.dsprelated.com/josimages_new/filters/img89.png)
![$ \left\vert p_2\right\vert>\left\vert p_1\right\vert$](http://www.dsprelated.com/josimages_new/filters/img789.png)
![$ p_2^n$](http://www.dsprelated.com/josimages_new/filters/img790.png)
Going back to Eq.(6.14), we have
![]() |
(7.15) |
Setting
![$ p_1=p_2=p$](http://www.dsprelated.com/josimages_new/filters/img792.png)
![]() |
(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 , 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]
$](http://www.dsprelated.com/josimages_new/filters/img795.png)
![$\displaystyle h_1(n)\eqsp \sum_{m=0}^n 1\cdot 1 \eqsp n+1
$](http://www.dsprelated.com/josimages_new/filters/img796.png)
![\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*}](http://www.dsprelated.com/josimages_new/filters/img797.png)
Alternate Stability Criterion
In §5.6 (page ), a filter was defined to
be stable if its impulse response
decays to 0 in
magnitude as time
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
, where
is some
finite-order polynomial in
, and
is the
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
plane, i.e., when
. 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,
is
irreducible, while
is reducible, since
there is the common factor of
in the numerator and
denominator. Using this terminology, we may state the following
stability criterion:
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}}
$](http://www.dsprelated.com/josimages_new/filters/img805.png)
![]() |
(7.17) |
for
![$ M<N$](http://www.dsprelated.com/josimages_new/filters/img689.png)
for
![$ M\geq N$](http://www.dsprelated.com/josimages_new/filters/img690.png)
![$ z^{-(K+1)}$](http://www.dsprelated.com/josimages_new/filters/img807.png)
- When
, perform a step of long division to obtain an FIR part
and a strictly proper IIR part
.
- Find the
poles
,
(roots of
).
- If the poles are distinct, find the
residues
,
from
- If there are repeated poles, find the additional residues via
the method of §6.8.5, and the general form of the PFE is
wheredenotes the number of distinct poles, and
denotes the multiplicity of the
th pole.
In step 2, the poles are typically found by factoring the
denominator polynomial . This is a dangerous step numerically
which may fail when there are many poles, especially when many poles
are clustered close together in the
plane.
The following matlab code illustrates factoring
to
obtain the three roots,
,
:
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}}
$](http://www.dsprelated.com/josimages_new/filters/img817.png)
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
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
![]() |
(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.
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 -8Thus, this example finds that
![$ H(z)$](http://www.dsprelated.com/josimages_new/filters/img308.png)
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
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 -8That this must work can be seen by looking at Eq.
![$ \,$](http://www.dsprelated.com/josimages_new/filters/img94.png)
![$ n=2$](http://www.dsprelated.com/josimages_new/filters/img821.png)
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