[5.10 MATLAB: DISCRETE-TIME](#page-11-0) IIR FILTERS

← Back to LINEAR SYSTEMS AND SIGNALS Overview

5.10 MATLAB: DISCRETE-TIME IIR FILTERS

Recent technological advancements have dramatically increased the popularity of discrete-time filters. Unlike their continuous-time counterparts, the performance of discrete-time filters is not affected by component variations, temperature, humidity, or age. Furthermore, digital hardware is easily reprogrammed, which allows convenient change of device function. For example, certain digital hearing aids are individually programmed to match the required response of a user.

Typically, discrete-time filters are categorized as infinite-impulse response (IIR) or finite-impulse response (FIR). A popular method to obtain a discrete-time IIR filter is by transformation of a corresponding continuous-time filter design. MATLAB greatly assists this process. Although discrete-time IIR filter design is the emphasis of this section, methods for discrete-time FIR filter design are considered in Sec. 9.7.

5.10-1 Frequency Response and Pole-Zero Plots

Frequency response and pole-zero plots help characterize filter behavior. Similar to continuous-time systems, rational transfer functions for realizable LTID systems are represented in the z-domain as

$$H[z] = \frac{Y[z]}{X[z]} = \frac{B[z]}{A[z]} = \frac{\sum_{k=0}^{N} b_k z^{-k}}{\sum_{k=0}^{N} a_k z^{-k}} = \frac{\sum_{k=0}^{N} b_k z^{N-k}}{\sum_{k=0}^{N} a_k z^{N-k}}$$

(5.54)

When only the first (N¹ + 1) numerator coefficients are nonzero and only the first (N² + 1) denominator coefficients are nonzero, Eq. (5.54) simplifies to

$$H[z] = \frac{Y[z]}{X[z]} = \frac{B[z]}{A[z]} = \frac{\sum_{k=0}^{N_1} b_k z^{-k}}{\sum_{k=0}^{N_2} a_k z^{-k}} = \frac{\sum_{k=0}^{N_1} b_k z^{N_1 - k}}{\sum_{k=0}^{N_2} a_k z^{N_2 - k}} z^{N_2 - N_1}$$

(5.55)

The form of Eq. (5.55) has many advantages. It can be more efficient than Eq. (5.54); it still works when N¹ = N² = N; and it more closely conforms to the notation of built-in MATLAB discrete-time signal-processing functions.

The right-hand side of Eq. (5.55) is a form that is convenient for MATLAB computations. The frequency response H[ej] is obtained by letting z = ej, where has units of radians. Often, = ωT, where ω is the continuous-time frequency in radians per second and T is the sampling period in seconds. Defining length-(N² + 1) coefficient vector A = [a0,a1,…,aN² ] and length-(N1+1) coefficient vector B = [b0,b1,…,bN¹ ], program CH5MP1 computes H[ej] by using Eq. (5.55) for each frequency in the input vector .

function [H] = CH5MP1(B,A,Omega);
% CH5MP1.m : Chapter 5, MATLAB Program 1
% Function M-file computes frequency response for LTID systems
% INPUTS: B = vector of feedforward coefficients
% A = vector of feedback coefficients
% Omega = vector of frequencies [rad], typically -pi<=Omega<=pi
% OUTPUTS: H = frequency response
N_1 = length(B)-1; N_2 = length(A)-1;
H = polyval(B,exp(1j*Omega))./polyval(A,exp(1j*Omega)).*...
   exp(1j*Omega*(N_2-N_1));

Note that owing to MATLAB’s indexing scheme, A(k) corresponds to coefficient ak−¹ and B(k) corresponds to coefficient bk−1. It is also possible to use the signal-processing toolbox function freqz to evaluate the frequency response of a system described by Eq. (5.55). Under special circumstances, the control system toolbox function bode can also be used.

Program CH5MP2 computes and plots the poles and zeros of an LTID system described by Eq. (5.55), again using vectors B and A.

function [p,z] = CH5MP2(B,A);
% CH5MP2.m : Chapter 5, MATLAB Program 2
% Function M-file computes and plots poles and zeros for LTID systems
% INPUTS: B = vector of feedforward coefficients
% A = vector of feedback coefficients
N_1 = length(B)-1; N_2 = length(A)-1;
p = roots([A,zeros(1,N_1-N_2)]); z = roots([B,zeros(1,N_2-N_1)]);
ucirc = exp(1j*linspace(0,2*pi,200)); % Compute unit circle for plot
plot(real(p),imag(p),'xk',real(z),imag(z),'ok',real(ucirc),imag(ucirc),'k:');
xlabel('Real'); ylabel('Imag');
ax = axis; dx = 0.05*(ax(2)-ax(1)); dy = 0.05*(ax(4)-ax(3));
axis(ax+[-dx,dx,-dy,dy]); axis equal;

The right-hand side of Eq. (5.55) helps explain how the roots are computed. When N¹ = N2, the term zN2−N¹ implies additional roots at the origin. If N¹ > N2, the roots are poles, which are added by concatenating A with zeros(N_1-N_2,1); since N² − N¹ ≤ 0, zeros(N_2-N_1,1) produces the empty set and B is unchanged. If N² > N1, the roots are zeros, which are added by concatenating B with zeros(N_2-N_1,1); since N¹ − N² ≤ 0, zeros(N_1-N_2,1) produces the empty set and A is unchanged. Poles and zeros are indicated with black x’s and o’s, respectively. For visual reference, the unit circle is also plotted. The last two lines in CH5MP2 expand the plot axis box so that root locations are not obscured and also ensure that the real and imaginary axes are drawn to the same scale.

5.10-2 Transformation Basics

Transformation of a continuous-time filter to a discrete-time filter begins with the desired continuous-time transfer function

$$H(s) = \frac{Y(s)}{X(s)} = \frac{B(s)}{A(s)} = \frac{\sum_{k=0}^{M} b_{k+N-M} s^{M-k}}{\sum_{k=0}^{N} a_k s^{N-k}}$$

As a matter of convenience, H(s) is represented in factored form as

$$H(s) = \frac{b_{N-M}}{a_0} \frac{\prod_{k=1}^{M} (s - z_k)}{\prod_{k=1}^{N} (s - p_k)}$$

(5.56)

where zk and pk are the system poles and zeros, respectively.

A mapping rule converts the rational function H(s) to a rational function H[z]. Requiring that the result be rational ensures that the system realization can proceed with only delay, sum, and multiplier blocks. There are many possible mapping rules. For obvious reasons, good transformations tend to map the ω axis to the unit circle, ω = 0 to z = 1, ω = ∞ to z = −1, and the left half-plane to the interior of the unit circle. Put another way, sinusoids map to sinusoids, zero frequency maps to zero frequency, high frequency maps to high frequency, and stable systems map to stable systems.

568 CHAPTER 5 DISCRETE-TIME SYSTEM ANALYSIS USING THE Z-TRANSFORM

Section 5.9 suggests that the z-transform can be considered to be a Laplace transform with a change of variable z = esT or s = (1/T)lnz, where T is the sampling interval. It is tempting, therefore, to convert a continuous-time filter to a discrete-time filter by substituting s = (1/T)lnz into H(s), or H[z] = H(s)|^s=(1/T)ln*^z*. This approach is impractical, however, since the resulting H[z] is not rational and therefore cannot be implemented by using standard blocks. Although not considered here, the so-called matched-z transformation relies on the relationship z = esT to transform system poles and zeros, so the connection is not completely without merit.

5.10-3 Transformation by First-Order Backward Difference

Consider the transfer function H(s) = Y(s)/X(s) = s, which corresponds to the first-order continuous-time differentiator

$$y(t) = \frac{d}{dt}x(t)$$

An approximation that resembles the fundamental theorem of calculus is the first-order backward difference

$$y(t) = \frac{x(t) - x(t - T)}{T}$$

For sampling interval T and t = nT, the corresponding discrete-time approximation is

$$y[n] = \frac{x[n] - x[n-1]}{T}$$

which has transfer function

$$H[z] = Y[z]/X[z] = \frac{1 - z^{-1}}{T}$$

This implies a transformation rule that uses the change of variable s = (1−z−¹)/T or z = 1/(1−sT). This transformation rule is appealing since the resulting H[z] is rational and has the same number of poles and zeros as H(s). Section 3.4 discusses this transformation strategy in a different way in describing the kinship of difference equations to differential equations.

After some algebra, substituting s = (1−z−¹)/T into Eq. (5.56) yields

$$H[z] = \left(\frac{b_{N-M} \prod_{k=1}^{M} (1/T - z_k)}{a_0 \prod_{k=1}^{N} (1/T - p_k)}\right) \frac{\prod_{k=1}^{M} \left(1 - \frac{1}{1 - Tz_k} z^{-1}\right)}{\prod_{k=1}^{N} \left(1 - \frac{1}{1 - Tp_k} z^{-1}\right)}$$

(5.57)

The discrete-time system has M zeros at 1/(1−Tzk) and N poles at 1/(1−Tpk). This transformation rule preserves system stability but does not map the ω axis to the unit circle (see Prob. 5.7-10).

MATLAB program CH5MP3 uses the first-order backward difference method of Eq. (5.57) to convert a continuous-time filter described by coefficient vectors A = [a0,a1,…,aN] and B = [bN−^M,bN−M+1,…,bN] into a discrete-time filter. The form of the discrete-time filter follows Eq. (5.55).

function [Bd,Ad] = CH5MP3(B,A,T); % CH5MP3.m : Chapter 5, MATLAB Program 3 % Function M-file first-order backward difference transformation

% of a continuous-time filter described by B and A into a discrete-time filter.
% INPUTS: B = vector of continuous-time filter feedforward coefficients
% A = vector of continuous-time filter feedback coefficients
% T = sampling interval
% OUTPUTS: Bd = vector of discrete-time filter feedforward coefficients
% Ad = vector of discrete-time filter feedback coefficients
z = roots(B); p = roots(A); % s-domain roots
gain = B(1)/A(1)*prod(1/T-z)/prod(1/T-p);
zd = 1./(1-T*z); pd = 1./(1-T*p); % z-domain roots
Bd = gain*poly(zd); Ad = poly(pd);

5.10-4 Bilinear Transformation

The bilinear transformation is based on a better approximation than first-order backward differences. Again, consider the continuous-time integrator

$$y(t) = \frac{d}{dt}x(t)$$

Represent signal x(t) as

$$x(t) = \int_{t-T}^{t} \frac{d}{d\tau} x(\tau) d\tau + x(t-T)$$

Letting t = nT and replacing the integral with a trapezoidal approximation yield

$$x(nT) = \frac{T}{2} \left[ \frac{d}{dt} x(nT) + \frac{d}{dt} x(nT - T) \right] + x(nT - T)$$

Substituting y(t) for (d/dt)x(t), the equivalent discrete-time system is

$$x[n] = \frac{T}{2}(y[n] + y[n-1]) + x[n-1]$$

From z-transforms, the transfer function is

$$H[z] = \frac{Y[z]}{X[z]} = \frac{2(1 - z^{-1})}{T(1 + z^{-1})}$$

The implied change of variable s = 2(1−z−1)/T(1+z−1) or z = (1+sT/2)/(1−sT/2) is called the bilinear transformation. Not only does the bilinear transformation result in a rational function H[z], the ω axis is correctly mapped to the unit circle (see Prob. 5.6-18a).

After some algebra, substituting s = 2(1−z−¹)/T(1+z−¹) into Eq. (5.56) yields

$$H[z] = \left(\frac{b_{N-M} \prod_{k=1}^{M} (2/T - z_k)}{a_0 \prod_{k=1}^{N} (2/T - p_k)}\right) \frac{\prod_{k=1}^{M} \left(1 - \frac{1 + z_k T/2}{1 - z_k T/2} z^{-1}\right)}{\prod_{k=1}^{N} \left(1 - \frac{1 + p_k T/2}{1 - p_k T/2} z^{-1}\right)} (1 + z^{-1})^{N-M} \tag{5.58}$$

In addition to the M zeros at (1+zkT/2)/(1−zkT/2) and N poles at (1+pkT/2)/(1−pkT/2), there are N−M zeros at minus 1. Since practical continuous-time filters require M ≤ N for stability, the number of added zeros is thankfully always nonnegative.

MATLAB program CH5MP4 converts a continuous-time filter described by coefficient vectors A = [a0,a1,…,aN] and B = [bN−^M,bN−M+1,…,bN] into a discrete-time filter by using the bilinear transformation of Eq. (5.58). The form of the discrete-time filter follows Eq. (5.55). If available, it is also possible to use the signal-processing toolbox function bilinear to perform the bilinear transformation.

function [Bd,Ad] = CH5MP4(B,A,T);
% CH5MP4.m : Chapter 5, MATLAB Program 4
% Function M-file bilinear transformation of a continuous-time filter
% described by vectors B and A into a discrete-time filter.
% Length of B must not exceed A.
% INPUTS: B = vector of continuous-time filter feedforward coefficients
% A = vector of continuous-time filter feedback coefficients
% T = sampling interval
% OUTPUTS: Bd = vector of discrete-time filter feedforward coefficients
% Ad = vector of discrete-time filter feedback coefficients
if (length(B)>length(A)),
    disp('Numerator order must not exceed denominator order.');
    return
end
z = roots(B); p = roots(A); % s-domain roots
gain = real(B(1)/A(1)*prod(2/T-z)/prod(2/T-p));
zd = (1+z*T/2)./(1-z*T/2); pd = (1+p*T/2)./(1-p*T/2); % z-domain roots
Bd = gain*poly([zd;-ones(length(A)-length(B),1)]); Ad = poly(pd);
As with most high-level languages, MATLAB supports general if-structures:

if expression,
       statements;
elseif expression,
       statements;
else,
       statements;
end

In the program CH5MP4, the if statement tests M > N. When true, an error message is displayed and the return command terminates program execution to prevent errors.

5.10-5 Bilinear Transformation with Prewarping

The bilinear transformation maps the entire infinite-length ω axis onto the finite-length unit circle (z = ej) according to ω = (2/T)tan(/2) (see Prob. 5.6-18b). Equivalently, = 2arctan(ωT/2). The nonlinearity of the tangent function causes a frequency compression, commonly called frequency warping, that distorts the transformation.

To illustrate the warping effect, consider the bilinear transformation of a continuous-time lowpass filter with cutoff frequency ω*^c* = 2π3000 rad/s. If the target digital system uses a sampling rate of 10 kHz, then T = 1/(10,000) and ω*^c* maps to ^c = 2arctan(ωcT/2) = 1.5116. Thus, the transformed cutoff frequency is short of the desired ^c = ωcT = 0.6π = 1.8850.

Cutoff frequencies are important and need to be as accurate as possible. By adjusting the parameter T used in the bilinear transform, one continuous-time frequency can be exactly mapped to one discrete-time frequency; the process is called prewarping. Continuing the last example, adjusting T = (2/ωc)tan(c/2) ≈ 1/6848 achieves the appropriate prewarping to ensure ω*^c* = 2π3000 maps to ^c = 0.6π.

5.10-6 Example: Butterworth Filter Transformation

To illustrate the transformation techniques, consider a continuous-time 10th-order Butterworth lowpass filter with cutoff frequency ω*^c* = 2π3000, as designed in Sec. 4.12. First, we determine continuous-time coefficient vectors A and B.

>> omega_c = 2*pi*3000; N=10;
>> poles = roots([(1j*omega_c)^(-2*N),zeros(1,2*N-1),1]);
>> poles = poles(find(poles<0));
>> B = 1; A = poly(poles); A = A/A(end);

Programs CH5MP3 and CH5MP4 are used to perform first-order forward difference and bilinear transformations, respectively.

>> Omega = linspace(0,pi,200); T = 1/10000; Omega_c = omega_c*T;
>> [B1,A1] = CH5MP3(B,A,T); % First-order backward difference transformation
>> [B2,A2] = CH5MP4(B,A,T); % Bilinear transformation
>> [B3,A3] = CH5MP4(B,A,2/omega_c*tan(Omega_c/2)); % Bilinear with prewarping

Magnitude responses are computed using CH5MP1 and then plotted.

>> H1mag = abs(CH5MP1(B1,A1,Omega));
>> H2mag = abs(CH5MP1(B2,A2,Omega));
>> H3mag = abs(CH5MP1(B3,A3,Omega));
>> plot(Omega,(Omega<=Omega_c),'k',Omega,H1mag,'k-.',...
>> Omega,H2mag,'k--',Omega,H3mag,'k:');
>> axis([0 pi -.05 1.5]);
>> xlabel('\Omega [rad]'); ylabel('Magnitude Response');
>> legend('Ideal','FOBD','BLT','Prewarp BLT','location','best');

The result of each transformation method is shown in Fig. 5.31, where FOBD and BLT stand for first-order backward difference and bilinear transformation, respectively.

Although the first-order backward difference results in a lowpass filter, the method causes significant distortion that makes the resulting filter unacceptable with regard to cutoff frequency. The bilinear transformation is better, but, as predicted, the cutoff frequency falls short of the desired value. Bilinear transformation with prewarping properly locates the cutoff frequency and produces a very acceptable filter response.

Figure 5.31 Comparison of various transformation techniques.

5.10-7 Problems Finding Polynomial Roots

Numerically, it is difficult to accurately determine the roots of a polynomial. Consider, for example, a simple polynomial that has a root at minus 1 repeated four times, (s + 1)⁴ = s⁴ + 4s³ +6s² +4s+1. The MATLAB roots command returns a surprising result:

>> roots([1464 1])'
   ans = -1.0002 -1.0000-0.0002i -1.0000+0.0002i -0.9998

Even for this low-degree polynomial, MATLAB does not return the true roots.

The problem worsens as polynomial degree increases. The bilinear transformation of the 10th-order Butterworth filter, for example, should have 10 zeros at minus 1. Figure 5.32 shows that the zeros, computed by CH5MP2 with the roots command, are not correctly located.

When possible, programs should avoid root computations that may limit accuracy. For example, results from the transformation programs CH5MP3 and CH5MP4 are more accurate if the true transfer function poles and zeros are passed directly as inputs rather than the polynomial coefficient vectors. When roots must be computed, result accuracy should always be verified.

5.10-8 Using Cascaded Second-Order Sections to Improve Design

The dynamic range of high-degree polynomial coefficients is often large. Adding the difficulties associated with factoring a high-degree polynomial, it is little surprise that high-order designs are difficult.

As with continuous-time filters, performance is improved by using a cascade of second-order sections to design and realize a discrete-time filter. Cascades of second-order sections are also more robust to the coefficient quantization that occurs when discrete-time filters are implemented on fixed-point digital hardware.

To illustrate the performance possible with a cascade of second-order sections, consider a 180th-order transformed Butterworth discrete-time filter with cutoff frequency ^c = 0.6π ≈ 1.8850. Program CH5MP5 completes this design, taking care to initially locate poles and zeros without root computations.

Figure 5.32 Pole-zero plot computed by using roots.

% CH5MP5.m : Chapter 5, MATLAB Program 5
% Script M-file designs a 180th-order Butterworth lowpass discrete-time filter
% with cutoff Omega_c = 0.6*pi using 90 cascaded second-order filter sections.
omega_0 = 1; % Use normalized cutoff frequency for analog prototype
psi = [0.5:1:90]*pi/180; % Butterworth pole angles
Omega_c = 0.6*pi; % Discrete-time cutoff frequency
Omega = linspace(0,pi,1000); % Frequency range for magnitude response
Hmag = zeros(90,1000); p = zeros(1,180); z = zeros(1,180); % Pre-allocation
for stage = 1:90,
    Q = 1/(2*cos(psi(stage))); % Compute Q for stage
    B = omega_0^2; A = [1 omega_0/Q omega_0^2]; % Compute stage coefficients
    [B1,A1] = CH5MP4(B,A,2/omega_0*tan(0.6*pi/2)); % Transform stage to DT
    p(stage*2-1:stage*2) = roots(A1); % Compute z-domain poles for stage
    z(stage*2-1:stage*2) = roots(B1); % Compute z-domain zeros for stage
    Hmag(stage,:) = abs(CH5MP1(B1,A1,Omega)); % Compute stage mag response
end
ucirc = exp(j*linspace(0,2*pi,200)); % Compute unit circle for pole-zero plot
figure;
plot(real(p),imag(p),'kx',real(z),imag(z),'ok',real(ucirc),imag(ucirc),'k:');
axis equal; xlabel('Real'); ylabel('Imag');
figure; plot(Omega,prod(Hmag),'k'); axis([0 pi -0.05 1.05]);
xlabel('\Omega [rad]'); ylabel('Magnitude Response');

The figure command preceding each plot command opens a separate window for each plot.

The filter’s pole-zero plot is shown in Fig. 5.33, along with the unit circle, for reference. All 180 zeros of the cascaded design are properly located at minus 1. The wall of poles provides an amazing approximation to the desired brick-wall response, as shown by the magnitude response in Fig. 5.34. It is virtually impossible to realize such high-order designs with continuous-time filters, which adds another reason for the popularity of discrete-time filters. Still, the design is not

Figure 5.33 Pole-zero plot for 180th-order discrete-time Butterworth filter.

Figure 5.34 Magnitude response for a 180th-order discrete-time Butterworth filter.

trivial; even functions from the MATLAB signal-processing toolbox fail to properly design such a high-order discrete-time Butterworth filter.