Chapter 7

IIR Digital Filter Design

We have focused our discussion so far on FIR or nonrecursive filter design. There are, however, some disadvantages in using nonrecursive filters:

(i) Get ripples. Usually need long filters to suppress them.

(ii) Must have smooth transitions since it is a sum of smooth functions. Difficult to get sharp transitions.

(iii) For real-time operation, it is necessary to introduce a delay approximately equal to half of the filter length (see P.S. #6). If the filter length is long, this may not be acceptable.

In this chapter, we study the use of recursive filters of the form

y_{n} = ∑_{k=1}^{N} a_{k}y_{n−k} + ∑_{k=0}^{M} b_{k}u_{n−k}

(Note: this is causal since responds to past values. Noncausal recursive filters are possible but almost never used).

Since frequency response is

H̃(e^{i2πf}) = \frac{∑_{k=0}^{M} b_{k}e^{−2πikf}}{1 − ∑_{k=1}^{N} a_{k}e^{−2πikf}}

can get narrow transitions as denominator comes near zero. Also recursive filters generally have infinite memory, and may thus need fewer coefficients.

Fourier approach is clearly not applicable. How do we go about determining a_{k}, b_{k}? We proceed indirectly and take advantage of known results from analog filter design. It is known that the following class of filters has properties that allow them to give good l.p. filters.

7.1 Butterworth filters

H_{a}(ω)H_{a}(−ω) = \frac{1}{1 + (\frac{ω}{ω_{c}} )^{2N}}

ω_{c} \frac{1}{\sqrt{2}}  3 db pt |H_{a}(ω)|

Filter parameters are ω_{c} and N. Butterworth filters are also called maximally-flat filters as their amplitude characteristics have the largest number of derivatives with respect to ω vanishing at 0.

Suppose we want to design an analog l.p. Butterworth filter

\frac{1}{A^{2}} 1 \frac{1}{1 + ϵ^{2}}  ω_{p}  ω_{s}

Then

\frac{1}{1 + ε^{2}} = \frac{1}{1 + (\frac{ω_{p}}{ω_{c}} )^{2N}} ⇒ (\frac{ω_{p}}{ω_{c}} )^{2N} = ε^{2}

and

\frac{1}{A^{2}} = \frac{1}{1 + (\frac{ω_{s}}{ω_{c}} )^{2N}} ⇒ (\frac{ω_{s}}{ω_{c}} )^{2N} = A^{2} − 1

⇒ (\frac{ω_{p}}{ω_{s}} )^{2N} = \frac{ε^{2}}{A^{2} − 1}

Solving this equation for N will generally not yield an integer. Since N is an integer, we use the ceil function where ceil(x) is the smallest integer ≥ x.

N = ceil (\begin{matrix}log & \frac{ε}{\sqrt{A^{2} − 1}} \\ log & \frac{ω_{p}}{ω_{s}} \end{matrix} )

For the critical frequency ω_{c}, there are two choices:

ω_{c} = \frac{ω_{p}}{(ε^{1/N})}

or

ω_{c} = \frac{ω_{s}}{(A^{2} − 1)^{1/2N}}

As shown in P.S. #8, the first choice gives a filter which meets the passband specs exactly, and exceeds the stopband specs, which the second choice gives a filter which meets the stopband specs exactly, and exceeds the passband specs.

Once N and ω_{c} are found, we have determined H_{a}(iω)H_{a}(−iω). Since the analog filter has transfer function H(s), we still need to obtain H(s) from H_{a}(iω)H_{a}(−iω) to complete the analog filter design.

Example:

ω_{p} = 10^{4}
ω_{s} = 10^{5}

ε^{2} = \frac{1}{999} ⇒ \frac{1}{1 + ε^{2}} = \frac{999}{1000}

A^{2} = 10^{3}

Since N is an integer, we take

N = ceil (\frac{1}{log 10^{−1}} log ( \frac{\frac{1}{999}}{999} )^{1/2}) = ceil(log 999) = 3

Two choices for the critical frequency ω_{c}:

ω_{c1} = \frac{ω_{p}}{(ε^{2})^{1/2N}}
= \frac{10^{4}}{(0.999)^{\frac{1}{6}}} = 1.00017 × 10^{4}

or

ω_{c2} = \frac{ω_{s}}{(A^{2} − 1)^{1/2N}}
= \frac{10^{5}}{(999)^{\frac{1}{6}}} = 3.1628 × 10^{4}

As mentioned before, using ω_{c1} gives a filter whose passband specs are exactly met and stopband specs exceeded, while using ω_{c2} gives a filter whose stopband specs are exactly met and passband specs exceeded.

7.2 Factorization

First we substitute s for iω to get H(s)H(−s) from H(iω)H(−iω). For fixed ω_{c} and N, we then have

H(s)H(−s) = \frac{1}{1 + (\frac{s}{iω_{c}} )^{2N}}

The factorization problem is to factor the above product to get a stable H(s).

First, let us note the following property of the distribution of poles of H(s)H(−s): If s = σ + iω is a pole, then −s = −(σ + iω) is also a pole. Furthermore, since H(s) is a rational function with real coefficients, s^{∗} = σ − iω is also a pole. But then −s^{∗} = −σ + iω is also a pole. This shows that complex poles of H(s)H(−s) with nonzero real and imaginary parts must occur in a 4-pole pattern symmetric with respect to the origin. Real poles occur in pairs symmetric with respect to the imaginary axis while purely imaginary poles occur in pairs symmetric with respect to the real axis.

Pole Patterns for H(s)H(−s)

Poles of H(s)H(−s) are at

s = (−1)^{1/2N}(iω_{c})

= ω_{c}e^{i( \frac{π}{2} + \frac{(2k + 1)}{2N}π)} k = 0, …, 2N − 1

(i) No pole on imaginary axis, since this would require \frac{(2k + 1)}{N} to be even (so that \frac{2k + 1}{2N} is an integer), impossible.

(ii) One on real axis for N odd (corresponds to k = \frac{N − 1}{2}, which is an integer for N odd), but not for

N even (\frac{2k + 1}{N} cannot be an integer if N even).

For example N = 3, 1^{st} pole at e\frac{π}{2} i + i\frac{π}{6} . It is easy to see that the poles corresponding to k = 0, 1, …, N − 1 are poles located in the left half plane. In fact, once the poles in the second quadrant are determined, the rest of the stable poles can be determined through symmetry.

\frac{π}{6}  ω_{c}

For this example, the stable poles are ω_{c}(−0.5 + i0.866), −ω_{c}, and ω_{c}(−0.5 − i0.866).

Once the poles of H(s) have been found, the complete transfer function can be written as

H(s) = ∏_{j=1}^{N} \frac{1}{1 + \frac{s}{s_{j}}}

where −s_{j}, j = 1, …, N are the N stable poles of H(s). Later we shall describe some shortcuts in writing down the complete transfer function H(s) using the symmetry properties of the poles.

Now that the analog filter has been found to satisfy the lowpass filter specs, the next task is to translate the analog design into a digital filter which satisfies the digital lowpass filter specs. As in analog filters, we require digital IIR filters to be stable.

Stability of recursive filters:

Equivalent to stability of difference equation

y_{n} = ∑_{k=1}^{N} a_{k}y_{n−k} + ∑_{k=0}^{M} b_{k}u_{n−k}

Readily seen that if the impulse response is h_{k} then system stable (BIBO) ⇔

∑_{k=0}^{∞} |h_{k}| < ∞

⇔ H(z) has no poles on or outside the unit circle.

Now recursive filter has transfer function

H(z) = \frac{z^{N−M}(b_{0}z^{M} + b_{1}z^{M−1} + … + b_{M})}{z^{N} − a_{1}z^{N−1} − … − a_{N}}

Hence it is stable ⇔ the polynomial z^{N} − a_{1}z^{N−1} − … − a_{N} has no roots on or outside the unit circle.

Since the analog Butterworth filter design gives a stable analog filter, it is important to ensure that going from analog to digital gives stable recursive digital filters.

7.3 Getting recursive digital filters from analog designs

Consider the mapping

z = \frac{1 + s}{1 − s} or s = \frac{z − 1}{z + 1}

This mapping is called the bilinear transformation.

For

s = σ_{a} + iω_{a}
z = \frac{1 + σ_{a} + iω_{a}}{1 − σ_{a} − iω_{a}}

If

σ_{a} < 0 |z| < 1 σ_{a} > 0 ⇒ |z| > 1
σ_{a} = 0 |z| = 1

This means that the open left half plane for s is mapped onto the open unit disk for z, the imaginary axis for s is mapped onto the unit circle for z, and the open right half plane for s is mapped onto the exterior of the unit circle for z. This means that the bilinear transformation preserves stability.

What happens to frequencies? If we denote the digital frequency by ω

iω_{a} = \frac{e^{iω} − 1}{e^{iω} + 1} = \frac{i sin \frac{ω}{2}}{cos \frac{ω}{2}}
⇒ ω_{a} = tan \frac{ω}{2} or 2 tan^{−1} ω_{a} = ω

πω  ω = 2 tan^{−1} ω_{a} ω_{a} ω^{as} ω^{ap}

The compression of the analog frequency ω_{a} due to the bilinear transformation is called frequency warping. Knowing that we are going to use the bilinear transformation allows us to pre-compensate for the effects of frequency warping. This procedure is called frequency prewarping.

Given digital frequency specs with passband frequency ω_{p} and stopband frequency ω_{s}, translate these specs into those for the analog filter using

\begin{matrix}ω_{s} & to & ω_{as} = tan \frac{ω_{s}}{2} \\ ω_{p} & to & ω_{ap} = tan \frac{ω_{p}}{2} \end{matrix} \} prewarping

By attaining specs for H_{a} (tan \frac{ω_{s}}{2} ) etc. for the analog filter, the bilinear transformation will then yield a

digital filter that will attain the specs for H_{d}(e^{iω}).

Example: Design a digital IIR filter satisfying the following specs:

0.89125 ≤ |H(e^{iω})| ≤ 1 0 ≤ ω ≤ 0.2π  (7.1)
|H(e^{iω})| ≤ 0.17783 0.3π ≤ ω ≤ π  (7.2)

Step 1. Prewarp the digital filter specs to the following analog filter specs: ω_{ap} = tan \frac{0.2π}{2} = 0.3249, and ω_{as} = tan \frac{0.3π}{2} = 0.5095. The analog filter specs corresponding to 7.1 and 7.2 become

0.89125 ≤ |H(iω_{a})| ≤ 1 0 ≤ ω_{a} ≤ 0.3249  (7.3)
|H(iω_{a})| ≤ 0.17783 0.5095 ≤ ω_{a} < ∞  (7.4)

Step 2. Determine the analog Butterworth filter parameters N and ω_{c}:

\frac{1}{1 + (\frac{ω_{ap}}{ω_{c}} )^{2N}} = 0.89125^{2} = 0.7943

\frac{1}{1 + (\frac{ω_{as}}{ω_{c}} )^{2N}} = 0.17783^{2} = 0.03162

Combining these equations gives

(\frac{ω_{as}}{ω_{ap}} )^{2N} = \frac{\frac{1}{0.03162} − 1}{\frac{1}{0.7943} − 1} = 118.265

Therefore

N = ceil (\frac{log(118.265)}{2 log(\frac{0.5095}{0.3249})}) = ceil(5.30428) = 6

Of the 2 choices for ω_{c}, for this example we take

ω_{c} = \frac{0.5095}{\frac{1}{0.17783^{2}} − 1} = \frac{0.5095}{1.329925} = 0.383115

Step 3. Factorize H(s)H(−s): The poles in the second quadrant are located at these angles on the circle of

radius ω_{c}

(\frac{π}{2} + \frac{π}{12}), (\frac{π}{2} + \frac{3π}{12} ), (\frac{π}{2} + \frac{5π}{12} )

Their real parts are given by −0.09916, −0.2709, and −0.37, respectively. They are paired with their complex conjugates to make up the 6 poles needed for H(s).

Step 4. Now note that for every complex number s_{1} lying on the circle of radius ω_{c}, we can express the product

\frac{1}{1 + \frac{s}{s_{1}}} \frac{1}{1 + \frac{s}{s_{1}^{∗}}} = \frac{|s_{1}|^{2}}{s^{2} + 2 Re(s_{1})s + |s_{1}|^{2}}

= \frac{ω_{c}^{2}}{s^{2} + 2 Re(s_{1})s + ω_{c}^{2}}

where Re denotes the real part. Thus knowing the real part of the complex conjugate pair of poles will enable us to write down the quadratic factor immediately. In the above product, the poles are identified as −s_{1} and −s_{1}^{∗}. Thus, if −s_{1} is in the second quadrant, Re(s_{1}) > 0. We can now write down the complete analog filter transfer function

H(s) = \frac{ω_{c}^{6}}{(s^{2} + 0.198315s + ω_{c}^{2})(s^{2} + 0.54181s + ω_{c}^{2})(s^{2} + 0.7401s + ω_{c}^{2})}

Step 5. Apply the bilinear transformation to get the digital filter: We take the bilinear transformation in the

form

s = \frac{1 − z^{−1}}{1 + z^{−1}}

It is then easy to verify that an analog quadratic factor transforms in the following manner:

\frac{b}{s^{2} + as + b} → \frac{b(1 + z^{−1})^{2}}{(1 + a + b) + (2b − 2)z^{−1} + (1 − a + b)z^{−2}}  (7.5)

while a linear factor transforms in the following manner:

\frac{ω_{c}}{s + ω_{c}} → \frac{ω_{c}(1 + z^{−1})}{1 + ω_{c} + (ω_{c} − 1)z^{−1}}  (7.6)

For this example, we can then write down the IIR digital filter transfer function

H(z) = \frac{0.003162(1 + z^{−1})^{6}}{(1.3451 − 1.7064z^{−1} + 0.94846z^{−2})(1.6886 − 1.7064z^{−1} + 0.60497z^{−2})(1.8869 − 1.7064z^{−1} + 0.40665z^{−2})}

It is customary to normalize the denominator polynomial to have leading coefficient 1, which may be achieved by normalizing each quadratic factor. This gives

H(z) = \frac{0.00073783(1 + z^{−1})^{6}}{(1 − 1.2686z^{−1} + 0.7051z^{−2})(1 − 1.0106z^{−1} + 0.3583z^{−2})(1 − 0.9044z^{−1} + 0.2155z^{−2})}

This completes the digital lowpass Butterworth filter design procedure.

For implementation of the digital filter, we do not need to multiply out the denominator to get the sixth degree polynomial. It is usually implemented as a cascade of quadratic sections, with the output from the first section acting as input to the second quadratic section, etc.

So far we have only discussed Butterworth filters in the design of lowpass filters. Other families include Chebyshev filters I and II, and elliptic filters. Their descriptions are beyond the scope of this course.

7.4 Optional Material: Frequency transformations

To be added later!