ECE 348: Digital Signal Processing La
Lab 4 (Spring 2020)
Waheed U. Bajwa
March 20, 2020
General Instructions
Please submit a report that carefully answers the problems posed in here, and include all graphs and Matla
code. Examples of all the graphs that you need to generate in this lab and include in your report are contained
in this handout. The reports must be uploaded to Sakai under Assignments within the specified time frame.
The written report will count only if you had attended and worked on the lab during the full double-period
lab session in which you are registered.
1 Filter Realization Examples
Consider the second-order IIR filter:
H(z) =
0 + b1z
−1 + b2z
−2
1 + a1z−1 + a2z−2
. (1)
The time-domain operation of this filter can be structured in different ways co
esponding to different block-
diagram realizations. See for example Figures 8.13 and 8.14 of the text. In the Direct Form I realization, the
calculation of the output signal y[n] from the input signal x[n] is implemented by the difference equation:
y[n] = b0x[n] + b1x[n− 1] + b2x[n− 2]− a1y[n− 1]− a2y[n− 2]. (2)
Assuming zero initial conditions, x[−1] = x[−2] = y[−1] = y[−2] = 0, the first two output samples
are computed by
y[0] = b0x[0],
y[1] = b0x[1] + b1x[0]− a1y[0]. (3)
Afterward, (2) can be iterated for n ≥ 2.
Another realization is the transposed (of the Direct Form II) realization, which is used internally by the
MATLAB function filter, and is based on the system of difference equations, for n ≥ 0,
y[n] = b0x[n] + v1[n]
v1[n+ 1] = b1x[n]− a1y[n] + v2[n]
v2[n+ 1] = b2x[n]− a2y[n], (4)
where v1[n] and v2[n] are internal state-variables. Note that (4) is usually initialized at v1[0] = v2[0] = 0.
The built-in function filter also keeps track of a similar internal state-variable vector and its more general
usage is as follows:
1
[y,vout] = filter(b,a,x,vin),
where vin and vout are the initial and final state vectors; that is, if (4) was iterated for 0 ≤ n ≤ N − 1, then
the last computed output and states were y[N − 1], v1[N ], v2[N ], so that
vin = [v1[0]; v2[0]] , vout = [v1[N ]; v2[N ]] .
Problem XXXXXXXXXXpoints). Write two functions, direct.m (7.5 points) and tran.m (7.5 points), that implement
(2) and (4) with usages:
y = direct(b,a,x); % Direct Form I
y = tran(b,a,x); % Transposed
[y,vout] = tran(b,a,x,vin); % with a
itrary vin
[y,vout] = tran(b,a,x); % with vin = [0; 0]
Problem XXXXXXXXXXpoints). Test your functions on the following example:
x = [1, 1, 2, 1, 2, 2, 1, 1]′ , H(z) =
4 + 2.4z−1 − 1.6z−2
1− 0.5z−1 + 0.6z−2
,
and compare their outputs to that of filter (2.5 points). Assuming vin = [0; 0], verify (2.5 points) that the
output state vector vout is the same for tran and filter.
Problem XXXXXXXXXXpoints). Using the same example, run your function tran on a sample-by-sample basis
(i.e. within a for-loop), printing with fprintf the quantities x[n], y[n], v1[n], v2[n], and generating a table as
shown below, where you must replace the * entries by their computed values.
n x y v1 v2
------------------------------------------
0 1 *.***** XXXXXXXXXX
1 1 *.***** *.****** *.******
2 2 **.***** *.****** *.******
3 1 *.***** *.****** *.******
4 2 *.***** *.****** *.******
5 2 *.***** *.****** *.******
6 1 *.***** *.****** *.******
7 1 *.***** *.****** *.******
XXXXXXXXXX485510
2
2 Notch Filter Design by Pole-Zero Placement
A simple way to design a filter with a notch at some frequency ω0 (rads/sample) is to construct a transfe
function with a zero at z0 = ejω0 , and a pole behind the zero at p0 = Rejω0 , with 0 < R < 1. The complex
conjugates z∗0 , p
∗
0 must also be included to make the filter coefficients real-valued; see Fig. 1.
Fig. 1: Notch filter design by pole-zero placement.
Thus, the transfer function is given by:
H(z) = G
(z − ejω0)(z − e−jω0)
(z −Rejω0)(z −Re−jω0)
= G
1− 2 cos(ω0) z−1 + z−2
1− 2R cos(ω0) z−1 +R2z−2
. (5)
Note that the gain factor G may be adjusted to normalize the transfer function to unity gain at DC.
The basic tradeoff for such design is that as the pole radius R gets closer to the unit circle (from inside),
the notch becomes sharper, but at the expense of a longer transient response. This tradeoff can be quantified
y the following approximate formulas for the 3-dB width ∆f in Hz and the effective 40-dB time constant
Ï„eff in seconds, where fs is the sampling rate in Hz:
∆fa =
fs
Ï€
(1−R) , τeff =
1
fs
ln(0.01)
ln(R)
. (6)
In this lab you will design two such notch filters and study their time-domain filtering properties, as well
as demonstrate the above tradeoff and the validity of the 3-dB width approximation formula. During this lab,
you may find the following anonymous MATLAB functions useful, where we note that freqz was avoided
ecause it cannot compute at a single frequency point, requiring always a frequency vector of length at least
two:
fresp = @(b,a,w) polyval(flip(b),exp(-j*w))./polyval(flip(a),exp(-j*w));
mag = @(b,a,w) abs(fresp(b,a,w));
Problem XXXXXXXXXXpoints). Design two notch filters (15 points) operating at a rate of fs = 200 Hz with a
notch frequency at f0 = 4 Hz, for the following values of R:
R = XXXXXXXXXXfilter 1),
R = XXXXXXXXXXfilter 2).
With the help of fprintf, make a table of filter coefficients (5 points) as follows, replacing the * entries
with the computed values:
3
filter 1
--------------------------------------
= [ *.****** *.****** *.****** ]
a = [ *.****** *.****** *.****** ]
filter 2
--------------------------------------
= [ *.****** *.****** *.****** ]
a = [ *.****** *.****** *.****** ]
Make a plot of the magnitude responses |H(f)| of these filters over the range 0 ≤ f ≤ 10 Hz (10 points).
See example graphs in Fig. 2, where additional features have been added as asked in Problem 2.2.
(a) (b)
XXXXXXXXXX10
0
0.2
0.4
0.6
0.8
1
magnitude response
f (Hz)
|
H
(f
)|
magnitude
f
1
, f
0
, f
2
3−dB width
XXXXXXXXXX10
0
0.2
0.4
0.6
0.8
1
magnitude response
f (Hz)
|
H
(f
)|
magnitude
f
1
, f
0
, f
2
3−dB width
Fig. 2: Plots for Problem 2.1 and Problem 2.2.
Problem XXXXXXXXXXpoints). Since the responses were normalized to unity gain at DC, the left and right 3-dB
frequencies fL, fR can be obtained by solving the following equation with respect to f :
|H(f)|2 = 1
2
⇒ |H(f)| = 1√
2
. (7)
For each filter, solve this equation with the help of fzero (5 points). You may choose the initial search point
for fzero to be the approximate left and right 3-dB frequencies obtained using the approximation of (6), that
is,
∆fa =
fs
Ï€
(1−R) ⇒ fLa = f0 −
1
2
∆fa , fRa = f0 +
1
2
∆fa.
Hint: Pass the following function handle into fzero, assuming that b, a, fs were already defined:
F = @(f) mag(b,a,2*pi*f/fs) - 1/sqrt(2).
4
For each filter, calculate the exact and approximate 3-dB frequencies in Hz and the co
esponding per-
centage e
or (5 points), and make a table (2.5 points) like the one shown below (for the case of filter #2):
exact approx
--------------------
fL = XXXXXXXXXX
fR = XXXXXXXXXX
-----------------------
percent e
or = 0.0011%
The percent e
or may be defined in terms of the L1 norms:
e
or = 100
|fL − fLa|+ |fR − fRa|
|fL|+ |fR|
.
Moreover, add the exact 3-dB frequency points to the graphs of Problem 2.1 at the 3-dB level (5 points) (see
Fig. 2).
Problem XXXXXXXXXXpoints). In this part we study the time properties of the filters. In order to observe the
notching behavior as well as the transients of the filters, we construct an input signal consisting of a 4-sec
segment of the sinusoid f0, sandwiched between 4-sec segments of two sinusoids of frequencies f1 = 2 Hz,
and f2 = 6 Hz. Because f1, f2 lie in the passband of the filters, they will pass through, whereas the middle
f0 segment will be notched out after the transients die out.
First, add the two frequency points f1, f2 to the magnitude response graphs of Problem 2.1 in Fig. 2
(2.5 points). Then construct the following signal of duration of 12 sec:
x(t) =

cos(2πf1t), 0 ≤ t < 4 sec,
cos(2πf0t), 4 ≤ t < 8 sec,
cos(2πf2t), 8 ≤ t ≤ 12 sec.
and sample it a rate of fs, that is, replace t by tn = nTs, n = 0, 1, 2 . . . , where Ts = 1/fs (5 points). Plot
x(t) versus 0 ≤ t ≤ 12 sec (2.5 points).
Then, using your function tran of Problem 1.1, filter this signal through both notch filters (5 points) and,
on two separate graphs, plot the co
esponding output signals (2.5 points). Moreover, calculate the effective
time constants given in (6) and discuss whether they are consistent with what you observe in your output
graphs (7.5 points). Make a table of calculated values as shown below (2.5 points):
R t_eff
--------------------
0.980 *.**** (sec)
0.995 *.**** (sec)
5
(a)
XXXXXXXXXX
t (sec)
-2
-1
0
1
2
x
(t
)
Input Signal
(b) (c)
XXXXXXXXXX
t (sec)
-2
-1
0
1
2
y
(t
)
Notch Filter Output, R = 0.980
XXXXXXXXXX
t (sec)
-2
-1
0
1
2
y
(t
)
Notch Filter Output, R = 0.995
Fig. 3: Plots for Problem 2.3.
6