IIR Filter Design

3.11 IIR Filter Design

There are two main categories of IIR filter design: direct and indirect methods. In the latter category, one first designs a continuous-time system function $H (s)$ , and in a second step converts $H (s)$ into a discrete-time filter $H (z)$ . In contrast, the direct methods obtain $H (z)$ directly from the specifications.

3.11.1 Direct IIR filter design

The direct methods for IIR design are tightly related to parametric power spectrum estimation. For instance, the Yule-Walker IIR filter design method is related to the techniques discussed in Section 4.8.

The Yule-Walker method is implemented in the Matlab function yulewalk. Its input is a simplified description of the desired discrete-time frequency response magnitude $| H (e^{j Ω}) |$ , which can be arbitrary. In other words, $| H (e^{j Ω}) |$ is not constrained to be the standard lowpass, highpass, bandpass or bandstop responses, which is an advantage of this method. Yule-Walker basically consists in discretizing the specification of $| H (e^{j Ω}) |$ to obtain $| H [k] |$ , finding the inverse DFT of $| H [k] |$ and solving modified Yule-Walker equations to obtain the filter coefficients [FP84]. The phase information cannot be specified. The code below illustrates how Matlab can design a digital filter with two passbands and order 16:

1N = 16; %filter order 
2f = [0 0.3 0.3 0.6 0.6 0.8 0.8 1]; %frequencies 
3m = [1 1   0   0   1   1   0   0]; %magnitudes (linear scale) 
4[B,A] = yulewalk(N,f,m); %get H(z)=B(z)/A(z) 
5[h,w] = freqz(B,A,128); %frequency response with 128 points 
6plot(w/pi,20*log10(abs(h))) %plot magnitude in dB

3.11.2 Indirect IIR filter design

When the desired filter should have standard lowpass, highpass, bandpass or bandstop responses, indirect design methods are typically adopted. The reason is that there are many well-established recipes for designing analog filters $H (s)$ and also methods for mapping $H (s)$ into a discrete-time $H (z)$ filter.

The indirect IIR filter design has the following steps:

f 1.: convert the filter specification from discrete to continuous-time,
f 2.: design the continuous-time $H (s)$ ,
f 3.: convert $H (s)$ into $H (z)$ .

Some practical analog filter designs (the step that obtains $H (s)$ ) are discussed in Example 3.14. The emphasis in the next paragraphs is the last step: conversion of $H (s)$ into $H (z)$ .

Note that, after $H (s)$ is converted into $H (z)$ , the frequency response $H (ω)$ of the former is mapped to $H (e^{j Ω})$ . When the goal of $H (z)$ is to perform filtering, it is often of interest to compare $H (e^{j Ω})$ and $H (ω)$ . In this case, as discussed in Section 3.5.7.0, the D/S conversion of the discrete-time impulse response corresponding to $H (e^{j Ω})$ can be used to obtain $H (e^{jω T_{s}})$ in rad/s, which can then be more directly compared with the original $H (ω)$ that is also in rad/s. In such situation, it may be worth to use subscripts $H_{s} (ω)$ and $H_{z} (e^{jω T_{s}})$ to identify the original and D/S converted frequency responses, respectively. Given that $H_{z} (e^{jω T_{s}})$ is periodic, a reasonable goal may be to seek a good approximation over the first Nyquist zone, i. e.

H_{z} (e^{jω T_{s}}) \approx H_{s} (ω), over 0 \leq ω < π F_{s} .

(3.56)

Five methods to convert a continuous-time $H (s)$ into a discrete-time $H (z)$ transfer function are briefly discussed in the sequel with an emphasis on the bilinear method, which is the most used in designing standard filters (lowpass, highpass, etc.).

3.11.3 Methods to convert continuous into discrete-time system functions

Matched Z-transform or pole-zero matching

The matched Z-transform method, also called pole-zero matching method, consists in simply using Eq. (2.59) to map all finite²¹ poles and zeros of $H (s)$ into new poles and zeros, in the z-plane, which then compose $H (z)$ . In other words, assuming $s_{i}$ is the $i$ -th zero or pole of $H (s)$ , it is mapped into a corresponding zero or pole of $H (z)$ given by $z_{i} = e^{s_{i} T_{s}}$ , where $T_{s}$ is the sampling interval.

After the mapping is performed, a scaling factor is used such that the magnitude of $H (z)$ is the same of $H (s)$ for a specific frequency of interest. For example, for lowpass systems, the DC gain is adjusted to have $H (z) |_{z = 1} = H (s) |_{s = 0}$ .

Example 3.18. Matched Z-transform of a first-order system. The matched Z-transform of $H (s) = a ∕ (s + a)$ is discussed in this example. The pole at $s = - a$ is mapped to a pole at $z = e^{- a T_{s}}$ such that $H (z) = 1 ∕ (z - e^{- a T_{s}})$ . Note that $H (s)$ at DC is 1, while $H (z)$ at DC is $H (z) |_{z = 1} = 1 ∕ (1 - e^{- a T_{s}}) \neq 1$ . Now, one can choose a scaling factor $g$ to create a new $H (z) = g ∕ (z - e^{- a T_{s}})$ that has $H (z) |_{z = 1} = H (s) |_{s = 0} = 1$ .

In this case, the scaling factor is $g = 1 - e^{- a T_{s}}$ and the new $H (z)$ is

H (z) = \frac{1 - e^{- a T_{s}}}{z - e^{- a T_{s}}} = \frac{(1 - e^{- a T_{s}}) z^{- 1}}{1 - e^{- a T_{s}} z^{- 1}} .

(3.57)

The extra delay $z^{- 1}$ imposed by the numerator of $H (z)$ in Eq. (3.57) is often eliminated and the final result is

H (z) = \frac{1 - e^{- a T_{s}}}{1 - e^{- a T_{s}} z^{- 1}},

(3.58)

which minimizes the delay in the corresponding difference equation. $□$

When pole-zero matching (matched Z-transform) starts from direct specification of zeros and poles, it is sometimes called pole-zero placement. The following example of designing a notch filter illustrates this situation.

Example 3.19. Using pole-zero placement, design a 2nd-order notch filter to minimize the impact of interference at a given frequency. The value $f_{0}$ of the frequency used by electrical power systems depends on the country and is typically $f_{0} = 50$ or 60 Hz. The power line signal can cause strong interference at this specific frequency and its harmonics. For instance, a circuit that digitizes an ECG signal may have to filter out the interference at $f_{0}$ to obtain better results. Here we are interested in using pole-zero placement to design a notch filter to reduce the power at $f_{0}$ . Four project guidelines will be used:

The filter will have two complex-conjugate zeros $z_{1}$ and $z_{2}$ and two complex-conjugate poles $p_{1}$ and $p_{2}$ .
The digital frequencies $Ω$ (the angles) of poles and zeros are the same. Given Eq. (1.36), which can be written as $2 πf = Ω F_{s}$ , the frequency $f_{0}$ specified in Hz can be mapped to the digital domain as $Ω_{0} = 2 π f_{0} ∕ F_{s}$ , in radians.
The zeros have magnitudes equal to one because they must be located on the unit circle.

The poles magnitude $r$ depend on $BW$ , which is the 3 dB filter bandwidth. More specifically, $r$ can be obtained from Eq. (3.35) as

r = e^{- BW π ∕ F_{s}},

or from Eq. (3.36), which leads to

r \approx 1 - \frac{BW}{F_{s}} π .

The latter is a good approximation when $r$ is in the range $[0.9, 1 [$ .

These guidelines lead to the following pairs of zeros: $z_{1} = e^{j Ω_{0}}$ and $z_{2} = e^{- j Ω_{0}}$ , while the poles are $p_{1} = r e^{j Ω_{0}}$ and $p_{2} = r e^{- j Ω_{0}}$ .

Listing 3.19 provides an example of notch filter design assuming $f_{0} = 60$ Hz.

Listing 3.19: MatlabOctaveCodeSnippets/snip_systems_notch.m. [ Python version]

1f_notch=60.0; % notch frequency in Hz 
2BW_3dB=4;  % 3 dB bandwidth in Hz 
3Fs=400; % sampling frequency in Hz 
4Omega_0=2*pi*(f_notch/Fs); % digital angular frequency 
5r=exp(-(BW_3dB/Fs)*pi); % poles radius 
6notch_zeros = [1*exp(1i*Omega_0) 1*exp(-1i*Omega_0)]; 
7notch_poles = [r*exp(1i*Omega_0) r*exp(-1i*Omega_0)]; 
8Bz_notch=poly(notch_zeros); % numerator B(z) 
9Az_notch=poly(notch_poles); % denominator A(z) 
10freqz(Bz_notch, Az_notch) % show frequency response

The frequency response shows that this filter can attenuate the signal by more than 20 dB around the frequency $f_{0}$ . However, when one decreases $BW$ to avoid filtering out components other than $f_{0}$ , the poles get closer to the unit circle, which may lead to unstable filters when implemented in an embedded system. For instance, if $BW = 4$ Hz, then $r = 0.9691$ . If $BW$ is decreased to $BW = 0.1$ Hz, then $r = 0.9992$ and the filter attenuation is rather small because the poles are almost canceling the zeros. $□$

Similar to the matched Z-transform, the next filter design techniques assume $H (s)$ is available.

Impulse invariance

The impulse invariant transformation method is based on the criterion of sampling the continuous-time impulse response such that $h [n] = T_{s} h (n T_{s})$ , where $h (t)$ and $h [n]$ are the impulse responses of $H (s)$ and $H (z)$ , respectively, with

H (z) = T_{s} Z {h (t) |_{t = n T_{s}}}

(3.59)

and $Z {\cdot}$ denoting the Z-transform.

The impulse invariant method is restricted to strictly proper $H (s)$ (degree of numerator smaller than the denominator’s). Besides, special care must be exercised when $H (s)$ has multiple-order poles.²²

Example 3.20. Example of $H (s)$ to $H (z)$ conversion using impulse invariance. Assume that $T_{s} = 0.5$ and $H (s) = 4 ∕ [(s + 2) (s + 3)]$ . Performing a partial fraction expansion leads to

H (s) = \frac{4}{s + 2} + \frac{- 4}{s + 3},

which, assuming a causal system (to define the ROC), corresponds to

h (t) = L^{- 1} {H (s)} = 4 [e^{- 2 t} - e^{- 3 t}] u (t) .

The impulse invariance is obtained after a S/D process. Using the simplified notation of Example 1.33, one can write

h [n] = T_{s} h (t) |_{t = n T_{s}} = 0.5 \times 4 [{(e^{- 2 T_{s}})}^{n} - {(e^{- 3 T_{s}})}^{n}] u [n]

and then, using a procedure similar to the one in Example 2.23, leads to

\begin{align} H (z) = Z {h [n]} & = \frac{2}{1 - e^{- 2 T_{s}} z^{- 1}} - \frac{2}{1 - e^{- 3 T_{s}} z^{- 1}} = \frac{2}{1 - 0.3679 z^{- 1}} - \frac{2}{1 - 0.2231 z^{- 1}} \\ = \frac{0.2895 z^{- 1}}{1 - 0.5910 z^{- 1} + 0.0821 z^{- 2}} . & (3.60) \end{align}

Note that the poles $s = - 2$ and $s = - 3$ in the s-plane were converted to $z = e^{- 2 T_{s}}$ and $z = e^{- 3 T_{s}}$ in z-plane, respectively, which corresponds to using Eq. (2.59) of the matched Z-transform method. And the numerators of the parcels of $H (z)$ are the pole residues scaled by $T_{s}$ . Listing 3.20 illustrates the use of impinvar in this example. In this specific case, the impulse invariant and matched Z-transform methods lead to the same poles of $H (z)$ .

Listing 3.20: MatlabOctaveCodeSnippets/snip_systems_impulseinvariance

1Ts=0.5; %sampling period 
2Bs=4; As=poly([-2 -3]); % Define H(s)=Bs/As 
3[Bz,Az] = impinvar(Bs,As,1/Ts);%impulse invariance conversion to H(z) 
4t=0:Ts:10; ht=4*(exp(-2*t)-exp(-3*t)); %create h(t) to compare 
5hn=impz(Bz,Az,length(ht)); %discrete-time h[n] from H(z)=Bz/Az 
6h_Error=Ts*ht(:) - hn(:); %error from column vectors 
7disp(['mean square error (MSE) = ' num2str(mean(h_Error.^2))]); 
8plot(t,Ts*ht); hold on; plot(t,hn,'rx'); %compare curves

The result of Listing 3.20 is a MSE = 7.7342e-34, which indicates that $h [n] \approx T_{s} h (n T_{s})$ , as desired. $□$

When $H (s)$ has only single-order poles, the impulse invariant transformation method can be summarized as:

f 1.: a partial fraction expansion of $H (s)$ to obtain the residue $r_{i}$ of each pole $p_{i}$ and
f 2.: create the parcels of $H (z)$ with poles $z = e^{p_{i} T_{s}}$ and respective numerators $T_{s} r_{i}$ .

Step invariance

The step invariant transformation method is similar to the impulse invariant but uses the step responses in continuous $q (t)$ and discrete-time $q [n]$ instead of the impulses responses. Similar to the impulse response $h (t)$ , which is the system output when the input is $δ (t)$ , the step response $q (t)$ is the system output when the input is $u (t)$ . The goal of this transformation is to achieve

q [n] = q (t) |_{t = n T_{s}},

(3.61)

which can be written as

Z^{- 1} {H (z) \frac{1}{1 - z^{- 1}}} = L^{- 1} {{H (s) \frac{1}{s}} |}_{t = n T_{s}},

(3.62)

where the S/D process is denoted in the simplified notation of Example 1.33.

Eq. (3.62) used the fact that $1 ∕ s$ and $1 ∕ (1 - z^{- 1})$ are the Laplace and Z transforms of $u (t)$ and $u [n]$ , respectively. Taking the Z-transform of both sides of Eq. (3.62) leads to

H (z) = Z {L^{- 1} {{H (s) \frac{1}{s}} |}_{t = n T_{s}}} (1 - z^{- 1}) .

(3.63)

The script snip_systems_stepinvariance.m illustrates the step invariance method for a single pole $H (s) = a ∕ (s + a)$ .

Backward difference and forward difference

The backward difference and the forward difference methods consist in substituting $s$ by $(1 - z^{- 1}) ∕ T_{s}$ and $(1 - z^{- 1}) ∕ (T_{s} z^{- 1})$ , respectively. They are motivated by the solution of differential equations using finite differences. For digital filter design these methods are outperformed by others, and have more pedagogical than practical importance.

Bilinear or Tustin’s

In the bilinear transform (also called Tustin’s) method, $s$ is substituted by $(2 ∕ T_{s}) (z - 1) ∕ (z + 1)$ .

Before discussing the bilinear transformation in more details, the next subsection provides an example of each of the major IIR design techniques.

3.11.4 Summary of methods to convert continuous-time system function into discrete-time

In previous sections, several methods to convert $H (s)$ into $H (z)$ were discussed. Table 3.6 summarizes their application to a first order analog $H (s) = a ∕ (s + a)$ , with $a \in ℝ$ . Recall from Figure 3.19, that for $H (s) = a ∕ (s + a)$ , the pole center frequency is $ω_{0} = 0$ , and the natural and cutoff frequencies $ω_{n} = ω_{c} = a$ coincide.

The matched Z-transform in Table 3.6 used Eq. (3.58).

Table 3.6: Methods to convert

H (s)

into

H (z)


Method	Mapping	$H (z)$ when
		$H (s) = \frac{a}{s + a}$

Matched Z-transform	Map poles or zeros at	$\frac{(1 - e^{- a T_{s}}) z^{- 1}}{1 - e^{- a T_{s}} z^{- 1}}$
	$s = a$ to $z = e^{a T_{s}}$

Impulse-invariant	$H (z) = T_{s} Z {L^{- 1} {{H (s)} \|}_{t = n T_{s}}}$	$\frac{T_{s} a}{1 - e^{- a T_{s}} z^{- 1}}$

Step-invariant	$H (z) = Z {L^{- 1} {{H (s) \frac{1}{s}} \|}_{t = n T_{s}}} (1 - z^{- 1})$	$\frac{(1 - e^{- a T_{s}}) z^{- 1}}{1 - e^{- a T_{s}} z^{- 1}}$

Backward difference	$s = \frac{1 - z^{- 1}}{T_{s}}$	$\frac{a}{\frac{1 - z^{- 1}}{T_{s}} + a}$

Bilinear	$s = \frac{2}{T_{s}} (\frac{1 - z^{- 1}}{1 + z^{- 1}})$	$\frac{a}{\frac{2}{T_{s}} (\frac{1 - z^{- 1}}{1 + z^{- 1}}) + a}$

Due to its importance, the bilinear transformation is further discussed in the next sections. Section 3.12 concerns the definition of the bilinear transformation and its main properties. The discussion in Section 3.12 considers that $H (s)$ and $F_{s}$ have already been chosen and focuses on finding $H (z)$ through the bilinear transformation, and understanding this mapping. Section 3.13 complements this discussion by detailing the design of system functions $H (z)$ using the bilinear transformation. This includes choosing $F_{s}$ , designing $H (s)$ and finally obtaining $H (z)$ .

²¹ As discussed in [HCJK92], a variation of the method maps poles and zeros at $\infty$ to $z = - 1$ .

²² See [Cav96, Nel07], which discuss pitfalls of Matlab’s impinvar function. You can also try using As=poly([-2 -2 -3]) in Listing 3.20 to observe that impinvar does not work well with a multiple-order pole ( $s = - 2$ in this case).