EL7133 Exercises
EL7133 Exercises
EL7133 Exercises
HW 1
The DFT: 1.1, 1.2, 1.3, 1.5, 1.6, 1.9, 1.13, 1.50, 1.51 1.52, 1.53
Speech filtering: 14.1,
Ivan W. Selesnick
April 11, 2013
HW 2
Filters: 3.17, 3.20
The DFT: 1.11, 1.12, 1.16, 1.18, 1.28, 1.35, 1.61 1.64, 1.66,
The FFT: 2.1, 2.3, 2.4, 2.6
Speech filtering: 14.2,
Contents
1 The Discrete Fourier Transform
18
3 Filters
19
30
5 Windows
39
51
55
8 Spectral Factorization
57
59
65
11 Multirate Systems
69
12 Quantization
75
13 Spectral Estimation
76
14 Speech Filtering
83
15 More Exercises
87
16 Old Exercises
92
HW 3
The DFT: 1.22, 1.27, 1.33, 1.55
Filters: 3.3, 3.6
Linear-phase FIR filters: 4.1, 4.2, 4.3, 4.5, 4.7, 4.8, 4.9, 4.10, 4.15
Speech filtering: 14.3
HW 4
DFT: 1.44, 1.45
Filters: 3.4, 3.9, 3.29
Linear-phase FIR filters: 4.13, 4.14,
Windows: 5.1, 5.2, 5.3, 5.5, 5.6, 5.7, 5.8, 5.15, 5.16, 5.17, 5.18, 5.20, 5.24
Speech filtering: 14.4
HW 5
Least-square FIR filter design: 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9
HW 6
Least-square FIR filter design: 6.14
Filter specifications: 3.32
Minimax FIR filter design: 7.1, 7.2, 7.4, 7.6, 7.7, 7.8
HW 7
Spectral factorization: 8.1, 8.4
IIR filters: 10.1, 10.2, 10.4, 10.5, 10.6, 10.7, 10.8, 10.10, 10.12
HW 8
Minimum-phase FIR filter design: 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.10a
HW 9
Multirate systems: 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.15, 11.16, 11.24
HW 10
Multirate systems: 11.8, 11.9, 11.13
Finite precision effects: 12.1, 12.3, 12.4, 12.5
1
1.1 Compute the DFT of the 2-point signal by hand (without a calculator or
computer).
HW 11
Spectral estimation: 13.1, 13.2, 13.5, 13.8, 13.9, 13.10, 13.11, 13.12
Fractional delay systems: 15.5, 15.6
Other: 15.2 (optional problem)
x = [20,
5]
2,
5,
1]
1.3 The even samples of the DFT of a 9-point real signal x(n) are given by
X(0) = 3.1,
X(2) = 2.5 + 4.6 j,
X(4) = 1.7 + 5.2 j,
X(6) = 9.3 + 6.3 j,
X(8) = 5.5 8.0 j,
Determine the missing odd samples of the DFT. Use the properties of the
DFT to solve this problem.
1.4 The DFT of a 5-point signal x(n), 0 n 4 is
X(k) = [5,
6,
1,
2,
9],
0 k 4.
0 n 4.
What are the DFT coefficients G(k) of the signal g(n), for 0 k 4?
1.5 Compute by hand the circular convolution of the following two 4-point
signals (do not use MATLAB, etc.)
g = [1, 2, 1, 1]
h = [0, 1/3, 1/3, 1/3]
1.6 What is the circular convolution of the following two sequences?
x = [1 2 3 0 0 0 0];
h = [1 2 3 0 0 0 0];
1.7 What is the circular convolution of the following two sequences?
2
This vectors X1, X2, X3, X4 are shown below out of order. Using properties
of the DFT, match them to vectors A, B, C, D, by completing the table.
Explain how you get your answers. Do not use MATLAB for this problem
and do not explicitly compute the DFT; instead use the properties of the
DFT.
x = [1 2 3 0 0 0 0];
h = [0 0 0 0 1 2 3];
Signal
x1
x2
x3
x4
1.9 What is y after running the following MATLAB commands? Do not use
MATLAB to get your answer and do not explicitly compute the DFT;
instead use the properties of the DFT.
A:
clear
x = [1 2 0 -1];
g = [1 0 0 2 1];
X = fft([x 0 0 0]);
G = fft([g 0 0]);
Y = X.*G;
y = ifft(Y);
C:
16.0000
-4.5489
0.1920
-0.1431
-0.1431
0.1920
-4.5489
+
+
+
[1
[1
[1
[0
2
2
2
2
3
3
3
3
X1
X2
X3
X4
=
=
=
=
fft(x1)
fft(x2)
fft(x3)
fft(x4)
4
4
4
4
-12.4480i
+ 4.9582i
- 4.8440i
+ 4.8440i
- 4.9582i
+12.4480i
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
-12.4480i
+ 4.9582i
- 4.8440i
+ 4.8440i
- 4.9582i
+12.4480i
1.12 Let x(n), for 0 n N 1, be a real N -point signal. The DFT coefficients
are X(k), for 0 k N 1.
(a) Show that X(0) is a real number.
D:
0
0
0
0
0
0
0
clear
x = [1 2 3 4];
g = [5 6 7 8];
X = fft([x 0 0 0]);
G = fft([g 0 0]);
Y = X.*G;
y = ifft(Y);
19.0000
-5.0489
-0.3080
-0.6431
-0.6431
-0.3080
-5.0489
2.1906i
0.2408i
0.6270i
0.6270i
0.2408i
2.1906i
B:
1.10 What is y after running the following MATLAB commands?
x1
x2
x3
x4
DFT
3 2 1];
4 3 2];
-4 -3 -2];
-4 -3 -2];
1.18 What is the inverse DFT of the following N -point discrete-time sequence
X(k)?
(a) [3, 2, 1, 0, 0, 0, 0, 2, 1]
(b) [3, 2, 1, 0, 0, 0, 0, 2, 1]
(c) [3, 2, 1, 0, 0, 0, 0, 2, 1]
(d) [0, 2, 1, 0, 0, 0, 0, 2, 1]
X(k) = [0, ej , 0, 0, . . . , 0, ej ],
|
{z
}
(e) [0, 2, 1, 0, 0, 0, 0, 2, 1]
0 k N 1.
N 3 values
(f) [3, 2, 1, 0, 0, 0, 0, 1, 2]
Find the N -point signal x(n). Simplify your answer. Is x(n) real? If so,
express your answer without using j.
(g) [3, 2, 1, 0, 0, 0, 0, 1, 2]
(h) [0, 2, 1, 0, 0, 0, 0, 1, 2]
(i) [0, 2, 1, 0, 0, 0, 0, 1, 2]
Which of these signals have a real-valued 9-point DFT? Which of these
signals have an imaginary-valued 9-point DFT? Do not use MATLAB or
any computer to solve this problem and do not explicitly compute the DFT;
instead use the properties of the DFT.
X(N k) = X(k),
k = 0, 1, . . . , N 1.
x(n)
1
0.5
(1)
0.5
X(k) = [0
0],
k = 0, . . . , 12
10
15
20
(a) The 23-point DFT of x(n) is computed. The DFT coefficients are
denoted X(k). Accurately sketch |X(k)| for 0 k N 1.
1.17 What is the inverse DFT of the following N -point discrete-time sequence
X(k)?
X(k) = [0, 1, 0, 0, . . . , 0, 1],
|
{z
}
k = 0, 1, . . . , N 1.
N 3 terms
(a) What is the 30-point signal x(n) that has the following DFT? (Provide
a mathematical formula.)
(
1, k = 4
for 0 k < 30
X(k) =
0, k 6= 4
x = sin(2*pi*[0:7]/8);
and find the DFT of x. Do not use direct computation of the DFT.
where N = 20.
(a) Roughly sketch the signal for 0 n 19. Do not explicitly calculate
the values.
(b) Find the DFT coefficients of the signal. That means, find X(k) for
0 k 19. Show the derivation of your answer. You should not use
direct computation of the DFT.
x(n) = cos
n
n = 0, . . . , N 1
N
4
0,
1,
{z
repeats
0,
}
...,
1,
0,
1,
0]
The sequence is made of K periods of the 4-point sequence (1, 0, -1, 0).
The length of the sequence is N = 4K. Simplify your answer.
1.26 Sinusoids.
5
1.28 The following MATLAB commands define two ten-point signals and the
DFT of each.
x1
x2
X1
X2
=
=
=
=
cos([0:9]/9*2*pi);
cos([0:9]/10*2*pi);
fft(x1);
fft(x2);
DFT
X1
X2
(a) Roughly sketch each of the two signals, highlighting the distinction
between them.
(b) Which of the following four graphs illustrates the DFT |X1 (k)|? Explain your answer. Which graph illustrates the DFT |X2 (k)|?
GRAPH A
6
x1
x2
x3
X1
X2
X3
0
GRAPH C
6
(a) Roughly sketch each of the three signals, highlighting the values of the
first and last values of each signal
GRAPH A
GRAPH F
sin([0:9]/10*pi);
sin([0:9]/9*2*pi);
sin([0:9]/10*2*pi);
fft(x1);
fft(x2);
fft(x3);
=
=
=
=
=
=
(b) Which of the following graphs illustrates the DFT |X(k)| of each signal?
GRAPH E
6
GRAPH D
1.29 The following MATLAB commands define three ten-point signals and the
DFT of each.
GRAPH B
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9
GRAPH D
GRAPH E
GRAPH F
0
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9
GRAPH G
GRAPH H
GRAPH I
0
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9
0
2
0 1 2 3 4 5 6 7 8 9
GRAPH H
GRAPH C
GRAPH B
GRAPH G
Graph
0
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9
Give your answer by completing the table and provide an explanation for your answer. Note that there are more graphs shown than
needed. Your answer will use only three of the graphs. Provide a brief
explanation for your answer.
DFT
w(n)
1
0.8
0.6
Graph
0.4
X1
X2
X3
0.2
0
N1
f(n)
1.30 What is the result of the following MATLAB commands? Explain your
answer.
1
0.8
0.6
>> x = ones(1,10);
>> X = fft(x)
0.4
0.2
0
N1
n
1.31 Let N be an even integer. Let the N -point signal x(n) be defined for
0 n N 1 as
(
1, n even
x(n) =
0, n odd.
(a) Find the N -point DFT W (k) of the sequence w(n). Also, sketch W (k).
(b) Find the N -point DFT F (k) of the sequence f (n). Also, sketch |F (k)|.
(c) The signal x(t) = cos(2t) is sampled at interval T = 0.01 second for
N = 100 samples. The frequency of the signal is then measured using
a windowed DFT (with no zero padding). The window w(n) is used.
Find the resulting N -point DFT.
Find the DFT of the N -point discrete-time signal x(n). Your answer should
be simplified. If X(k) is real, write it without using j.
1.32 What is the result of the following MATLAB command? Do not use direct
computation of the DFT. Show your derivation or explain your answer.
>> fft([1
-1
-1
-1
0])
1.34 Let N be an even integer. What is the DFT of the following N -point
discrete-time signal x?
1.33 Two N -point sequences are created in Matlab using the commands:
n = 0:N-1;
w = 0.5 - 0.5 * cos(2*pi*n/N);
f = 0.5 - 0.5 * cos(2*pi*(n+0.5)/N);
x = [0, 0, . . . , 0, 1, 0, 0, . . . , 0]
|
{z
}
|
{z
}
N
2
terms
N
2
1 terms
1.35 DFT Matching. Match each discrete-time signal with its DFT. You
should be able to do this problem with out using a computer.
SIGNAL 1
DFT 1
DFT 2
30
30
20
20
10
10
SIGNAL 2
1.5
1.5
0.5
0.5
0.5
0.5
0
0
10
20
30
10
DFT 3
1.5
20
30
20
30
20
30
20
30
DFT 4
30
30
20
20
10
10
1.5
0
10
20
30
SIGNAL 3
10
20
30
SIGNAL 4
1.5
1.5
0.5
0.5
0.5
0.5
0
0
10
20
30
10
DFT 5
1.5
DFT 6
30
30
20
20
10
10
1.5
0
10
20
30
SIGNAL 5
10
20
30
SIGNAL 6
1.5
1.5
0.5
0.5
0.5
0.5
0
0
10
20
30
10
DFT 7
1.5
DFT 8
30
30
20
20
10
10
1.5
0
10
20
30
SIGNAL 7
1.5
0.5
0.5
0.5
0.5
1.5
1.5
10
20
30
SIGNAL 8
1.5
10
20
30
0
0
10
20
30
10
20
30
10
1.36 DFT Matching. Match each discrete-time signal with its DFT by filling
out the following table. You should be able to do this problem with out
using a computer.
DFT 1
DFT 2
25
25
20
20
15
15
10
10
0
0
10
15
20
25
DFT 3
SIGNAL 1
SIGNAL 2
1.5
1.5
0.5
0.5
0
0.5
0.5
25
20
20
15
15
10
10
20
25
20
25
20
25
20
25
0
0
10
15
20
25
DFT 5
1.5
15
DFT 4
25
0
0
10
10
15
DFT 6
25
25
20
20
15
15
10
10
1.5
0
10
15
20
25
SIGNAL 3
10
15
20
25
SIGNAL 4
1.5
1.5
0.5
0.5
0
0.5
0.5
0
0
10
15
20
25
DFT 7
1.5
10
15
DFT 8
25
25
20
20
15
15
10
10
1.5
0
10
15
20
25
SIGNAL 5
10
15
20
25
SIGNAL 6
1.5
1.5
0.5
0.5
0.5
0.5
1.5
1.5
0
10
15
20
25
SIGNAL 7
10
15
20
25
20
25
SIGNAL 8
1.5
1.5
0.5
0.5
0.5
0.5
1.5
1.5
0
10
15
20
25
0
0
10
15
10
15
20
25
10
15
1.37 DFT Matching. Match each of the following 28-point discrete-time signals with its DFT by completing a table. (The DFT plots show the magnitude of the complex DFT coefficients.)
SIGNAL 1
SIGNAL 2
1.38 DFT Matching. Match each of the following 28-point discrete-time signals with its DFT by completing a table. (The DFT plots show the magnitude of the complex DFT coefficients.)
SIGNAL 1
SIGNAL 3
SIGNAL 2
SIGNAL 3
1.5
1.5
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1.5
1.5
1.5
1.5
1.5
1.5
10
15
20
25
10
SIGNAL 4
15
20
25
10
SIGNAL 5
15
20
25
10
15
20
25
10
SIGNAL 4
SIGNAL 6
15
20
25
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1
5
10
15
20
1.5
0
10
SIGNAL 7
15
20
1.5
1.5
25
25
10
SIGNAL 8
15
20
10
15
20
25
10
15
20
25
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1
10
15
20
10
DFT 1
15
20
25
10
10
10
0
10
15
10
20
25
15
20
15
20
25
20
25
1.5
0
DFT 1
10
15
20
25
10
15
20
10
15
20
25
10
15
DFT 2
DFT 3
25
0
DFT 5
10
15
20
25
10
15
20
25
10
15
20
25
15
20
25
15
20
25
DFT 6
10
1.5
10
0
0
DFT 4
2
25
1.5
0
25
20
DFT 3
15
DFT 2
15
1.5
1.5
25
15
1.5
5
15
1.5
0.5
10
SIGNAL 9
SIGNAL 8
1.5
1.5
25
1.5
0
SIGNAL 7
SIGNAL 9
20
1.5
0
25
15
1.5
0.5
1
10
SIGNAL 6
1.5
1.5
SIGNAL 5
DFT 4
DFT 5
DFT 6
6
1
4
4
0.5
0
0
10
15
20
25
0
0
10
15
20
25
10
15
20
25
0
DFT 7
DFT 8
10
15
20
25
10
15
20
25
10
DFT 9
20
12
20
DFT 7
DFT 8
DFT 9
10
15
8
15
10
10
4
5
0
0
10
15
20
25
0
0
10
15
20
25
10
15
20
25
0
10
10
15
20
25
10
15
20
25
10
1.39 DFT Matching. Match each of the following 28-point discrete-time signals with its DFT by filling out a table. (The DFT plots show the magnitude of the complex DFT coefficients.)
SIGNAL 1
SIGNAL 2
1.40 DFT Matching. Match each of the following 28-point discrete-time signals with its DFT by completing a table. (The DFT plots show the magnitude of the complex DFT coefficients.)
SIGNAL 3
SIGNAL 1
SIGNAL 2
SIGNAL 3
1.5
1.5
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1.5
1.5
1.5
1.5
1.5
1.5
10
15
20
25
10
SIGNAL 4
15
20
25
10
SIGNAL 5
15
20
25
10
SIGNAL 6
15
20
25
10
SIGNAL 4
15
20
25
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1.5
0
10
15
20
1.5
25
10
SIGNAL 7
15
20
1.5
25
10
SIGNAL 8
15
20
10
SIGNAL 9
15
20
10
SIGNAL 7
15
20
25
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1.5
10
15
20
1.5
25
10
DFT 1
15
20
1.5
25
10
DFT 2
15
20
10
DFT 3
20
15
20
25
25
10
15
20
25
DFT 3
15
15
10
10
15
10
10
0
20
20
20
15
15
20
10
10
10
DFT 2
20
10
30
15
20
30
25
1.5
25
DFT 1
20
20
1.5
25
15
1.5
0.5
10
SIGNAL 9
1.5
SIGNAL 8
25
1.5
25
1.5
1.5
20
1.5
25
15
1.5
0.5
1
10
SIGNAL 6
1.5
1.5
SIGNAL 5
10
DFT 4
15
20
25
10
DFT 5
15
20
25
5
0
0
10
DFT 6
15
20
25
0
0
10
DFT 4
15
20
25
12
10
15
20
25
20
25
20
25
DFT 6
15
DFT 5
12
15
10
10
6
8
10
4
10
2
2
0
0
0
10
15
20
25
0
0
10
DFT 7
15
20
25
10
DFT 8
30
15
20
25
10
DFT 9
15
20
25
1.5
0.5
10
15
20
25
20
25
10
15
DFT 9
10
1
5
0
0
15
1.5
10
15
10
10
DFT 8
15
20
0
0
DFT 7
0
0
10
15
20
25
0.5
0
0
10
15
20
25
0
0
11
10
15
20
25
0
0
10
15
20
25
10
15
Do not use direct computation of the DFT. Show your use of the appropriate
DFT properties.
1.41 What is the result of the following MATLAB command? Explain your
answer.
>> fft([0 1 1 1 1 1])
Do not use direct computation of the DFT. Show your use of the appropriate
DFT properties.
1.45 What is the result of the following MATLAB commands? Explain your
answer.
N
1
X
exp j
n=0
>>
>>
>>
>>
2
nk ,
N
k Z.
x
X
G
g
=
=
=
=
[8 7 6 5 4 3 2 1];
fft(x);
[X(5:8) X(1:4)];
ifft(G)
Do not use direct computation of the DFT. Show your use of the appropriate
DFT properties.
1.46 What is the result of the following MATLAB commands? Explain your
answer.
>>
>>
>>
>>
x
X
G
g
=
=
=
=
[6 5 4 3 2 1];
fft(x);
X .* exp(-j*[0:5]*2*pi/6);
ifft(G)
1.43 What is the result of the following MATLAB commands? Explain your
answer.
>>
>>
>>
>>
x
X
G
g
=
=
=
=
[6 4 3 2 1];
fft(x);
conj(X);
ifft(G)
1.47 What is the result of the following MATLAB commands? Explain your
answer.
>>
>>
>>
>>
Do not use direct computation of the DFT. Show your use of the appropriate
DFT properties.
1.44 What is the result of the following MATLAB commands? Explain your
answer.
>>
>>
>>
>>
x
X
G
g
=
=
=
=
x
X
G
g
=
=
=
=
[6 5 4 3 2 1];
fft(x);
X .* (-1).^[0:5];
ifft(G)
Do not use direct computation of the DFT nor MATLAB. Show your use
of the appropriate DFT properties.
[6 5 4 3 2 1];
fft(x);
[X(1) X(end:-1:2)];
ifft(G)
12
5,
4,
3,
2,
1,
2,
3,
4,
1.53 DFT and linear convolution. Write a Matlab function that uses the
DFT (fft) to compute the linear convolution of two sequences that are
not necessarily of the same length. (Use zero-padding.) Verify that it
works correctly by comparing the results of your function with the Matlab
command conv.
5]
is
X = [35,
10.4721,
0,
1.5279,
0,
1,
0,
1.5279,
0,
10.4721]
Using properties of the DFT, what is the DFT of the following signal g?
g = [1,
2,
3,
4,
5,
6,
5,
4,
3,
2]
Do not use direct computation of the DFT. Show your use of the appropriate
DFT properties.
1.49 A 6-point signal x(n) has the DFT:
X(k) = [4
5],
k = 0, . . . , 5
1.54 Match each MATLAB code fragment to the figure it produces by completing
the table.
Find the DFT G(k) of the 6-point signal g(n) defined as g(n) = x(n) (1)n .
Do not use direct computation of the DFT. Show your use of the appropriate
DFT properties.
Code
1.50 DFT and deconvolution. (Porat 4.26) We are given the two sequences
x = [1, 3, 1, 5],
1
2
3
4
y = [7, 7, 9, 3].
y=h
x?
Here, h
x denotes circular convolution of h and x. If so, find h; if not,
prove that such h does not exist. (Use the DFT command in Matlab, fft.)
% CODE 2
h = [1 1 1 1 1 1 1 1];
H = fft(h,32);
stem(0:31,abs(fftshift(H)),.)
% CODE 3
h = [2 2 2 2 0 0 0 0];
H = fft(h,32);
stem(0:31,abs(H),.)
% CODE 4
h = [2 2 2 2 0 0 0 0];
H = fft(h,32);
stem(0:31,abs(fftshift(H)),.)
Hand in a hard copy of both functions, and an example verifying they give
the same results (you might use the diary command).
13
Figure
FIGURE 1
1.55 Match each MATLAB code fragment to the figure it produces. Provide a
brief explanation for your answers.
FIGURE 2
N = 40;
n = 0:N-1;
% CODE 4
h = [0.25 0.5 0.25];
H = fft(h,N);
stem(n,abs(H),.)
0
0
10
20
30
10
FIGURE 3
20
30
% CODE 1
h = hamming(N)/N;
H = fft(h,N);
stem(n,abs(H),.)
FIGURE 4
% CODE 2
h = [0.5 0.5];
H = fft(h,N);
stem(n,abs(H),.)
0
0
10
20
30
10
20
% CODE 5
h = [-0.25 0.5 -0.25];
H = fft(h,N);
stem(n,abs(H),.)
30
% CODE 6
h = ones(1,4)/4;
H = fft(h,N);
stem(n,abs(H),.)
% CODE 3
h = [0.5 -0.5];
H = fft(h,N);
stem(n,abs(H),.)
FIGURE 1
FIGURE 2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
10
20
30
40
10
FIGURE 3
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
10
20
30
40
10
FIGURE 5
1
1
0.8
0.6
0.6
0.4
0.4
0.2
0.2
14
10
20
30
40
20
30
40
30
40
FIGURE 6
0.8
20
FIGURE 4
30
40
10
20
1.56 Match each MATLAB code fragment to the figure it produces. Provide a
brief explanation for your answers.
% CODE 1
h = [1 1 1 1 1 1 1 1];
H = fft(h,40);
stem(0:39,abs(H),.)
% CODE 4
h = [2 2 2 2];
H = fft(h,40);
stem(0:39,abs(fftshift(H)),.)
% CODE 2
h = [1 1 1 1 1 1 1 1];
H = fft(h,40);
stem(0:39,abs(fftshift(H)),.)
% CODE 5
h = conv([1 1 1 1],[1 1 1 1])/2;
H = fft(h,40);
stem(0:39,abs(H),.)
% CODE 3
h = [2 2 2 2];
H = fft(h,40);
stem(0:39,abs(H),.)
% CODE 6
h = conv([1 1 1 1],[1 1 1 1])/2;
H = fft(h,40);
stem(0:39,abs(fftshift(H)),.)
FIGURE 1
1.57 Match each MATLAB code fragment to the figure it produces. Provide a
brief explanation for your answers.
N = 40;
n = 0:N-1;
% CODE 4
h = hamming(N)/N .* (-1).^n;
H = fft(h,N);
stem(n,abs(H),.)
% CODE 1
h = [1/3 1/3 1/3];
H = fft(h,N);
stem(n,abs(H),.)
% CODE 5
h = ones(1,N/4)/(N/4);
H = fft(h,N);
stem(n,abs(H),.)
% CODE 2
h = [1/3 -1/3 1/3];
H = fft(h,N);
stem(n,abs(H),.)
% CODE 6
h = ones(1,N/2)/(N/2);
H = fft(h,N);
stem(n,abs(H),.)
% CODE 3
h = hamming(N)/N;
H = fft(h,N);
stem(n,abs(H),.)
FIGURE 2
0.8
0.6
FIGURE 1
10
20
30
40
FIGURE 2
1
0.8
0.6
10
20
30
0.4
40
0.4
FIGURE 3
FIGURE 4
0.2
0.2
0
0
10
20
30
40
10
FIGURE 3
4
0.8
0.6
10
20
30
40
20
30
40
30
40
30
40
FIGURE 4
1
0.8
0.6
10
20
30
0.4
40
0.4
FIGURE 5
FIGURE 6
0.2
0.2
0
0
10
20
30
40
10
FIGURE 5
4
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
10
20
30
40
10
20
30
20
FIGURE 6
40
0
0
15
10
20
30
40
10
20
1.58 In this problem, simple sequences are zero-padded and their FFTs are computed. Some are displayed using fftshift. Match each MATLAB code
fragment to the figure it produces by completing a table.
% CODE 1
h = [1 1 1 1 1 1 1 1];
L = 2^8;
H = fft(h,L);
plot(0:L-1,abs(H))
% CODE 5
h = conv([1 1 1 1],[1 1 1 1])/2;
H = fft(h,L);
plot(0:L-1,abs(H))
FIGURE 1
8
100
150
200
250
50
FIGURE 3
8
50
100
150
200
250
50
100
150
1.61 An analog signal is sampled at 8192 Hz and 600 samples are collected. These
600 samples are available on the course webpage as the file signal1.txt.
Plot the signal versus time in seconds; and using the DFT, plot the spectrum
of the signal versus physical frequency in hertz with the DC component in
the center of the plot. If your computer has sound capability, use the
soundsc command to listen to the signal.
100
150
200
250
50
100
150
200
250
FIGURE 6
8
for n = 3, . . . , 3. g(n) is zero outside this range. Note that g(n) = g(n).
What effect does this have on your result?
FIGURE 5
8
g(n) = (1, 2, 4, 8, 4, 2, 1)
FIGURE 4
FIGURE 2
50
% CODE 6
h = conv([1 1 1 1],[1 1 1 1])/2;
H = fft(h,L);
plot(0:L-1,abs(fftshift(H)))
% CODE 3
h = [2 2 2 2];
H = fft(h,L);
plot(0:L-1,abs(H))
x(n) = (1, 2, 3, 2, 3, 4, 2)
% CODE 4
h = [2 2 2 2];
H = fft(h,L);
plot(0:L-1,abs(fftshift(H)))
% CODE 2
h = [1 1 1 1 1 1 1 1];
H = fft(h,L);
plot(0:L-1,abs(fftshift(H)))
1.60 DFT and samples of the DTFT. If a discrete-time signal x(n) is finite
in length, then samples of its DTFT can be computed using the DFT. Let
x(n) be given by
200
250
50
100
150
200
250
(a) Which DFT index k is nearest to 120 Hz, and what is its physical
frequency in hertz?
1.59 When using both windowing and zero padding for an N -point sequence
x(n), which is the proper procedure: (1) zero-pad first to length M , then
use a window of length M or (2) use a window of length N on x(n) and
then zero-pad to length M ? Give reasons.
(b) What is the minimum number of zeros we must pad onto the 980
samples to obtain a DFT value at 120 Hz exactly? What is the DFT
index k then corresponding to 120 Hz?
16
(b) The individual performing this operation expects that the DFT coefficients X(k) corresponding to the frequency range 30 Hz < f < 50 Hz
will be zero, because the analog signal xa (t) was band-limited to 30 Hz
and the signal is over-sampled (no aliasing takes place). However, the
DFT coefficients in that frequency range are not exactly zero. Explain
clearly why those DFT coefficients are not exactly zero.
The signal x(t) is sampled with a sampling rate of 50 Hz and 1000 samples
are collected. You then take the DFT of these 1000 samples. Which DFT
coefficients X d (k) are free of aliasing?
1.67 The analog signal x(t) is band-limited to 40 Hz. Suppose the signal is
sampled at the rate of 100 samples per second and that at this rate 200
samples are collected. Then 200 zeros are appended to the 200 samples to
form a 400-point vector. Then the 400-point DFT of this vector is computed
to get X(k) for 0 k 399.
(a) Which DFT coefficients are free of aliasing?
1.70 DFT. A 3 Hz continuous-time cosine signal with amplitude 1 is sampled
at a rate of 10 samples per second to produce a discrete-time signal x(n),
n Z.
(b) The DFT coefficient X(50) represents the spectrum of the analog signal at what frequency f ? (Give your answer in Hz).
1.68 An analog signal s(t) only contains frequencies between 30 Hz and 40 Hz.
(It is a bandpass signal). You want to know the value of the spectrum at
35 Hz.
The signal is sampled at a rate of 60 samples per second and 200 samples
are collected. Call them x(n) for n = 0, . . . , 199. You then take a 200-point
DFT of these 200 samples. Call the DFT values X(k), for k = 0, . . . , 199.
(c) The 10 point vector g(n) in part (b) can be extended with zeros outside
the range 0 n 9 to define the signal v(n), n Z,
(
g(n) 0 n 9
v(n) =
0
otherwise
1.69 An analog signal xa (t), band-limited to 30 Hz, is sampled 100 times per
second for 10 seconds. The DFT of the collected data is computed and
1000 DFT coefficients, X(k), are obtained.
(a) What frequency (in Hz) does the DFT coefficient X(100) correspond
to?
x = rand(1,126);
g = rand(1,126);
% METHOD 1
y = conv(x,g);
%
X
G
Y
y
METHOD 2
= fft([x zeros(1,120)]);
= fft([g zeros(1,120)]);
= X.*G;
= ifft(Y);
%
X
G
Y
y
METHOD 3
= fft([x zeros(1,125)]);
= fft([g zeros(1,125)]);
= X.*G;
= ifft(Y);
%
X
G
Y
y
METHOD 4
= fft([x zeros(1,130)]);
= fft([g zeros(1,130)]);
= X.*G;
= ifft(Y);
Comment on your results. For which values N is the MATLAB FFT algorithm most and least efficient?
On recent versions of MATLAB, the flops command is no longer available.
In this case just time the fft function with the MATLAB command etime
or with the commands tic and toc.
2.4 The following MATLAB code gives four methods to compute the linear convolution of x(n) and g(n). Identify which of the four methods are wrong!
Of the remaining methods, which one do you think requires the most additions and multiplications (the least efficient) and which method requires
the fewest (the most efficient). Explain!
2^(ceil(log2(n)))
Filters
In the past, we asked that you compare the efficiency of your new MATLAB
function with the conv function by plotting the flops of each verses N , where
N is the length of both x(n) and h(n). However, the flops command that
measures the number of flops is not available in recent versions of MATLAB.
ZEROS OF H(Z)
1
Imaginary Part
2
1.5
1
0.5
0
1
0.5
4
0
0.5
1
0.5
0
0.5
Real Part
y = Ax
(b) As part of an algorithm you need to compute the matrix vector product
g = A x where A is the following kind of matrix,
A=
b(2) b(1) a(0) a(1) a(2)
b(3) b(2) b(1) a(0) a(1)
b(4) b(3) b(2) b(1) a(0)
3/4
/2
/4
/4
/2
3/4
By extending A appropriately, you can compute the vector g by using circular convolution and the FFT. Describe how to extend A and x and how
to use the FFT to efficiently compute the matrix-vector product A x.
(The FFT method is more efficient for large N than direct matrix vector
multiplication, because it takes O(N log(N )) instead of O(N 2 ).)
(d) how would you classify the new filter (low-pass filter, high-pass, bandpass, band-stop, etc)?
19
ZEROS OF H(Z)
1
Imaginary Part
2
1.5
1
0.5
0
0.5
(b) What is the relationship between the zeros of G(z) and the zeros of
H(z)?
Based on the zero-diagram of H(z) shown below, sketch the zerodiagram of G(z).
0.5
1
0.5
0
0.5
Real Part
ZEROS OF H(Z)
1.5
0.4
0.3
0.2
3/4
/2
/4
/4
/2
3/4
1
Imaginary Part
0.1
0
1
0.5
3
0
0.5
1
1.5
0
Real Part
1
0.8
0.6
0.4
0.2
0
20
3/4
/2
/4
/4
/2
3/4
(c) What is the relationship between the zeros of G(z) and the zeros of
H(z)? Based on the zero-diagram of H(z) shown below, sketch the
zero-diagram of G(z).
(c) Find an expression for the frequency response Gf () of the new filter
g(n) in terms of the frequency response H f () of the filter h(n). Sketch
the frequency response magnitude |Gf ()| of the new filter.
(d) What is the relationship between the zeros of G(z) and the zeros of
H(z)? Based on the zero-diagram of H(z) shown below, sketch the
zero-diagram of G(z).
(e) How would you classify the new filter (low-pass filter, high-pass, bandpass, band-stop, etc)?
ZEROS OF H(Z)
1.5
0.4
0.2
0.1
0
1
0.5
0.5
0.4
1.5
1
0.6
0.5
1
0
Real Part
0.3
0.2
0.5
5
0
0.5
0.1
0
1
0
0.8
0.5
0
0.5
Real Part
0.6
2.5
0.4
0.2
1.5
ZEROS OF H(Z)
0.7
Imaginary Part
Imaginary Part
0.3
3/4
/2
/4
/4
/2
3/4
0.5
0
21
3/4
/2
/4
/4
/2
3/4
3.6 Very simple low-pass digital filtering can be performed by a running average
of N consecutive samples,
y(n) =
3.8 For the discrete-time LTI system having the impulse response h(n)
(
0.5 0 n 5,
h(n) =
0
otherwise.
N 1
1 X
x(n k)
N
k=0
where x(n) is the signal being smoothed by the filter. In the following,
assume N = 8.
(a) Sketch the impulse response of this low-pass filter.
(b) Sketch the frequency response magnitude |H(ej )|. Accurately indicate the values at = 0 and = . Also accurately indicate the nulls
of the frequency response.
(b) What is the dc gain of the filter? That means, find H f (0).
(c) Find the zeros of the transfer function and sketch the pole/zero diagram.
(d) Based on the zero diagram in (c), sketch the frequency response amplitude A() and the frequency response magnitude |H()|.
k=0
3.7 Very simple low-pass digital filtering can be performed by a running average
of N consecutive samples,
y(n) =
N 1
1
1 X
y(n 1) +
x(n k)
5
N
N 1
1 X
x(n k)
N
(b) Find the zeros of the transfer function and sketch the pole/zero diagram.
k=0
where x(n) is the signal being smoothed by the filter. In the following,
assume N = 7.
(c) Based on the zero diagram in (c), roughly sketch the frequency response ampltude A(). Explicitly indicate any frequency response
nulls if there are any.
(c) Find the zeros of the transfer function and sketch the pole/zero diagram.
(d) Based on the zero diagram in (c), sketch the frequency response ampltude A(). (This is a linear-phase FIR filter.)
N 1
1 X
x(n k)
N
k=0
(e) A simple filter to remove the dc component of a signal can be implemented using the equation
N 1
1 X
y(n) = x(M )
x(n k)
N
(b) Find the zeros of the transfer function and sketch the pole/zero diagram.
k=0
(c) Based on the zero diagram, roughly sketch the frequency response
ampltude A(). Explicitly indicate any frequency response nulls if
there are any.
22
3.11 Very simple high-pass digital filtering can be performed by the difference
equation:
y(n) =
(b) Based on the pole-zero diagram, roughly sketch the frequency response
magnitude |H()|. Explicitly indicate any frequency response nulls if
there are any.
N 1
1 X
(1)k x(n k)
N
k=0
3.14 The transfer function H(z) of an FIR filter has zeros in the z-plane as
illustrated.
WN
(d) Based on the zero diagram in (c), sketch the frequency response ampltude A() (this is a linear-phase FIR filter) and the frequency response magnitude |H f ()|.
N = 12
The zeros on the unit circle are at powers of W12 . The dc-gain of the filter
is 2. All the poles are at z = 0.
N 1
1 X
v(n) =
x(n k)
N
(b) Sketch the frequency response magnitude |H(ej )|. Accurately indicate the nulls of the frequency response.
k=0
3.15 The transfer function H(z) of an FIR filter has 13 zeros in the z-plane as
illustrated.
exp(j 6 )
(c) Find the zeros of the transfer function and accurately sketch the
pole/zero diagram.
0.5
2.0
(d) Based on the zero diagram in (c), sketch the frequency response magnitude |H()|.
The zeros on the unit circle are at powers of W12 . The dc-gain of the filter
is unity. All the poles are at z = 0.
23
3.17 Matching. The diagrams on the following three pages show the impulse
responses, pole-zero diagrams, and frequency responses magnitudes of 8
discrete-time causal LTI systems. But the diagrams are out of order. Match
each diagram by filling out the following table.
3.16 The transfer function H(z) of an FIR filter has 13 zeros in the z-plane as
illustrated.
ej 6
2.0
0.5
The zeros on the unit circle are at powers of W12 . The dc-gain of the filter
is unity. All the poles are at z = 0.
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.5
0.5
0.5
0.5
1
0
10
15
20
25
IMPULSE RESPONSE 3
15
20
25
20
25
20
25
20
25
IMPULSE RESPONSE 4
0.5
0.5
0.5
0.5
1
0
10
15
20
25
IMPULSE RESPONSE 5
10
15
IMPULSE RESPONSE 6
0.5
0.5
0.5
0.5
1
0
10
15
20
25
IMPULSE RESPONSE 7
10
15
IMPULSE RESPONSE 8
0.5
0.5
0.5
0.5
1
0
24
10
10
15
20
25
10
15
POLEZERO DIAGRAM 2
1
0.5
0.5
imag part
imag part
POLEZERO DIAGRAM 1
1
0.5
FREQUENCY RESPONSE 1
0.5
FREQUENCY RESPONSE 2
1.4
1.2
0.8
0.6
0.4
0.2
1
1
1.5
0.5
0
real part
0.5
1.5
1.5
POLEZERO DIAGRAM 3
0.5
0
real part
0.5
0
0.5
1.5
POLEZERO DIAGRAM 4
0.5
0.5
1
0
FREQUENCY (CYCLES/SECOND)
0.5
0
0.5
FREQUENCY RESPONSE 3
3.5
3
2.5
imag part
imag part
2
0
3
1.5
2
0.5
0.5
1
1
1.5
0.5
0
real part
0.5
0.5
FREQUENCY RESPONSE 4
6
5
0
FREQUENCY (CYCLES/SECOND)
1.5
1.5
POLEZERO DIAGRAM 5
0.5
0
real part
0.5
0
0.5
1.5
POLEZERO DIAGRAM 6
0.5
0.5
0.5
0
FREQUENCY (CYCLES/SECOND)
0.5
0
0.5
FREQUENCY RESPONSE 5
0
FREQUENCY (CYCLES/SECOND)
0.5
FREQUENCY RESPONSE 6
3.5
20
3
15
imag part
imag part
2.5
0
10
1.5
0.5
0.5
5
1
1.5
0.5
0
real part
0.5
1.5
1.5
POLEZERO DIAGRAM 7
0.5
0
real part
0.5
0.5
0.5
1.5
POLEZERO DIAGRAM 8
0.5
0.5
0
FREQUENCY (CYCLES/SECOND)
0.5
0
0.5
FREQUENCY RESPONSE 7
imag part
imag part
6
5
4
2
3
1.5
0.5
0.5
1
1
1.5
0.5
0
real part
0.5
1.5
0.5
FREQUENCY RESPONSE 8
3.5
2.5
0
FREQUENCY (CYCLES/SECOND)
1.5
0.5
0
real part
0.5
0.5
0.5
1.5
1
0
FREQUENCY (CYCLES/SECOND)
Impulse response
1
2
3
4
5
6
7
8
25
0.5
0
0.5
Pole-zero
0
FREQUENCY (CYCLES/SECOND)
Frequency response
0.5
FREQUENCY RESPONSE 1
FREQUENCY RESPONSE 2
1.5
1.5
0.5
0.5
0
0
Impulse Response
Frequency Response
Zero Diagram
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 3
0.4
0.6
0.8
FREQUENCY RESPONSE 4
1
2
3
4
5
6
0.2
2.5
2
1.5
1.5
1
1
0.5
0.5
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 5
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 6
2.5
1.5
1.5
1
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.5
0.5
0.5
0.5
0
0
0.5
0.5
1
2
IMPULSE RESPONSE 3
IMPULSE RESPONSE 4
0.5
0.5
0.5
0.5
1
0
IMPULSE RESPONSE 5
IMPULSE RESPONSE 6
0.5
0.5
0.5
0.5
1
0
0.4
0.6
/
1
0
0.2
26
0.8
0.2
0.4
0.6
/
0.8
ZERO DIAGRAM 1
Imaginary Part
Imaginary Part
0.5
5
0
0.5
3.21 Design a simple causal real discrete-time LTI system with the properties:
0.5
(b) The system should annihilate constant signals. That is, the frequency
response should have a null at dc.
0.5
0.5
0
0.5
Real Part
ZERO DIAGRAM 3
0.5
0
0.5
Real Part
ZERO DIAGRAM 4
0.5
5
0.5
Imaginary Part
1
Imaginary Part
(a) The system should exactly preserve the signal cos(0.5 n).
0.5
0.5
1
1
0.5
0
0.5
Real Part
ZERO DIAGRAM 5
0.5
0
0.5
Real Part
ZERO DIAGRAM 6
3.22 Transfer functions. The transfer function H(z) of a causal FIR filter has
zeros in the z-plane as illustrated.
0.5
5
0.5
1
Imaginary Part
1
Imaginary Part
ZERO DIAGRAM 2
0.5
5
2.0j
0.5
1
0.5
0
0.5
Real Part
0.5
0
0.5
Real Part
0.5j
H(z) = (1 z 1 )3 (1 + z 1 )3
(2)
(2)
0.5j
2.0j
(a) The filter completely rejects the frequencies 5 Hz and 2.5 Hz.
(b) Accurately sketch the impulse response of the filter (as far as it can
be determined).
(c) Sketch the frequency response magnitude |H(ej )|. Accurately indicate the nulls of the frequency response.
(a) Sketch the real and imaginary parts of the impulse response of the
new filter, g(n).
(b) What is the relationship between the zeros of G(z) and the zeros of
H(z)? Based on the zero-diagram of H(z) shown below, sketch the
zero-diagram of G(z).
(c) Find an expression for the frequency response Gf () in terms of the
frequency response H f (). Sketch |Gf ()|.
IMPULSE RESPONSE h(n)
ZEROS OF H(Z)
1.5
0.4
Imaginary Part
0.3
0.2
0.1
0
1
0.5
3
0
0.5
1
1.5
0
Real Part
1
0.8
0.6
0.4
0.2
0
28
3/4
/2
/4
/4
/2
3/4
Imaginary Part
1
0.5
What would you choose for the desired frequency response D() and weighting function W ()? Sketch those functions for 0 where is
normalized frequency (radians/sample).
0
0.5
1
1
0.5
0
0.5
Real Part
(b) Find and sketch the impulse response h(n). (The zeros of H(z) are
evenly spaced on the unit circle.)
(c) Define F (z) = H(z). Sketch the zero-diagram of F (z) and the impulse response f (n). Sketch the frequency response magnitude of the
filter, |F (ej )|.
(d) Define S(z) = H(z 2 ). Sketch the zero-diagram of S(z) and the impulse
response s(n). Sketch the frequency response magnitude of the filter,
|S(ej )|.
(e) Define G(z) := H(z) H(1/z). Sketch the zero-diagram of G(z) and
the impulse response G(n). Sketch the frequency response magnitude
of the filter, |G(ej )|.
3.30 Filter specifications. The spectrum of an analog bandpass signal is contained between 50 Hz and 80 Hz. The signal is corrupted by low-frequency
noise; the noise is bandlimited to 30 Hz. The noisy analog signal is sampled
at 200 samples/second. A digital filter is to be designed so as to remove
the noise from the signal. What would you choose for the desired frequency
response D() and weighting function W ()? Sketch those functions for
0 where is normalized frequency (radians/sample).
29
4.1 (Porat 9.4) The impulse response of a causal linear-phase FIR filter starts
with the values
h(0) = 1,
h(1) = 3,
h(2) = 2.
Find the shortest FIR impulse response h(n) for each of the four types.
(h(3) =?, h(4) =?, etc.)
4.2 (Based on Porat 9.7) When two linear-phase FIR filters h1 (n) and h2 (n)
are connected in cascade as shown the total filter linear-phase.
x(n)
What would you choose for the desired frequency response D() and weighting function W ()? Sketch those functions for 0 where is
normalized frequency (radians/sample). In normalized frequency, at what
frequency o should the null be located? (0 o ). Provide an explanation for your answer.
- H1 (z)
- H2 (z)
y(n)
What is the filter type of the total filter as a function of the types of h1 (n)
and h2 (n)? (Fill in the table.)
I
II
III
IV
II
III
IV
III
IV
IV
?
?
?
?
?
?
?
?
?
?
If two linear-phase FIR filters h1 (n) and h2 (n) are connected in parallel as
shown, is the total filter linear-phase?
30
4.7 (a) List all of the four types of linear-phase FIR filters can be used for the
implementation of a bandpass filter.
- H1 (z)
x(n)
- H2 (z)
(b) List all of the four types of linear-phase FIR filters can be used for the
implementation of a bandstop filter.
?
+j - y(n)
6
4.3 (Porat 9.5) Let H(z) be the transfer function of a linear-phase FIR filter
with real coefficients. The filter is known to have zeros in the following
locations:
n
o
z = 1, 0.5 ej/3 , 5, j
(a) Because H(z) is a linear-phase FIR filter it must have other zeros as
well, because the zeros of a linear-phase FIR must exist in particular
configurations. What are the other zeros of H(z)? Sketch the zero
diagram of this filter.
(b) What is the minimal length of the filter?
(c) What is the filter type of the minimal length filter? (I,II,III, or IV)
4.4 Consider a Type-I linear-phase FIR filter H(z) (with real-valued impulse
response). You are told that H(z) has zeros at z = 1 and z = 2 + 2j.
(a) What other zeros must H(z) have? Sketch the zero-diagram.
(b) Find the impulse response, h(n).
4.5 For the transfer function
H(z) = z 2 z 6
of an FIR linear-phase filter,
(a)
(b)
(c)
(d)
31
Test for Type III using a lowpass differentiator obtained by convolving the impulse response of a Type II lowpass filter and Type IV
differentiator.
h3 = conv(h2,h4);
where amplitude response A() is continuous and real-valued. Write a
Matlab function that uses the FFT to compute A() from the impulse
response h(n) at the frequency points
k =
2
k
L
0 k L 1.
4.9 Filter transformations. Let h(n) be the Type I linear-phase lowpass FIR
filter obtained in the previous problem. Consider the three transformations:
A) G(z) = H(z)
B)
G(z) = z M H(z)
C) G(z) = z M H(z)
Lets denote the amplitude response of h(n) by Ah () and the amplitude
response of g(n) by Ag (). For each of the filter transformations,
(a) Find an expression for g(n) in terms of h(n).
(b) Find an expression for Ag () in terms of Ah ().
Verify your Matlab function works correctly for each of the 4 filter types
by using it to plot A() for the following filters.
(c) Plot g(n), Ag () and use zplane to plot the zeros of G(z). (In the
command zplane(g), the impulse response g should be a row vector.
If it is a column vector you can type zplane(g).) Use the Matlab
command subplot to save paper. Also use orient tall so that when
printed out, the plots use the whole page.
(d) If H(z) is a low-pass filter with cutoff frequency o , what kind of filter
is G(z) and what is its cutoff frequency?
4.10 Matching. Match each impulse response with its frequency response and
zero diagram. You should do this problem with out using a computer.
32
IMPULSE RESPONSE A
FREQ RESP A
IMPULSE RESPONSE B
2
0
FREQ RESP B
2.5
2.5
1.5
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
/
IMPULSE RESPONSE C
FREQ RESP C
IMPULSE RESPONSE D
2
0
2.5
1.5
1.5
0.5
0.5
0.2
0.4
0.6
0.8
Impulse Response
0.5
0.5
0
0.5
1
0
0.5
1
0.5
0
0.5
Real Part
0.5
ZPLANE C
1
0.5
0.5
0
0.5
Real Part
ZPLANE D
Imaginary Part
Imaginary Part
A
B
C
D
ZPLANE B
1
Imaginary Part
Imaginary Part
ZPLANE A
1
0.5
1
0
0.5
1
0.5
0.5
0.5
0.5
33
0.8
0
0
0.8
FREQ RESP D
2.5
0.6
/
0.2
0.4
0.6
/
Zero Diagram
Frequency Response
4.11 Matching. Match each impulse response with its frequency response and
zero diagram by filling out the following table. You should do this problem
with out using a computer.
IMPULSE RESPONSE B
ZPLANE B
1
0.5
0.5
Imaginary Part
Imaginary Part
IMPULSE RESPONSE A
ZPLANE A
1
0
0.5
1
0.5
0
0.5
Real Part
0.5
ZPLANE C
0
0.5
Real Part
ZPLANE D
IMPULSE RESPONSE C
0.5
0.5
IMPULSE RESPONSE D
2
Imaginary Part
0
0.5
1
0.5
2
0
FREQUENCY RESPONSE A
A
B
C
D
FREQUENCY RESPONSE B
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE C
FREQUENCY RESPONSE D
0
0.2
0.4
0.6
/
0.8
0.2
0.4
0.6
0.8
34
0
0.5
1
0
0.5
Real Part
Impulse Response
2
2
0
0.5
1
Imaginary Part
Zero Diagram
0.5
0
0.5
Real Part
Frequency Response
FREQUENCY RESPONSE 1
FREQUENCY RESPONSE 2
1.5
1.5
0.5
0.5
0
0
Impulse Response
Frequency Response
Zero Diagram
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 3
0.4
0.6
0.8
FREQUENCY RESPONSE 4
1
2
3
4
5
6
0.2
2.5
2
1.5
1.5
1
1
0.5
0.5
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 5
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 6
2.5
1.5
1.5
1
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.5
0.5
0.5
0.5
0
0
0.5
0.5
1
2
IMPULSE RESPONSE 3
IMPULSE RESPONSE 4
0.5
0.5
0.5
0.5
1
0
IMPULSE RESPONSE 5
IMPULSE RESPONSE 6
0.5
0.5
0.5
0.5
1
0
0.4
0.6
/
1
0
0.2
35
0.8
0.2
0.4
0.6
/
0.8
ZERO DIAGRAM 1
Imaginary Part
Imaginary Part
0.5
5
0
0.5
H(z) = (1 z 1 )2 (1 + z 1 )2
0.5
5
0
0.5
1
1
1
0
Real Part
ZERO DIAGRAM 3
0.5
5
0.5
Imaginary Part
Imaginary Part
0
Real Part
ZERO DIAGRAM 4
0.5
5
4.14 Show a derivation for each of your answers. Do not use MATLAB to answer
any part of this question.
0.5
1
0
Real Part
ZERO DIAGRAM 5
0
Real Part
ZERO DIAGRAM 6
H1 (z) = (1 + z 1 )3
0.5
5
0.5
1
Imaginary Part
1
Imaginary Part
(d)
ZERO DIAGRAM 2
0.5
0.5
H2 (z) = H1 (z) (1 z 1 )2
1
1
0
Real Part
0
Real Part
For each of the following linear-phase filters: sketch the zero diagram, impulse response, amplitude response, and classify each filter as type I, II,
III, IV. Comment on your observations, especially with regard to what you
know about the zero locations of linear-phase filters.
(a)
(b)
H(z) = (1 z 1 )2 (1 + z 1 )
(c)
H(z) = (1 z 1 )(1 + z 1 )2
36
4.15 Design by DFT-based interpolation. Implement the DFT-based interpolation approach for Type II. (Modify the code in the notes for the Type I
case.) In particular, consider the design of a length N = 12 low-pass Type
II FIR filter where the amplitude function A() interpolates 1 and 0 as
follows.
2
1 k = 0, 1, 2
k =
A
0 k = 3, 4, 5, 6
N
-0.0318
-0.0494
0.0891
-0.0255
-0.1287
0.2892
0.6429
0.6429
0.2892
-0.1287
-0.0255
0, for other [, ]
0.0891
-0.0494
-0.0318
4.17 Optional : Mitra 4.19. (But use the frequencies 0.2, 0.4, 0.9)
4.18 Let h(n) be a (real-valued) FIR impulse response, not necessarily linearphase. Then the impulse response defined by convolving h(n) with its
time-reversed version is always a linear-phase FIR filter. In other words,
p(n) := h(n) h(n), is a linear-phase FIR impulse response. The question
is: which of the four linear-phase FIR filter types can p(n) be?
(d) The ideal band-pass Hilbert transformer (3) can not be realized. Consider the design of a real linear-phase FIR filter to approximate it.
Which of the four basic filters types are appropriate? Explain.
(e) Use the impulse response truncation method to design a real linearphase FIR filter approximating the band-pass Hilbert transform (3).
Provide an explicit formula for the impulse response h(n). Simplify so
that j does not appear. Which of the four types is your filter?
4.22 (a) Show that the transfer function H(z) of a Type I FIR filter will always
have a zero of even multiplicity at z = 1. One way to do this is by
contradiction. Assume the zero if of odd multiplicity: write H(z) as
h(n) = h(N 1 n)
can be written as
(2)
H(z) = H(1/z) z (N 1) .
|| <
(b) Suppose h(n) is a Type II FIR filter. We learned that the transfer
function H(z) must have a zero at z = 1. That is: H(1) = 0.
Classify the following statement as true or false. H(z) always has an
odd number of zeros at z = 1.
If true, give a derivation. If false, give an example of a Type II FIR
filter with an even number of zeros at z = 1.
(a) Which of the four basic types of FIR filter (I, II, III, IV) have A() =
A()? Which types have A() = A()? Show a clear derivation.
(b) The linear-phase filter (2) is cascaded with an ideal delay system D()
with delay samples, R.
x(n) H() D() y(n)
4.23 An FIR filter with a complex-valued impulse response can have linear phase.
Find the frequency response of each of the two FIR filters with the impulse
responses:
Does the total system have linear-phase? If so, what is the amplitude
response of the total system in terms of A()? Show a clear derivation.
38
h1 = [1 j 2 j 1];
h2 = [1 j 2 -j 1];
5.1 Linear-phase FIR filter design using windows (Based on Mitra M7.18,
M7.19)
One of these two FIR filters has linear phase. Which one? What is its the
phase response?
Use each of the Hamming, Hann, Blackman, and Kaiser windows to design
four linear-phase FIR digital lowpass filters. Each of the four filters should
be of the same length. It is desired that the filters meet the following specifications: passband edge fp at 2 Hz, stopband edge fs at 4 Hz, maximum
passband attenuation of 0.1 dB, minimum stopband attenuation of 40 dB.
The filter is to operate at a sampling frequency Fs of 20 Hz.
2-j
2+2j
2-2j
2+j
1-j];
h2 = [1+j
2-j
2+2j
2+2j
2-j
1+j];
h3 = [1+j
2-j
2+2j
h4 = [1+j
2-j
2+2j
5 2+2j
h5 = [1+j
2-j
2+2j
5j
2-2j
2+j
2-2j
First, determine the length of the filter using the MATLAB command
kaiserord and the specifications listed above. Use the length provided
by kaiserord for each of the four filters. You should find that the Kaiser
window leads to a filter that meets the specifications, but that the other
windows lead to filters that do not quite meet the specification.
1-j];
2-j
1+j];
2+j
Display the amplitude response and magnitude response in dB. You will use
the commands sinc, hamming, hanning, blackman, and kaiser. (Dont use
the fir1 command, although you may use it to check your work.)
1-j];
For those filters that you claim have linear phase, can you show why they
have linear phase?
4.25 Consider a Type I FIR filter of length N . The impulse response h(n) is
non-zero for 0 n N 1. We have learned how to compute the realvalued amplitude response A() on a uniform grid. Hans has suggested
another way. He defines the length L signal g(n)
h(n + M )
0nM
g(n) =
0
M +1nL1
where M = (N 1)/2 as usual. (This is simply the second half of h(n)
including the midpoint.) Then he computes:
C(k) = 2 Real{DFTL {g(n)}} h(M )
for 0 k L 1. Hans claims that
2
k .
C(k) = A
L
Is he correct? If so, derive the correctness of his formula. If not, show why
not.
To make it clear, suppose
h = [1, 2, 4, 2, 1]
Then in Matlab notation, Hans claims that
C = 2*real(fft([4 2 1 zeros(1,L-3)])) - 4;
will compute L samples of A() for =
2
L k,
Windows
0 k L 1.
39
5.2 Windows. Below are shown three different window functions, and the
discrete-time Fourier transform of each. Match each window function with
its Fourier transform.
W1()
Window function A
where d(n) is the inverse DTFT of the ideal low-pass frequency response.
The figures below illustrate two different sinc functions d(n) and 3 different
window functions w(n). In each case identify the frequency response of the
resulting filter h(n) by filling the table below.
30
1
0.8
20
0.6
0.4
10
0.2
0
10
20
30
0
1
0.5
0.5
0.5
Sinc Function
Window
1
1
1
2
2
2
1
2
3
1
2
3
W2()
Window function B
30
1
0.8
20
0.6
0.4
Frequency Response
10
0.2
0
10
20
30
0
1
SINC FUNCTION 1
0.5
0
W3()
Window function C
30
SINC FUNCTION 2
0.6
0.6
0.4
0.4
0.2
0.2
1
0.8
20
0.6
0.4
0.2
10
0.2
0
10
20
30
40
10
20
30
40
0.2
0
10
20
n
30
0
1
0.5
0
/
0.5
1
WINDOW 1
WINDOW 2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
10
20
30
40
30
40
WINDOW 3
1
0.8
0.6
0.4
0.2
0
0
40
10
20
10
20
30
40
FREQUENCY RESPONSE 1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
Normalized Frequency
FREQUENCY RESPONSE 3
1
1
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0.2
0.4
0.6
0.8
Normalized Frequency
0.4
0.6
0.8
Normalized Frequency
where w(n) is the window function and d(n) is the inverse DTFT of the
ideal frequency response. The figures below show the impulses responses
and frequency responses of six filters designed using the window method.
But they are out of order. Match each frequency response to its impulse
response by completing the table.
Impulse Response
FREQUENCY RESPONSE 5
0.2
0.4
0.6
0.8
Normalized Frequency
1
0.8
0.6
0.6
0.4
0.4
0.2
IMPULSE RESPONSE 1
0.2
0
0.2
0.4
0.6
0.8
Normalized Frequency
Frequency Response
1
2
3
4
5
6
FREQUENCY RESPONSE 6
0.8
0.2
FREQUENCY RESPONSE 4
0.8
5.4 Linear-Phase FIR Filter Design by Windows. The design of a linearphase FIR filter using the window method gives an impulse response of the
form
FREQUENCY RESPONSE 2
0.2
0.4
0.6
0.8
Normalized Frequency
IMPULSE RESPONSE 2
0.2
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0.05
0.05
0
10
20
30
40
IMPULSE RESPONSE 3
0.6
0.4
0.4
0.2
0.2
0
10
20
30
40
IMPULSE RESPONSE 5
30
40
10
20
30
40
IMPULSE RESPONSE 6
0.2
0.2
0.1
0.1
0.1
0.1
0.2
0.2
0
41
20
IMPULSE RESPONSE 4
0.6
10
10
20
30
40
10
20
30
40
FREQUENCY RESPONSE 1
FREQUENCY RESPONSE 2
0.5
where the frequencies fi are unknown. On the class webpage, a 100 point
signal (signal2.txt) of this form is available. Using appropriate windows,
estimate the three frequencies fi .
0
0
0.25
0.5
0.75
FREQUENCY RESPONSE 3
0.5
0.5
0.25
0.5
0.75
FREQUENCY RESPONSE 4
Suppose you want to use windowing and the FFT to estimate the frequencies of the sinusoidal components from the data (you do not know the
frequencies). For which of the following two signals is the estimation of the
frequencies more difficult? (MATLAB is not required for this problem.)
0
0
0.25
0.5
0.75
FREQUENCY RESPONSE 5
0.5
0.5
0.25
0.5
0.75
FREQUENCY RESPONSE 6
0
0
0.25
0.5
0.75
0.25
0.5
0.75
5.5 (Porat 6.1) We saw that the height of the largest side-lobe of the rectangular
window is about -13.5 dB relative to the main-lobe. What is the relative
height of the smallest side-lobe? You may assume the window length N is
odd. (Hint: At what frequency will it be located?).
(a) Explain why a rectangular window is not adequate for detecting the
second component. (Hint: the previous problem is relevant.)
(b) Of the Hann and Hamming windows, guess which one is better in this
case for detecting the second component, then test your guess on a
computer.
42
Suppose you want to use windowing and the FFT to estimate the frequencies of the sinusoidal components from the data (you do not know the
frequencies). For which of the following two signals is the estimation of the
frequencies more difficult? (MATLAB is not required for this problem.)
which can be written as f (n) = p(n) x(n). Sketch the DTFT, |F ()|,
for || . Show your work. Indicate the nulls of F () on your sketch.
At what frequencies (approximately) does |F ()| have its maximum
value?
(d) Similarly, consider the discrete-time signal
(
cos(0.5 n) 0 n 9
v(n) =
0
otherwise
are observed, and the zero-padded FFT is graphed, would you see two
distinct peaks or not? Explain.
(f) Based on the previous parts of the problem: If 10 samples are observed
of a signal that consists of two cosine signals, what should be the
minimum difference between their frequencies so that one can easily
distinguish them on a graph of the zero-padded FFT?
43
1.5
5.14 STFT. (Mitra 11.16) For this problem, it is convenient to use the following
definition of the STFT, which is only superficially different from the form
in the lecture notes. Show that the sampled STFT defined by
X d (k, n) =
SIGNAL 1
R1
X
0.5
0
0.5
1
m=0
1.5
20
40
60
80
100
120
140
160
180
200
20
40
60
80
100
120
140
160
180
200
20
40
60
80
100
120
140
160
180
200
20
40
60
80
100
120
140
160
180
200
1.5
SIGNAL 2
1
0.5
0
0.5
1
1.5
1.5
1
SIGNAL 3
5.15 Spectrogram. Look at the Matlab documentation for the specgram function (use the help command). Make a spectrogram plot of the signal in
the data file signal1.txt that was used in a previous problem. Try it
with several different block lengths, amount of overlap, etc, until you get
a spectrogram that best reveals the structure of the signal. Compare and
contrast the spectrogram with the spectrum of this signal you got already.
0
0.5
1
1.5
1.5
Signal
1
SIGNAL 4
0.5
x(t)
0.5
0.5
0
0.5
1
1.5
0.5
0.01
0.02
0.03
0.04
time (sec)
0.05
0.06
0.07
5.16 STFT. The following figures show four signals and their spectrograms.
Match each signal to its spectrogram.
44
SPECTROGRAM B
0.5
0.4
0.4
frequency
frequency
SPECTROGRAM A
0.5
0.3
0.2
SIGNAL X1
SIGNAL X2
0.5
0.5
0.3
0.5
0.5
0.2
0.1
0.1
1
0
20
40
60
80
100
50
100
150
50
100
time
time
SPECTROGRAM C
SPECTROGRAM D
0.5
0.5
0.4
0.4
20
150
frequency
frequency
0.3
50
100
150
100
80
100
1
20
40
60
80
100
20
40
60
SIGNAL X6
0.5
0.5
2
0
time
50
100
1
0
150
time
5.17 The Spectrogram. Shown below are 6 signals and 6 spectrograms. Match
each signal to its spectrogram. Give an explanation of your answer.
Signal
80
0.5
0.2
0.1
100
0
1
SIGNAL X5
0.1
80
0.2
60
0.5
0.3
40
SIGNAL X4
1
0
SIGNAL X3
Spectrogram
1
2
3
4
5
6
45
20
40
60
80
100
20
40
60
0.4
0.4
0.1
0
20
40
time
0.2
where
0.1
R = block length.
60
0.5
0.4
0.4
0.3
0.2
(A) SPECTROGRAM, R = ?, L = ?
0.1
0
20
40
time
60
20
40
time
60
SPECTROGRAM 6
0.4
0.4
frequency
0.5
0.3
0.2
40
time
60
0.8
0.6
0.4
0.6
0.4
0.2
1000
2000
3000
time
4000
1000
2000
3000
time
4000
0.2
(C) SPECTROGRAM, R = ?, L = ?
0
20
40
time
60
frequency
20
0.8
0.3
0.1
0
0.2
(D) SPECTROGRAM, R = ?, L = ?
0.8
0.8
0.6
0.4
0.2
0
0.6
0.4
0.2
1000
2000
time
3000
4000
1000
2000
time
3000
4000
If you like, you may listen to this signal with the soundsc command; the
data is in the file: signal3.txt. Here is a figure of the signal.
3
2
1
x(n)
0.1
(B) SPECTROGRAM, R = ?, L = ?
1
0.2
SPECTROGRAM 5
frequency
For each of the four spectrograms, indicate what you think R and L are.
Briefly explain your choices.
0.3
0.5
60
frequency
0.1
0
40
time
SPECTROGRAM 4
0.5
frequency
frequency
SPECTROGRAM 3
20
L {35, 250}
frequency
0.2
0.3
frequency
0.3
5.18 STFT. Shown below are four spectrograms of the same signal. Each spectrogram is computed using a different set of parameters.
SPECTROGRAM 2
0.5
frequency
frequency
SPECTROGRAM 1
0.5
0
1
2
3
1000
2000
3000
4000
5000
n
46
6000
7000
8000
9000
5.20 Shown below are eight spectrograms of the same signal. (The signal is the
sum of a cosine and an impulse.)
0.5
0
SPECTROGRAM A
0.5
0
50
100
150
TIME (SAMPLES)
200
250
frequency
SPECTROGRAM B
0.5
0.5
0.4
0.4
frequency
SIGNAL
0.3
0.2
0.1
0
0
150
SPECTROGRAM D
0.2
0.1
50
100
150
SPECTROGRAM F
50
100
150
time
200
250
0.2
50
100
150
time
200
SPECTROGRAM D
0.4
0.1
0.3
frequency
50
100
time
SPECTROGRAM G
0.2
150
50
100
time
SPECTROGRAM H
0.2
0.3
0.2
0.1
50
0.2
100
150
50
time
100
time
0.1
0
0
50
100
150
time
200
250
150
0.4
0.3
0.4
0.1
0.4
250
0.5
0.2
0.3
0.5
0.5
0.3
150
0.1
frequency
frequency
0.1
100
0.4
0.2
0.1
50
SPECTROGRAM E
frequency
0.1
0.3
SPECTROGRAM C
frequency
frequency
0.4
time
0.3
0.4
frequency
frequency
SPECTROGRAM B
0.5
0.2
time
0.4
0.5
0.3
0.1
0
0.2
150
0.4
0.3
100
SPECTROGRAM C
0.3
A
B
C
D
SPECTROGRAM A
50
time
frequency
frequency
100
time
0.4
50
0.5
Spectrogram
0.2
0.1
R = block length
0.3
50
100
150
time
200
250
R {18, 45},
L {1, 10}
where
R = block length.
L = time lapse between blocks.
47
N {45, 256}
150
5.21 The Spectrogram. Shown below are eight spectrograms of the same
signal.
For each of the eight spectrograms, indicate what you think R, L, and N
are, by filling out the table. Explain!
SPECTROGRAM A
A
B
C
D
E
F
G
H
0.4
frequency
0.4
0.3
0.2
0.1
0
0.3
0.2
0.1
20
40
time
60
SPECTROGRAM C
20
40
time
60
SPECTROGRAM D
0.5
0.4
frequency
frequency
0.4
0.3
0.2
0.1
0.3
0.2
0.1
0
0
20
40
time
60
SPECTROGRAM E
20
40
time
60
SPECTROGRAM F
0.5
frequency
0.4
0.4
frequency
0.3
0.2
frequency
0.3
0.2
0.1
0.1
0
0
20
40
time
SPECTROGRAM G
60
0.5
0.5
0.4
0.4
frequency
frequency
Spectrogram
SPECTROGRAM B
0.3
0.2
0.1
20
40
time
SPECTROGRAM H
60
20
60
0.3
0.2
0.1
0
0
20
40
time
60
40
time
L {1, 5}
where
R = block length.
L = time lapse between blocks.
48
N {31, 128}
Spectrogram
For each of the eight spectrograms, indicate what you think R, L, and N
are, by filling out the table. Explain!
Spectrogram
A
B
C
D
A
B
C
D
E
F
G
H
SPECTROGRAM B
0.5
0.4
0.4
Frequency (/2)
Frequency (/2)
SPECTROGRAM A
0.5
0.3
0.2
0.1
0
0.3
0.2
0.1
100
200
Time (sample)
300
SPECTROGRAM C
0.2
0.1
100
200
Time (sample)
300
SPECTROGRAM D
0.5
0.5
0.4
0.4
Frequency (/2)
Frequency (/2)
5.22 The following figures shows a digitized echolocation pulse emitted by the
Large Brown Bat, Eptesicus Fuscus. There are 400 samples; the sampling
period was 7 microseconds.
0.3
0.2
0.1
SIGNAL
0.3
0.2
0.1
0
0
0
0
0.1
100
200
Time (sample)
300
100
200
Time (sample)
300
0.2
0.3
0
50
100
150
200
SAMPLE
250
300
350
400
The spectrogram of the bat pulse was computed with parameter values
having a subset of the following values.
R {22, 50},
L {2, 20},
N {50, 256}
where
R = block length
L = time lapse between segments.
N = FFT length (nfft in the MATLAB specgram function). (Each block
is zero-padded to length N .)
For each of the spectrograms, indicate what you think R, L, and N are by
filling out the table, and explain your choices.
frequency ((radians/(2))/second)
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0
100
200
300
400
500
Which of the following three parameters would you suggest the student
modify to best improve the appearance of this spectrogram? What change
would you make to that parameter? The three parameters are:
R = block length.
L = time lapse between blocks.
N = FFT length. (Each block is zero-padded to length N .)
Explain your answer.
5.24 The Spectrogram. The signal x(n) = cos(o n) has the frequency o
radians/second. If o changes with time, then the signal has a time-varying
frequency. Here o = 0.01 n. It is just increasing. But the spectrogram
does not look like it is increasing it shows that the frequency increases,
then decreases, then increases, then decreases. Explain why.
5.25 Match each of the following Matlab commands to the figure it produces.
50
Figure 1
Figure 2
10
30
6.1 Redo the example of the weighted low-pass filter from the lecture notes,
but increase the stop-band weighting K from 10 to 40.
20
6
4
10
2
0
0.5
0.5
0
0.5
Figure 3
20
15
15
10
10
0.5
Figure 4
20
0
0.5
0.5
0
0.5
Figure 5
0.5
Figure 6
20
10
8
15
6.2 In the first design example of the notes we plotted the amplitude response
for N = 31, 61, 121, but we did not do this for the later design examples.
Repeat with longer lengths the design examples of the spline transition
low-pass filter and the zero-weighted transition-band low-pass filter. (The
spline transition is the straight line connecting the pass-band and stopbands.) Compare with the impulse response truncation method. Please
comment on your observations.
6
10
4
5
0
0.5
2
0
0.5
0
0.5
Figure 7
0.5
Figure 8
10
20
15
6
10
4
5
2
0
0.5
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
0.5
0
0.5
0.5
N1 = 16;
N2 = 30;
N3 = 48;
stem([0:N1-1]/N1-0.5,abs(fftshift(fft(rectwin(N1),N1))),.)
stem([0:N1-1]/N1-0.5,abs(fftshift(fft(hamming(N1),N1))),.)
stem([0:N2-1]/N2-0.5,abs(fftshift(fft(rectwin(N1),N2))),.)
stem([0:N2-1]/N2-0.5,abs(fftshift(fft(hamming(N1),N2))),.)
stem([0:N3-1]/N3-0.5,abs(fftshift(fft(rectwin(N1),N3))),.)
stem([0:N3-1]/N3-0.5,abs(fftshift(fft(hamming(N1),N3))),.)
stem([0:N2-1]/N2-0.5,abs(fftshift(fft(rectwin(N2),N2))),.)
stem([0:N2-1]/N2-0.5,abs(fftshift(fft(hamming(N2),N2))),.)
6.3 What is the effect of the width of the zero-weighted transition band? As
the zero-weighted transition band becomes wider, what happens to the
frequency response? Show examples using Matlab to design several filters
to illustrate what happens.
51
6.4 Four FIR lowpass filters of equal length are designed according to different
weighting functions W (), as discussed in the lecture notes.
6.5 Repeat the design example from the lecture notes of a low-pass filter with
a specified null, so that the zero is (a) second order, and (b) third order.
(Impose the constraint A0 (1 ) = 0 and A00 (1 ) = 0 in addition to the
constraint A(1 ) = 0.) Plot the amplitude response A(), impulse response
h(n), and zero-diagram.
Match each frequency response with the weighting function used for the
design.
FREQUENCY RESPONSE A
WEIGHTING FUNCTION 1
6.6 The weighted square error design of band-pass filters can be done with the
following functions D() and W (),
0 < < a
(stop-band 1)
0
1
a < < b
(pass-band)
D() =
0
b < <
(stop-band 2)
0.8
0.6
0.5
0.4
0.2
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE B
0.2
0.4
0.6
0.8
and
WEIGHTING FUNCTION 2
K1
0
K2
W () =
K3
1.5
1
0.5
0.5
0
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE C
0.2
0.4
0.6
0.8
0 < < 1
1 < < 2
2 < < 3
3 < < 4
4 < <
WEIGHTING FUNCTION 3
0.8
D()
0.6
0.5
0.4
0.2
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE D
0.2
0.4
0.6
0.8
K1
WEIGHTING FUNCTION 4
K3
K2
20
15
W ()
10
0.5
5
0
0
0
0.2
0.4
0.6
/
0.8
0.2
0.4
0.6
0.8
Find the formula for q(k) and b(k). Use the formula to set up the appropriate matrix equation and design band-pass filters of lengths 31 and 61. Use
1 = 0.28, 2 = 0.32, 3 = 0.58, and 4 = 0.62 and K1 = 1, K2 = 3,
K3 = 10. Plot the amplitude response A(), impulse response h(n), and
zero-diagram.
52
6.7 The design of an FIR digital notch filter is carried out as in the lecture notes
using the least-square method. But now, a weighting function is used. The
three different weighting functions are used:
0 1
0 1
100
1
0
1 < < 2
0
1 < < 2
W1 () =
W2 () =
1
2
100
2
0 1
10
0
1 < < 2
W3 () =
10
2
6.8 Extend the weighted least square error approach to Type II filters. Using a
Type II filter instead of Type I filter, repeat the example on pages 18 and 19
of the notes, the weighted low-pass example: Design a length 32 symmetric
FIR filter with band-edges p = 0.26, s = 0.34 and a stop-band weight
of K = 10.
Plot the impulse response, frequency response, and zero-diagram of the
filter.
6.9 If you try to design a Type II notch filter, similar to the one described in
the notes, what problem would arise?
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
NOTCH FILTER 2
|Hf()|
|Hf()|
NOTCH FILTER 1
0.8
0.2
0.4
0.6
0.8
You want to design an FIR digital filter to remove the noise. If you are
using the least-squares design method, what should be the desired frequency
response D() and the weighting function W ()? ( ) Make a
sketch of both of these functions. Explain!
NOTCH FILTER 3
1
|Hf()|
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
6.11 Optional. We saw that truncating the inverse DTFT of D() yields a
filter that minimizes the unweighted (integral) square error. (With D()
being symmetric, we get a Type I filter.) Similarly, show that truncating
the inverse DFT of D(2k/L) yields a Type I filter that minimizes the
unweighted discrete square error, computed over a uniform grid:
Match each filter to the weighting function used for its design and give an
explanation for your answer.
2 =
Filter
Weighting Function
L1
X
k=0
1
2
3
with
k =
2
k,
L
0 k L 1.
6.12 Very optional. Extend the least square error approaches to Type III and
IV FIR filters.
53
6.13 Consider a Type I FIR filter of length N . The impulse response h(n) is
non-zero for 0 n N 1. We have learned how to compute the realvalued amplitude response A() on a uniform grid. Hans has suggested
another way. He defines the length L signal g(n)
h(n + M )
0nM
g(n) =
0
M +1nL1
p1 = 0.3 ,
p2 = 0.6 ,
s2 = 0.8 .
Note that the width of the two transition-bands are not the same. Find a
formula for the Type I FIR filter of length N that minimizes the unweighted
square error,
Z
E2 =
(A() D())2 d.
6.15 The Savitsky-Golay FIR filters (or windows) are a family of lowpass digital
filters obtained by minimizing an unweighted square error subject to side
constraints.
Is he correct? If so, derive the correctness of his formula. If not, show why
not.
To make it clear, suppose
(4)
h = [1, 2, 4, 2, 1]
A(0) = 1.
C = 2*real(fft([4 2 1 zeros(1,L-3)])) - 4;
will compute L samples of A() for =
2
L k,
(5)
0 k L 1.
1
0.8
0.6
0.4
6.14 Consider the design of a band-pass filter, where the desired amplitude response has the form
0.2
0
0.2
1
D()
0.8
0.6
0.4
0.2
0
/
0.2
0.4
0.6
0.8
(a) How would you solve the constrained minimization problem stated in
equations (4) and (5)? (Write down the equations you need to solve.)
Solve the equations and determine h(n).
0.8
0.6
0.4
0.2
(b) You can see that the stop-band ripples are rather large. By using
only an error weighting function how could you reduce the sizes of the
stopband ripples without changing the length of the filter? (Explain.)
What happens to the other characteristics of the frequency response?
0
0.2
0.1
0.2
0.3
0.4
0.5
/
0.6
0.7
0.8
0.9
54
7.1 Consider the design of a type I FIR filter of minimal length satisfying the
following specifications:
1 p A() 1 + p ,
s A() s ,
Kp
W ()
0 p
where
0
p = 0.018,
s = 0.004,
p = 0.31,
s = 0.38.
p s
What should Kp and Ks be? You can estimate what the length of the
impulse response should be with the command remezord. Use fircheb or
remez to design the filter. For these specifications, does remezord correctly
indicate what the shortest Type I filter length is? Plot the amplitude
response A() and impulse response h(n).
7.2 Use N = 81, with Kp = 1, Ks = 4, and find the low-pass filter that
minimizes the weighted Chebyshev error with band edges:
(d) If the length of the filter is increased to 81 and the filter is redesigned,
how are the new pass-band and stop-band ripples related to the original ones (with Kp = 1, Ks = 10)?
Plot A() and h(n). What happens to A() and h(n) as you make the
transition band narrower? Why might one not want to use the last filter in
practice?
h1 (n) = [1, 2, 2, 1]
has very simple coefficients (is simple to implement) and has zeros exactly
on the unit circle that fall in the stop-band of a low-pass filter you are
designing. You will use it in cascade with another filter h2 (n) to reduce the
total number of multiplications. Design h2 (n) using the Remez algorithm
so that the total response
7.3 A low-pass length 51 Type I FIR filter is designed using the Remez (ParksMcClellan) algorithm with the following desired amplitude response D()
and weighting function W () with Kp = 1, Ks = 10. The filter has a
pass-band ripple of 0.015 (about 0.03 from peak to peak).
A() = A1 () A2 ()
of h(n) = h1 (n) h2 (n) minimizes the weighted Chebyshev error. The
length of h(n) should be 22. Use the following band edges and weightings.
p = 0.2, s = 0.3, and Kp = 1, Ks = 3. Plot A(), h(n) and the zero
diagram.
1
D()
You may want to use fircheb rather than remez because you will need to
use a non-uniform weighting function W () and desired response D().
55
7.5 A lowpass filter H(z) designed using the Remez algorithm will not in general
have the property that the dc gain is 1 exactly. (That is, H f ( = 0) is not
usually exactly 1.) However, suppose you want to design a filter according
to the minimax criterion that does have unity dc gain (H f ( = 0) = 1).
How could you use the Remez algorithm to design a filter satisfying this
constraint?
7.8 The following type of linear-phase FIR filter will be useful later on for signal
interpolation.
h(n) = [h0 , h1 , 0, h3 , h4 , 1/3, h4 , h3 , 0, h1 , h0 ]
This kind of filter is called a Nyquist-3 filter.
Use the linear programming approach to design a length-11 lowpass digital of this form so as to minimize the maximum error. Use the following
parameters,
7.6 Modify the program fircheb so that it can design type II FIR filters that
minimize the weighted Chebyshev error. The size of the reference set (R)
should be equal to the number of filter parameters plus 1. Use your program
to design a length 32 low-pass filter with p = 0.26, s = 0.34 and
Kp = 1, Ks = 4.
Kp = Ks = 1,
o = /3,
p = o 0.1 ,
s = o + 0.1 .
7.7 Linear Programming Method In the second example of linear programming, we had an example where constraints were imposed on the step response of the filter. Notice that the step response s(n) will equal a constant
after a certain point.
s(n) = c,
for all
n N 1.
Show that this constant c is given by c = A(0). Add the constraint A(0) = 1
to the design problem and redo the design. (Keep the other constraints on
s(n).) This will result in a filter where the final value of the step response is
1. Make plots of A(), h(n), and s(n). By adding the constraint A(0) = 1,
what is the increase in the value of the peak-error in the pass-band?
7.9 Constrained Minimax Filter Design: A Type II FIR filter has impulse
response
g = [1,
Note: To indicate with the lp that the first Ne constraints are equality
constraints, you can use the following syntax.
3,
3,
1]
x = lp(c,Q,b,[],[],[],Ne);
and
Q(Ne+1:end,:) * x <= b(Ne+1:end,:)
You should specify which of the constraints are equality constraints, and
which of the constraints are inequality constraints.
Also, could the Remez algorithm be used to solve this design problem?
56
Spectral Factorization
Find and sketch, h(n), a real-valued spectral factor of P (z). That is, find
h(n) such that
8.1 For the filter, H(z), whose zeros are shown in the following diagram, make
a sketch of the zero diagram of each spectral factorization of H(z).
ZPLANE(h)
1
8.3 The system P (z) has the zeros shown in the following diagram and unity
dc gain (meaning P (z = 1) = 1).
Imaginary Part
2
0.5
2
ZERO DIAGRAM
0.5
1.5
1
2
1
1
Real Part
IMAG
0.5
8.2 The sequence p(n) is symmetric (p(n) = p(n)). The zeros of P (z) are
shown in the following diagram. Also, P (z) has unity dc gain (meaning
P (z = 1) = 1).
0.5
ZERO DIAGRAM
1.5
2
1.5
0.5
REAL
0.5
1.5
0.5
Imaginary Part
Find and sketch the impulse response h(n) of a spectral factor of P (z).
0
0.5
You should be able to do this problem by hand (not using MATLAB, etc.).
(a) Sketch the impulse response p(n).
1.5
0.5
Real Part
0.5
57
(d) Perform spectral factorization on P (z). That means, find h(n) (realvalued) so that
X
p(n) =
h(k) h(k n)
(2)
(2)
0.5
2.0
5,
10,
10,
5,
The dc-gain of the filter is unity. (Note that the roots on the unit circle are
double zeros.)
1]
(b) Sketch the frequency response magnitude |H(ej )|. Accurately indicate the nulls of the frequency response.
(c) Find a spectral factorization of H(z). Accurately sketch the impulse
response of the spectral factor and its zero-diagram.
4,
6,
4,
1]
8.8 Spectral Factorization. The transfer function P (z) of a Type I FIR filter
has zeros in the z-plane as illustrated.
(2)
2,
(2)
3,
3,
2,
1]
(2)
2,
(2)
0.5
(2)
3,
2,
2.0
(2)
(2)
1]
8.7 The transfer function H(z) of a Type I FIR filter has zeros in the z-plane
as illustrated.
58
9.1 Use the lifting procedure to design a minimal-length minimum-phase lowpass FIR filter satisfying the specifications
||H(ej )| 1| p
0 p
(6)
for s
(7)
for
|H(ej )| s
|H(ej )| s1
||H(ej )| 1| p
j
|H(e )| s2
0 1
(8)
for 2 3
(9)
for 4
(10)
where
where
p = 0.02,
for
s = 0.01,
and p = 0.72,
s1 = s2 = 0.01,
s = 0.80.
and
Also design a linear-phase filter of the same length with the weight function
0 p
Kp
0
p < < s
W () =
Ks
s
1 = 0.30,
2 = 0.35,
3 = 0.60,
4 = 0.65.
Also design a linear-phase filter of the same length with the weight function
Ks1
0 < < 1
1 < < 2
0
Kp
2 < < 3
W () =
0
3 < < 4
Ks2
4 < <
Ks
Kp
W ()
p = 0.02
where
p s
Ks1 =
where
1
Kp =
,
p
1
,
s1
Kp =
1
,
p
Ks2 =
1
,
s2
(Note that the weight function will be different for the minimum-phase
filter. For the minimum-phase filter use the weight function derived in the
lecture notes.) Compare the minimum-phase and linear-phase filters to
each other in terms of the maximum pass-band and stop-band error. For
both filters, plot the impulse response, frequency response magnitude, and
zero diagram. Also plot the unwrapped phase and the group delay.
1
Ks =
,
s
(Note that the weight function will be different for the minimum-phase
filter. For the minimum-phase filter use the weight function derived in the
lecture notes.) Compare the minimum-phase and linear-phase filters to
each other in terms of the maximum pass-band and stop-band error. For
both filters, plot the impulse response, frequency response magnitude, and
zero diagram (be sure to use zplane correctly!). Also plot the unwrapped
phase and group delay; use the unwrap command.
and
s2 = 0.03.
s = 0.01,
and p = 0.12,
9.5 Matching. The diagrams on the following three pages show the impulse
responses, pole-zero diagrams, and frequency responses magnitudes of 4
discrete-time causal FIR filters. But the diagrams are out of order. Match
each set of diagrams by filling out the following table. NOTE: There are
only 2 distinct frequency response magnitudes among these 4 systems.
s = 0.20.
59
Impulse response
Pole-zero diagram
FREQUENCY RESPONSE 1
Frequency response
1
2
3
4
FREQUENCY RESPONSE 2
1.4
1.4
1.2
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0
1
0.2
0.5
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.6
0.4
0.5
0.3
0.4
0.2
0.3
0.1
0.2
0.1
0.1
0.2
0.1
0.3
0.2
0.4
0
10
15
20
25
30
10
IMPULSE RESPONSE 3
0.6
0.3
0.5
0.2
0.4
0.1
0.3
0.2
0.1
0.1
0.2
0.3
0.1
0.4
25
30
25
30
1.5
1.5
0.2
0
10
15
20
25
30
10
POLEZERO DIAGRAM 1
15
20
POLEZERO DIAGRAM 2
0.5
0.5
imag part
imag part
20
IMPULSE RESPONSE 4
0.4
30
30
0.5
0.5
1
1.5
0.5
0
real part
0.5
1.5
1.5
POLEZERO DIAGRAM 3
1
0.5
0.5
30
0.5
0
real part
0.5
POLEZERO DIAGRAM 4
imag part
imag part
15
30
0.5
0.5
1
1.5
0.5
0
real part
0.5
1.5
1.5
0.5
0
real part
0.5
60
0.5
0
1
0.5
0.5
9.6 Matching. The diagrams on the following three pages show the impulse
responses, pole-zero diagrams, and frequency responses magnitudes of 4
discrete-time causal FIR filters. But the diagrams are out of order. Match
each set of diagrams by filling out the following table.
Impulse response
Pole-zero diagram
FREQUENCY RESPONSE 1
Frequency response
1.2
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0
1
0.2
0.5
0.5
0
1
0.5
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.5
FREQUENCY RESPONSE 3
FREQUENCY RESPONSE 4
1.4
1.4
1.2
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.3
0.25
0.2
0.2
0.5
0
1
0.15
0.1
0
0.05
0.5
0.1
0
10
15
20
IMPULSE RESPONSE 3
10
15
20
15
20
IMPULSE RESPONSE 4
0.4
0.3
0.2
0.5
0.1
0
0.1
0.2
0.3
0.4
0.5
0
10
15
20
POLEZERO DIAGRAM 1
10
POLEZERO DIAGRAM 2
0.5
imag part
0.5
20
0.5
20
0.5
1
1.5
0.5
0
real part
0.5
1.5
1.5
POLEZERO DIAGRAM 3
0.5
0
real part
0.5
1.5
1.5
POLEZERO DIAGRAM 4
1
0.5
0.5
imag part
20
0.5
20
0.5
1
1.5
0.5
0
real part
0.5
1.5
0.2
0.5
0.05
imag part
1.4
0.2
1
2
3
4
imag part
FREQUENCY RESPONSE 2
1.4
1.5
0.5
0
real part
0.5
61
0.5
0
1
0.5
0.5
9.7 Matching. The diagrams on the following three pages show the impulse
responses, pole-zero diagrams, and frequency responses magnitudes of 8
discrete-time causal FIR filters. But the diagrams are out of order. Match
each set of diagrams by filling out the following table. NOTE: There are
only four distinct frequency response magnitudes among these 8 systems.
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.4
0.3
0.2
0.5
0.1
0
0.1
0.2
0.3
0.5
0.4
0
10
15
20
25
30
IMPULSE RESPONSE 3
10
15
20
25
30
25
30
25
30
25
30
IMPULSE RESPONSE 4
0.4
0.3
0.3
0.25
0.2
0.2
0.1
0.15
0.1
0.1
0.05
0.2
0.3
0.05
0.4
0.1
0
10
15
20
25
30
IMPULSE RESPONSE 5
10
15
20
IMPULSE RESPONSE 6
0.3
0.25
0.2
0.5
0.15
0.1
0.05
0
0.05
0.5
0.1
0
10
15
20
25
30
IMPULSE RESPONSE 7
1
0.5
0.5
Frequency response
1
2
3
4
5
6
7
8
62
20
0.5
0
Pole-zero
15
IMPULSE RESPONSE 8
0.5
Impulse response
10
10
15
20
25
30
10
15
20
POLEZERO DIAGRAM 1
FREQUENCY RESPONSE 2
1.4
1.4
1.2
1.2
0.5
0.5
imag part
imag part
FREQUENCY RESPONSE 1
POLEZERO DIAGRAM 2
30
30
0.5
0.5
0.8
0.8
0.6
0.6
0.4
0.4
0.2
1.5
0.5
0
real part
0.5
1.5
1.5
POLEZERO DIAGRAM 3
0.5
real part
0.5
0
1
1
0.5
imag part
imag part
30
0
0.5
30
0.5
0.5
1
1.5
0.5
0.5
real part
1.5
1.5
POLEZERO DIAGRAM 5
0.5
0
real part
0.5
1.4
1.2
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.5
imag part
30
0
0.5
30
0.5
1
2
1.5
0.5
0
real part
0.5
1.5
0.5
POLEZERO DIAGRAM 7
0.5
1
real part
1.5
POLEZERO DIAGRAM 8
0.5
imag part
0.5
30
0.5
30
0.5
1
1.5
0.5
0
real part
0.5
1.5
1.5
0.5
0
real part
0.5
1.5
63
0.5
0.2
0.5
POLEZERO DIAGRAM 6
0.5
0.5
FREQUENCY RESPONSE 4
1.4
0
1
1.5
0
1
0.2
imag part
FREQUENCY RESPONSE 3
POLEZERO DIAGRAM 4
0.5
0.5
imag part
0.2
0.5
0
1
0.5
0.5
9.8 Matching. The following diagrams show the impulse responses, pole-zero
diagrams, frequency responses magnitudes, and group delay responses of
four discrete-time causal FIR filters. But the diagrams are out of order.
Match each set of diagrams by filling out the following table. Each impulse
response is 31 samples long.
POLEZERO DIAGRAM 1
POLEZERO DIAGRAM 2
0.5
imag part
imag part
0.5
30
30
0.5
0.5
1
1
0.5
0.5
real part
1.5
1.5
POLEZERO DIAGRAM 3
imag part
imag part
30
0.5
1
1.5
1
2
3
4
0.5
0
0.5
real part
1.5
1.5
FREQUENCY RESPONSE 1
IMPULSE RESPONSE 1
30
Group delay
1.5
0.5
0.5
Freq. resp.
POLEZERO DIAGRAM 4
0.5
P/Z
0
0.5
real part
Impulse resp.
0.5
0.5
0
real part
0.5
1.5
FREQUENCY RESPONSE 2
1.4
1.4
1.2
1.2
IMPULSE RESPONSE 2
0.5
0.5
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.5
0
1
0.5
0.2
0.5
10
15
20
25
30
10
15
20
25
0.5
0
1
30
0.5
FREQUENCY RESPONSE 3
IMPULSE RESPONSE 3
0.5
1
0
FREQUENCY RESPONSE 4
1.4
1.4
1.2
1.2
IMPULSE RESPONSE 4
0.5
0.5
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
1
0.5
5
10
15
20
25
30
10
15
20
25
0.5
0.5
0
30
64
0.5
0
1
0.5
0.5
GROUP DELAY 1
(a) Design a notch filter with a single zero at the notch frequency (accordingly, the linear-phase filter must have a double zero at the notch
frequency).
GROUP DELAY 2
30
30
25
25
20
20
15
15
10
10
(b) Design a notch filter with a double zero at the notch frequency (accordingly, the linear-phase filter must have a zero at the notch frequency
of multiplicity 4).
0
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
GROUP DELAY 3
GROUP DELAY 4
30
30
25
25
20
20
15
15
10
10
0.8
Compare your notch-filter with the linear-phase ones given in the earlier
notes. For both filters, plot the impulse response, frequency response magnitude, and zero diagram. Also plot the unwrapped phase and the group
delay.
10
0
0
0.2
0.4
0.6
/
0.8
0.2
0.4
0.6
0.8
10.1 Consider the design of an analog filter to meet the following requirements.
9.9 Consider the design of a minimum-phase FIR low-pass filter with a specified
null in the stop-band, as we have seen in previous notes on linear-phase FIR
filters.
(a) What problem arises when you try to use the lifting procedure to
design this type of filter?
for || 4
for || 4.6
(b) (Optional.) An alternative to the lifting procedure is to use an approach based on a linear program formulation. Develop an approach
to design minimum-phase FIR filters based again on the spectral factorization of a Type I FIR filter. Design a minimum-phase FIR filter
with the same specification and length as the example described in the
notes on Chebyshev linear-phase FIR design and compare your result
to that filter.
9.10 In previous notes we designed linear-phase FIR notch filters with a specified notch frequency. For this problem, you will design minimum-phase
FIR notch filters. This can be done by performing spectral factorization of
a suitably designed Type I linear-phase FIR notch filter. It will not be necessary to lift the frequency response of the notch filter because it already
has A() 0. Recall the examples in previous notes. Specifically, design
a minimum-phase notch filter of length 51 with notch frequency n = 0.6
and band-edges 0 = 0.55, 1 = 0.65 where
0 0
1
0
0 < < 1
D() = 1,
and
W () =
1
1 .
Ga (s) = Ha (C s)
This can be done by rescaling the poles and zeros by a suitable scaling
constant,
z = z/C;
p = p/C;
where z and p contain the zeros and poles of the transfer function Ha (s).
You will need to determine the appropriate value of C based on the desired
band-edges and the filter type.
10.4 (From Mitra M7.5.) Design a digital Chebyshev-I lowpass filter operating
at a sampling rate of 80 kHz with a passband edge frequency at 4 kHz, a
passband ripple of 0.5 dB, and a minimum stopband attenuation of 45 dB
at 20 kHz using the bilinear transformation method. Determine the order
of the analog prototype using the command cheb1ord and then design the
analog prototype using cheb1ap. Transform the analog filter into a digital
one using the bilinear command. Plot the frequency response, pole-zero
diagram and impulse response of the digital filter. Show the steps of the
design procedure.
10.5 BLT of the integrator. Apply the bilinear transformation to the ideal
continuous-time integrator
Ha (s) =
1
s
(b) Suppose you have an unstable causal analog filter. Is it possible that
the bilinear transformation will produce a stable causal digital IIR
filter?
s=C
z1
z+1
1
.
Ha (s) =
(s + 1)(s + 2)
(d) Sketch the pole-zero plot of H(z) and comment on its stability.
10.6 Impulse invariance method. The bilinear transformation (BLT), for
converting an analog filter into a digital filter, uses a transformation of
variables,
10.3 Consider the design of a recursive digital filter whose frequency response
satisfies:
1
0 0.3
f
H ()
0
0.5
s=C
z1
.
z+1
Let us use C = 2/T where T is the sampling interval for the digital filter.
That means, the passband edge is p = 0.3 and the stopband edge is
s = 0.5 . Suppose that the recursive digital filter is to be obtained from
an analog filter Ha (s) using the bilinear transformation with
s=
Another way to convert an analog filter into a digital filter is the impulse
invariant method. This problem compares the two methods applied to the
following analog filter,
z1
.
z+1
Ha (s) =
Then what should be the passband and stopband edges of the analog filter
Ha (s) so that after the bilinear transformation the resulting digital filter
will have the desired band edges?
1
s+3
(a) Use the BLT, to transform the analog filter into a digital filter. Find
the transfer function and sketch its pole/zero diagram.
10
(b) The impulse invariance method samples the analog impulse response,
to obtain the impulse resonse of a digital filter,
IMAGINARY
h(n) = ha (nT ),
Applying this method to the analog filter above, find the transfer function of the resulting digital filter and sketch its pole/zero diagram.
10
(c) For the two digital filters obtained in parts (a) and (b), what is the
difference between the frequency responses?
(d) For a higher order analog system Ha (s), the impulse invariance method
can be applied to first order terms in a partial fraction expansion of
Ha (s). If the impulse invariance method is applied to a stable analog
filter, is it guaranteed that the digital filter will be stable?
0
REAL
(a) Express in hertz the frequencies that are annihilated (blocked) by the
system.
(b) Find the transfer function H(s) of the system. Simplify.
(c) Sketch the frequency response magnitude |H(j)|. Explain how it can
be sketched directly from the pole-zero diagram rather than direct
numerical calculation. What kind of filter is the system?
10.7 BLT and FIR filters. The bilinear transformation (BLT), for converting
an analog filter into a digital filter, uses a transformation of variables,
s=
Part 2. Design a second-order discrete-time system, operating at a sampling frequency of 30 samples/second, that approximates the continuoustime system in Part 1. Your discrete-time system should annihilate (block)
the same frequencies as the continuous-time system. Clearly describe your
approach and your steps.
z1
.
z+1
Usually, the BLT produces a recursive digital filter. This problem investigates this issue.
(a) What are the poles and zeros of H(z) of your discrete-time system?
Sketch the pole-zero diagram.
5
,
6
1
]
6
(c) Give the difference equation that implements your discrete-time system.
(b) Use the BLT in reverse, to transform the FIR digital filter into an
analog filter. Find the transfer function Ha (s) of the analog filter.
10.9 The following figure shows the pole-zero plots of 6 different digital IIR
filters. By inspecting the pole-zero plots,
10.8 Part 1. A continuous-time LTI system with transfer function H(s) has the
No
pole-zero diagram shown. The system has a dc gain of unity.
Matlab
needed.
67
FILTER B
1
0.5
0.5
Imaginary Part
Imaginary Part
FILTER A
1
0.5
1
0.5
1
0.5
0
0.5
Real Part
0.5
0.5
0.5
0
0.5
1
0
0.5
1
0.5
0
0.5
Real Part
0.5
FILTER E
0
0.5
Real Part
FILTER F
0.5
0.5
Imaginary Part
Imaginary Part
FILTER D
Imaginary Part
Imaginary Part
FILTER C
0
0.5
Real Part
0
0.5
1
0.5
1
0.5
0
0.5
Real Part
0.5
0
0.5
Real Part
68
2
y(n 1)
3
is used in forwards-backwards filtering, then what is the transfer function of the total system being implemented? Sketch its pole-zero diagram.
(b) Sketch the frequency response magnitude for the system implemented
using forwards-backwards filtering.
(c) Find the impulse response of the system implemented using forwardsbackwards filtering.
11
Multirate Systems
x(n) = {. . . , 0, 0, 1, 2, 3, 2, 1, 0, 0, . . . } ,
h(n) = {1, 2}
is used with the MATLAB filtfilt command, then what is the transfer function of the system being implemented? Sketch its pole-zero
diagram.
(c) Sketch the frequency response magnitude for the system in (b).
x(n) = {. . . , 0, 0, 1, 2, 1, 0, 1, 0, 0, . . . }
10.13 If the causal filter
1
1 c z 1
11.5 If h(n) is the impulse response of a linear-phase FIR filter, are the two
polyphase components of h(n) linear-phase as well? Consider separately
the case when the length N is even and odd.
The path from the input signal to each of the four output signals can be
rewritten using the Noble identities as
11.6 (From Mitra 10.20) The running-sum filter, also called the boxcar filter,
H(z) =
N
1
X
z n
for 1 k 4.
n=0
H(z) = (1 + z 1 ) (1 + z 2 ) (1 + z 4 ) (1 + z 2
(b) Given the frequency response of h and g shown below, sketch the
frequency responses of the four transfer functions Fk (z).
11.7 In the following discrete-time multirate system the filter h is a lowpass
filter, and the filter g is a highpass filter. The system produces four subband
signals, we can call them s1 (n), . . . , s4 (n).
(c) Classify each of the four transfer functions Fk (z) as lowpass, bandpass,
bandstop, or highpass.
70
(b)
H FREQUENCY RESPONSE
2 2
1
0.8
(c)
0.6
0.4
5 10 2
0.2
0
0.6
0.4
0.4
0.6
G FREQUENCY RESPONSE
1
0.8
0.6
0.4
0.2
0.6
0.4
0.4
0.6
11.8 Can you change the order of an up-sampler and a down-sampler with out
change the total system? In other words, are the following two systems
equivalent?
11.11 Rate changing. A discrete-time signal x(n) having a rate of 100 samples
per second must be converted into a new discrete-time signal y(n) having
a rate of 90 samples per second.
System A:
M L
(a) Sketch a multirate system that performs the necessary sampling rate
conversion. Sketch the frequency responses of the required filter(s).
System B:
(b) Because ideal brick-wall filters can not be realized, how would you
proceed with the design of filters(s) for this multirate system?
L M
Consider the following cases:
(a) M = 2, L = 2.
(b) M = 2, L = 3.
(c) M = 2, L = 4.
Try an example in each case with a simple test input signal.
Sketch a multirate system that performs the necessary sampling rate conversion. Sketch the frequency responses of the required filter(s).
71
||
,
For (b) and (c), the Nyquist property can be included in the design by
adding constraints on h(n).
for ||
Suppose we generate the sequences y(n) and s(n) from x(n) with the following system
11.16 Let h(n) be a low-pass half-band filter with real coefficients. Then show
that
s(n)
H f () + H f ( ) = 2.
H f () =
1, || < /2
0, /2 || <
(That means the ripple in the pass-band of H f () is the same as the ripple
in the stop-band of H f ()). You can assume the Nyquist filter is centered
at no = 0.
Here, the DC gain will be about 2, rather than 1 as is usual for a low-pass
filter. What must the cut-off frequency o be? (H f (o ) = 1) Does this
explain the terminology half-band ?
H(z) = Q(z) (1 + z 1 + z 2 )3 .
(a) Find the impulse response h(n) of the filter.
s(n)
(c) What particular properties does the filter have, when it is used for
interpolation?
1, || < /2
0, /2 || <
f
11.18 Half-band filter design. Design a Type I linear-phase FIR digital halfband filter of minimal length having a transfer function H(z) of the form
H(z) = Q(z) (1 + z 1 )2 (1 + z 1 + z 2 ).
11.15 Linear-phase FIR Nyquist-L filter design. Using Matlab, design a
Type I linear-phase Nyquist-3 FIR filter of length 29. A Nyquist-3 filter
h(n) is one for which every third coefficient is zero, except for one such
value. (See the notes on multirate systems.) The cut-off frequency o
should be /3. Use the following design methods.
20
, s = o +
11.19 Half-band filter design. Design a Type I linear-phase FIR digital halfband filter of minimal length having a transfer function H(z) of the form
20 .
Use Kp = Ks = 1 and p = o 20
, s = o + 20
.
H(z) = Q(z) (1 + 2z 1 + z 2 ) (1 + 3z 1 + z 2 ).
72
(d) What particular properties does the filter have, when it is used for
interpolation?
for some integers K and L? If so, show a derivation and identify K and L.
You can assume the Nyquist filter is centered at n = 0, i.e. h(4n) = (n).
Do this by downsampling the output y(n) of the 2X-interpolator and showing that it is the same as x(n) when a halfband filter is used.
0
0.5
1
(a) What condition should the impulse response of the interpolation filter
satisfy?
10
15
20
25
30
35
(a) h(n) = [1
1]/2
(b) h(n) = [1
1]/4
(c) h(n) = [1
1]/8
That Matlab code to genearate the test signal shown above is available on
the course webpage.
73
x(n)
H(z)
H(z)
y(n)
x(n)
G(z)
y(n)
(d) Given an expression for the impulse response h(n) of the total system in terms of f (n) and g(n). Explain why this system is called an
Interpolated FIR filter.
x(n)
AVE/
DIFF
AVE/
DIFF
d(n)
AVE/
DIFF
d2 (n)
c3 (n)
d3 (n)
H(z) = H1 (z M ) H2 (z).
Sometimes an IFIR filter can meet given specifications with a lower implementation complexity than a generic FIR filter.
AVE/
DIFF
c(n)
d(n)
H(z) = H1 (z 4 ) H2 (z)
(a) Draw a block diagram for obtaining c(n) from x(n). The block diagram
may combine an up-sampler, a down-sampler, and an LTI filter; but
not more than one of each. What is the impulse response of the LTI
filter? Sketch the zeros of the filter.
Similarly, draw a block diagram for obtaining d(n) from x(n).
and that H1 (z) is the ideal lowpass filter with cut-off frequency 0.5 . That
means, H1 (z) is given, and H2 (z) is to be designed so that the total filter
H(z) is a lowpass filter with cut-off frequency 0.125 = /8.
(b) Draw a block diagram expressing the system between x(n) and d2 (n)
in the 3-level system above. The block diagram may combine upsamplers, down-samplers, and LTI filters; but it should not have more
than one of each. (Hint: use part a and the Noble Identities.) Sketch
the impulse response and the zeros of the LTI system. Based on the
zero diagram, roughly sketch the frequency response of the LTI system.
Show your work.
11.27 The following system consists of a rate changer in sequence with itself.
74
12
Quantization
(a) Choose x2 = 0.01 and use 10000 samples. Measure and plot the SNR
in dB as a function of the number of bits b for 1 b 12. Check to
see if the plot you obtain agrees with the expression obtained in the
lecture notes,
!
2
R
f
s
SNR = 10 log10 x2 + 20 b log10 2 10 log10
12
12.1 Scaling. Suppose the signal x(n) is bounded by 0.5, i.e. |x(n)| < 0.5 for
all n. The impulse response of an FIR digital filter is
h(n) = {1.5, 2.0, 0.5, 0.3, 0, 0, . . . }
To ensure that the output signal y(n) does not exceed 1 in absolute value,
how should x(n) be scaled before filtering? Show your work.
by plotting both curves on the same graph. (Note that fxquant uses
Rf s = 2.) Does the SNR improve by 6 dB per bit?
Repeat with x2 = 0.1. Explain any deviation you find from the SNR
formula. What assumptions are violated?
12.2 Scaling. The DTFT X f () of signal x(n) is bounded by 2, i.e. |X f ()| < 2
for all . It is filtered with an LTI system with frequency response
1, || 0.3
H f () =
0, 0.3 < ||
To ensure that the output signal y(n) does not exceed 1 in absolute value,
how should x(n) be scaled before filtering? Show your work.
12.3 You have developed a system using an 8 bit uniform quantizer. For this
system, the SNR (ratio of signal power to quantization noise power) turns
out to be 10 dB, although the target SNR is 20 dB. How many additional
bits do you expect is needed in the uniform quantizer in order to meet that
goal?
12.5 Quantization. If a white Gaussian random discrete-time signal with variance 2 is quantized using a uniform quantizer with N bits, then under
standard assumptions, the quantization error is uniformly distributed.
The following figure shows quantization error histograms for three different
combinations of N and . The quantizer provides an output in the range
[1, 1). Identity which histogram results from each of the three combinations by completing the table.
12.4 SNR of a quantizer. In this problem we assume the original signal x(n) is
Gaussian (normally distributed). We will measure the SNR
SNR = 10 log10
x2
e2
Explain in each case how the shape of the histogram is influenced by N and
.
under different conditions. We will use the Matlab command randn to generate random numbers according to the Gaussian bell shape distribution.
In each case, is the true SNR likely to be higher than, lower than, or in
agreement, with the SNR predicted using the assumption that the quantization error is uniformly distributed?
Number of bits
2
4
4
0.1
0.1
1.0
Histogram
13.2 Power Spectral Densities. The figure below shows the power spectral
density (PSD) of four stationary distrectre-time random processes. The
following figures also shows a realization (signal) generated using each of
the four PSDs, but they are out of order. Match each signal to its most
likely PSD by completing the table.
HISTOGRAM 1
x 10
4
3
2
1
0
1
0.8
0.6
0.4
0.2
0
err/Q
0.2
0.4
0.6
0.8
HISTOGRAM 2
x 10
6
PSD
SIGNAL
0
1
0.8
0.6
0.4
0.2
0
err/Q
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
1
2
3
4
HISTOGRAM 3
x 10
10
0
1
0.8
0.6
0.4
0.2
0
err/Q
13
Spectral Estimation
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
h(n)
/4
/2
3/4
y(n)
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
76
/4
/2
3/4
/2
3/4
/4
/4
/2
3/4
13.3 Power Spectral Densities. The figure below shows the power spectral
density (PSD) of four different stationary discrete-time random processes.
The following figure also shows a realization (signal) generated using each
of the four PSDs, but they are out of order. Match each signal to its most
likely PSD by completing the table. Explain your answers.
SIGNAL 1
1
0
1
2
20
40
60
80
100
120
140
160
180
200
SIGNAL 2
1
0
PSD
1
2
20
40
60
80
100
120
140
160
180
200
20
40
60
80
100
120
140
160
180
200
SIGNAL
1
2
3
4
SIGNAL 3
1
0
1
2
SIGNAL 4
1.2
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
20
40
60
80
100
120
140
160
180
200
/4
/2
3/4
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
77
/4
/2
3/4
/2
3/4
/4
/4
/2
3/4
13.4 Power spectral densities. The figures below show the power spectral
density (PSD) of four stationary distrectre-time random processes. The
following figures also show a realization (signal) generated using each of
the four PSDs, but they are out of order. Match each signal to its most
likely PSD by completing the table.
SIGNAL 1
1
0
1
2
20
40
60
80
100
120
140
160
180
200
SIGNAL 2
1
0
PSD
1
2
20
40
60
80
100
120
140
160
180
200
20
40
60
80
100
120
140
160
180
200
Signal
1
2
3
4
SIGNAL 3
1
0
1
2
SIGNAL 4
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
20
40
60
80
100
120
140
160
180
200
0
0
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
0.5
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
78
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
13.5 Power spectral densities. The figures below show the power spectral
density (PSD) of four stationary distrectre-time random processes. The
following figures also show a realization (signal) generated using each of
the four PSDs, but they are out of order. Match each signal to its most
likely PSD by completing the table.
SIGNAL 1
1
0
1
2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
SIGNAL 2
1
0
PSD
1
2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
Signal
1
2
3
4
SIGNAL 3
1
0
1
2
1.2
SIGNAL 4
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
0
0
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
0.5
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
79
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
13.6 Power spectral densities. The figures below show the power spectral
density (PSD) of four stationary discrete-time random processes and a realization (signal) generated using each of the four PSDs, but they are out
of order. Match each signal to its most likely PSD.
SIGNAL 1
1
0
1
2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
SIGNAL 2
1
0
1
2
SIGNAL 3
1
0
1
2
1.2
SIGNAL 4
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
0
0
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
0.5
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
80
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
0.1
0.2
0.3
0.4
NORMALIZED FREQUENCY
0.5
13.7 Power spectral densities. The figures below show the power spectral
density (PSD) of four stationary discrete-time random processes and a realization (signal) generated using each of the four PSDs, but they are out
of order. Match each signal to its most likely PSD by completing the table.
SIGNAL 1
1
0
1
SIGNAL 2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
PSD
1
2
3
4
0
1
2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
Signal
SIGNAL 3
1
0
1
2
1.2
SIGNAL 4
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
20
40
60
80
100
120
TIME (SAMPLE)
140
160
180
200
0
0
/4
/2
3/4
/2
3/4
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
81
/4
/4
/2
3/4
/4
/2
3/4
SIGNAL 1
1
0
1
2
0
20
40
60
80
100
120
140
160
180
200
SIGNAL 2
13.9 Power spectral densities. The discrete-time signal x(n) is a white stationary random sequence with zero-mean and unit variance. The signal
x(n) is filtered with a discrete-time LTI system to produce output signal
y(n). What is the variance of y(n) in each of the following cases?
1
0
1
2
0
20
40
60
80
100
120
140
160
180
200
H() = 1,
|| .
SIGNAL 3
1, || 0.4
H() =
0, 0.4 < || .
1
0
1
2
0
20
40
60
80
100
120
140
160
180
200
|| 0.4
1,
H() =
SIGNAL 4
1
0
1
2
0
20
40
60
80
100
120
140
160
180
200
H(z)
x(n)
The autocorrelation function rx (k) of the generated random signal x(n) has
the values:
rx (0) = 3,
82
rx (1) = 2,
rx (2) = 1.
14
Speech Filtering
14.1 Speech Recording. In this exercise, you are to record yourself saying
yes and no. You can create your recordings using a computer or other
digital recording device (some mp3 players have a recording function). In
case it is not available in the operating system, there are many utilities for
both MS and Mac systems. For example, Audacity is a free recording and
audio editing software system for both Mac and MS.
Audacity: a free cross-platform audio recorder and editor.
http://audacity.sourceforge.net/
Using the recording utility, you can create a wav file; most utilities can
export to a wav file. You can then read the wav file into MATLAB using
the wavread command. Your recordings should be at 8000 samples per
second. If they are not at that sampling rate, then you can change them
after you load them into MATLAB using the resample command.
Using MATLAB, you can form one vector of 4000 samples (half second)
for the yes. You should form a second vector also of 4000 samples of the
no. Note that the recording you originally make does not need to be a
half second in duration; you can trim the signal down to 4000 samples after
you read the signal into MATLAB.
(b) Display the PSD of the signal (computed using the filter you use).
(c) Estimate the PSD from the signal using the periodogram.
(d) Estimate the PSD from the signal using Welchs periodogram.
You can save your 4000-point signal as a simple data file using the save
command. For example:
(e) Estimate the PSD from the signal using AR spectral estimation.
To verify that your data file works correctly, you can use the following
commands to listen to the audio and plot the signal.
>> clear
>> load no_data.txt
>> whos
Name
Size
no_data
4000x1
>>
>>
>>
>>
13.12 The signal in the file cat_8KHz.wav is recorded using a laptop microphone.
There is noise present in the signal. Using the first 10,000 samples of the
signal (during which there is no speech) estimate the PSD of the noise.
Use the periodogram, Welchs periodogram, and AR spectral estimation.
Display the estimated PSD on linear and on logarithmic scales. Comment
on your observations.
fs = 8000;
plot( (1:4000)/fs, no_data );
xlabel( Time (sec) )
sound( no_data, fs )
Bytes
32000
Class
double
Attributes
% sampling rate
Everyones signal will be different. To give you a sense of what you might
get, my data is shown in the figure.
83
To turn in: For each of your yes and no signals, you should turn in
plots of your speech signal, the selected 50 millisecond segment, and the
spectrum on linear and dB scales. For example, see the following figures
for my speech signal.
0.05
0.1
0.15
0.2
14.3 Notch filter. In the previous speech exercises, you have seen that the
spectrum of a voiced segment of your speech signal contains peaks which are
approximately equally spaced. In this exercise you will design a notch filter
and apply it to your speech signal so as to eliminate one of the prominent
peaks in the spectrum (without affecting the adjacent peaks).
0.25
0.3
0.1
0.2
0.3
0.4
0.5
Time (sec)
no at 8000 samples/sec
0.15
0.1
0.05
0
0.05
0.1
0.15
0.2
H(z) =
0
0.1
0.2
0.3
0.4
0.5
b0 + b1 z 1 + b2
2
1 + a1 z 1 + a2 z 2
Time (sec)
14.2 Frequency Analysis of Speech Segments. In the previous speech exercise, you recorded yourself saying yes and no. In this exercise you
will perform frequency analysis on your two recordings. Specifically, you
will compute and display the spectrum of one segment of each of your two
signals.
(b)
(c)
(d)
(e)
84
of a second order digital notch filter, that has zeros at the selected
frequencies and corresponding poles. The zeros should be at ej2fn
and the poles should be at r ej2fn where r is slightly less than 1.
You can try different values of r.
Plot the pole-zero diagram of your filter (use the Matlab command
zplane) in Matlab. Verify that the poles and zero match where they
were designed to be.
Plot the frequency response magnitude of your filter |H f (f )| versus
physical frequency (the frequency axis should go from zero to half the
sampling frequency). You can use the Matlab command freqz to
compute the frequency response. Verify that the frequency response
has a null at the intended frequency.
Plot the impulse response h(n) of your filter. You can create an impulse signal and then use the Matlab command filter to apply your
filter to the impulse signal.
Apply your filter to your speech signal (use Matlab command filter).
Extract a short segment of your speech signal before and after filtering. Plot both the original speech waveform x(n) and filtered speech
no (8000 samples/sec)
0.5
0.5
0.5
0.5
0.05
0.1
0.15
0.2
0.25
0.3
Time (sec)
50 msec frame
0.35
0.4
0.45
0.5
0.5
0.5
0.5
0.5
1
0.05
0.055
0.06
0.065
0.07
0.075
Time (sec)
0.08
0.085
0.09
0.095
1
0.1
0.1
0.05
0.1
0.15
0.2
0.105
0.11
0.115
0.12
0.25
0.3
Time (sec)
50 msec frame
0.125
Time (sec)
0.13
0.35
0.4
0.45
0.5
0.135
0.14
0.145
0.15
60
50
40
40
30
20
20
10
0
4000
3000
2000
1000
0
Frequency (Hz)
1000
2000
3000
0
4000
4000
3000
2000
20
20
3000
2000
1000
0
Frequency (Hz)
1000
0
Frequency (Hz)
1000
2000
3000
4000
2000
3000
4000
40
20
4000
1000
2000
3000
20
4000
4000
85
3000
2000
1000
0
Frequency (Hz)
1000
waveform y(n) (you might try to plot them on the same axis). Also
plot the difference between these two waveforms d(n) = y(n) x(n).
What do you expect the signal d(n) to look like?
For the short segment you extract from the original and filtered speech
waveforms, compute and plot the spectrum. Were you able to eliminate the spectral peak that you intended to?
To turn in: plots of your filter (frequency response, impulse response, polezero diagram), plots of your speech signal (short segment) before and after
filtering, plots of the spectrum (of the short segment) before and after
filtering, discussion.
86
15
More Exercises
15.1 Several linear-phase FIR filters are shown in the following figure. Match
each impulse response with its zero diagram and frequency response. Also
classify each as Type I, II, III, IV.
In this exercise, use the same second order recursive filter and the same
speech waveform, but this time apply the filter to the data using forwardbackward filtering with the MATLAB command filtfilt. This implements the filter as a non-causal zero-phase filter. As before, to examine the
phase-distortion, compare the filtered signal with the original signal (subtract one from the other). Also, plot the impulse response of the filter. Because the filter is non-causal, to compute the impulse response, you should
use as the input to the filter an impulse signal of the form [zeros(1,L) 1
zeros(1,L)] or such (the 1 should not be at the beginning or end of the
input signal because in that case you will not see both sides of the two-sided
impulse response).
1 A
|Hf()|
ZERO DIAGRAM
IMPULSE RESPONSE
H1
1
0
0.5
0.5
0.5
0.5
0.5
0.5
/
H2
1 B
1
0
1 C
H3
0
2
1
0
1 D
H4
0
2
1
0
1 E
H5
0
2
1
0
2
0
1
0
87
H6
1 F
ZERO DIAGRAM 1
ZERO DIAGRAM 2
1
14
15
1
1
Find and sketch the DTFT X (). Hint: use multirate concept.
4
ZERO DIAGRAM 3
ZERO DIAGRAM 4
1
15.3 Match each impulse response with its zero diagram and frequency response.
1
0.5
IMPULSE RESPONSE 1
IMPULSE RESPONSE 2
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.2
14
15
0
0.5
1
1
1
ZERO DIAGRAM 5
0.4
5
10
15
IMPULSE RESPONSE 3
10
0.6
0.4
0.4
0.2
0.2
0.2
0.2
ZERO DIAGRAM 6
15
0.5
IMPULSE RESPONSE 4
0.6
0.4
0
14
14
0
0.5
1
1
3
ZERO DIAGRAM 7
0.4
ZERO DIAGRAM 8
0.4
0
10
15
IMPULSE RESPONSE 5
10
15
IMPULSE RESPONSE 6
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.2
0.4
5
10
15
IMPULSE RESPONSE 7
10
15
IMPULSE RESPONSE 8
0.6
0.4
0.4
0.2
0.2
0.2
0.2
0.4
0.4
0
10
15
10
15
88
14
0.6
14
0.4
0
FREQUENCY RESPONSE 1
FREQUENCY RESPONSE 2
1.5
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 3
0.2
0.4
0.6
0.8
D() = ej ,
FREQUENCY RESPONSE 4
1.5
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 5
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 6
1.5
|| <
0.8
0.6
15.6 Fractional delay system. The ideal fractional delay system has the frequency response
0.4
0.5
0.2
0
D() = ejd ,
0
0
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 7
0.2
0.4
0.6
0.8
FREQUENCY RESPONSE 8
1.5
1.5
0.5
0.5
Consider the design of a real FIR filter to implement an approximate fractional delay system with delay d = 9.3 samples and with impulse response
length N = 21. For a real FIR filter, the impulse response must be realvalued. Consider the design of such a filter by interpolation using the DFT.
In this case, the N -point impulse response h(n) is given by the inverse DFT
of an N -point sequence, the later sequence representing N equally spaced
values of the frequency response between = 0 and 2.
0
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
|| <
K
X
N
j
d
H
h
=
=
=
=
=
21;
sqrt(-1);
9.3;
exp(-j*d*2*pi*(0:N-1)/N);
ifft(H);
k=0
The code almost works. However, the impulse response produced by this
code is not real-valued. The impulse response (complex-valued) is shown
in the following figure:
The following is an EOG signal measured from the head. The EOG signal measures an electrical signal related to eye movement. The positive
(negative) spikes correspond to a sudden movement of the eye to the left
(right).
0.5
0
0.5
EOG SIGNAL
6
1
0
10
n
12
14
16
18
20
0.5
4
6
1000
2000
3000
4000
5000
TIME (SAMPLE)
6000
7000
8000
0.5
1
0
10
n
12
14
16
18
For this exercise, design a lowpass filter of your own choosing to remove as
much of the noise as possible, while maintaining the spikes in the signal.
This signal is available on the course webpage (data.m or data.mat they
contain the same data).
20
(a) Explain why the resulting impulse response is not real valued.
Provide plots of signals and the filters you used and comment on your
observations.
0 n 5.
If the input
x(n) = C cos(o n)
produces the output
y(n) = D cos(o (n no )).
What is no ?
90
15.11 Total variation filtering. Given a noisy signal y = [y(1), . . . , y(N )], the
method of total variation (TV) filtering attempts to remove the noise by
finding the signal x = [x(1), . . . , x(N )] that minimizes the cost function
F (x):
In some bio-signal processing tasks it is necessary to automatically determine when an event occurs. Many algorithms for this detection problem
begin by differentiating the signal (and proceed by thresholding). Because
differentiation is a linear time-invariant process, it can be viewed as an LTI
system (as convolution). See the text book for a description of the design
of filters for digital differentiation. The Matlab remez program can be used
for example. Note that a linear-phase FIR differentiator must by either a
Type III or Type IV (why?).
F (x) =
N
X
|y(n) x(n)|2 +
n=1
N
X
n=2
For the EOG signal in the previous problem, perform digital differentiation
using an appropriate linear-phase FIR filter. How well does this work? (It
should work poorly.)
The figure below shows a noisy signal and the result of TV filtering using
two different values of : = 0.3 and = 3. Which signal was produced
by each of the two values of ? Explain the reason.
Noisy signal
3
2
Provide plots of signals and the filters you used and comment on your
observations.
0
0
20
40
60
80
100
120
140
160
180
200
160
180
200
160
180
200
15.10 You are asked to design a digital filter to perform integration (as opposed
to differentiation).
20
40
60
80
100
120
140
1
0
0
20
40
60
80
100
120
140
15.12 STFT. Assume the window function is used in both the forward and inverse
STFT. Which of the following 8-point window functions (if any) satisfy
the perfect reconstruction condition for the STFT implemented with 50%
overlapping? Explain.
91
(a) w = [1
1]/sqrt(2)
(b) w = [0.6
0.6
0.8
0.8
0.8
0.8
0.6
0.6]
(c) w = [0.2
0.4
0.6
0.8
0.8
0.6
0.4
0.2]
15.13 Short-Time Fourier Transform. The inverse of the STFT can be computed using the same window that is used in the forward STFT provided
the window satisfies appropriate properties.
16
(a) When 50% overlapping is used, what property should the window satisfy?
Old Exercises
(b) Given an example of such a window (other than the one used in the
class handout).
X(f ) = 0,
(c) When the overlapping is by 2/3 (66.6%) what property should the
window satisfy?
for |f | > 50 Hz
You need to filter the signal with a lowpass filter with a cut-off frequency
at 30 Hz. You will do this by sampling the analog signal and applying a
digital lowpass filter. You will sample the signal with a sampling rate of 150
Hz. In terms of normalized frequency, what should be the cut-off frequency
o of the digital filter? (0 < o < ).
92