Lab    #8:    Frequency    Response

ECESe 352    

Lab    #8:    Frequency    Response:    Bandpass    &    Nulling    Filters
(Lab Report due at beginning of next lab)
Formal    Lab    Report:        You    must    write    a    formal    lab    report    that    describes    your    
approach. You    should    read    the    Pre-Lab    section    of    the    lab    and    do    all    the    
exercises    in    the    Pre-Lab    section    before    your    assigned    lab    time.    
Important:        When    it    instructs    you    to    get    an    updated    matlab    file    (like    
specgram.m),        download    
http://dspfirst.gatech.edu/matlab/toolbox/
Introduction
The    goal    of    this    lab    is    to    study    the    response    of    FIR    filters    to    inputs    such    as    complex    
exponentials    and    sinusoids.    In    the    experiments    of    this    lab,    you    will    use    firfilt(),    or    
conv(),    to    implement    filters    and    freqz()    to    obtain    the    filter’s    frequency    response.1    As    
a    result,    you    should    learn    how    to    characterize    a    filter    by    knowing    how    it    reacts    to    
different    frequency    components    in    the    input.    
This    lab    also    introduces    two    practical    filters:    bandpass    filters    and    nulling    filters.    
Bandpass    filters    can    be    used    to    detect    and    extract    information    from    sinusoidal    
signals,    such    as    tones    in    a    touch-tone    telephone dialer.    Nulling    filters    can    be    used    to    
remove    sinusoidal    interference,    such    as    jamming    signals    in    a    radar.    
1.1 Frequency Response of FIR Filters
The    output    or    response    of    a    filter    for    a    complex    sinusoid    input,    e^{jwn},    depends    on    
the    frequency,    w.    Often    a    filter    is    described    solely    by    how    it    affects    different    input    
frequencies—this    is    called    the    frequency    response.
For    example,    the    frequency    response    of    the    two-point    averaging    filter:
y[n]=1/2*(x[n]+x[n+1])
can be found by using a general complex exponential as an input and observing the output or response.
x[n] = Aej(⌅ˆn + ⇤) (1)
y[n] = 1
2Aej(⌅ˆn + ⇤) + 1
2Aej(⌅ˆ(n  1) + ⇤) (2)
= Aej(⌅ˆn + ⇤) 1
2
⇥
1 + ej⌅ˆ⇤
(3)
In (3) there are two terms, the original input, and a term that is a function of ⌅ˆ. This second term is the
frequency response and it is commonly denoted by H(ej⇥ˆ ).
2
H(ej⇥ˆ ) = 1
2
⇥
1 + ej⌅ˆ⇤
(4)
Once the frequency response, H(ej⇥ˆ ), has been determined, the effect of the filter on any complex exponential may be determined by evaluating H(ej⇥ˆ ) at the corresponding frequency. The output signal y[n],
will be a complex exponential whose complex amplitude has a constant magnitude and phase. The phase
describes the phase change of the complex sinusoid and the magnitude describes the gain applied to the
complex sinusoid.
For a general FIR linear time-invariant system
y[n] =
M
k=0
bkx[n  k]
the frequency response is
H(ej⇥ˆ ) =
M
k=0
bkej⇥ˆk (5)
1.1.1 MATLAB Function for Frequency Response
MATLAB has a built-in function for computing the frequency response of a discrete-time LTI system. The
following MATLAB statements show how to use freqz to compute and plot both the magnitude (absolute
value) and the phase of the frequency response of a two-point averaging system as a function of ⌅ˆ in the
range ⇥ ⌅ ⌅ˆ ⌅ ⇥:
bb = [0.5, 0.5]; %-- Filter Coefficients
ww = -pi:(pi/100):pi; %-- omega hat
H = freqz(bb, 1, ww); %<--freekz.m is an alternative
subplot(2,1,1);
plot(ww, abs(H))
subplot(2,1,2);
plot(ww, angle(H))
xlabel(’Normalized Radian Frequency’)
For FIR filters, the second argument of freqz must always be equal to 1, as we have in freqz(bb,1,ww).
3
Also, the frequency vector ww should cover an interval of length 2⇥ for ⌅ˆ, and its spacing must be fine
enough to give a smooth curve for H(ej⇥ˆ ). Note: we will always use capital H for the frequency response.
2
The notation H(ejˆ ) is used in place of the notation H(⇥ˆ) in Chapter 6 of the text for the frequency response because we will
eventually connect this notation with the z-transform, H(z), in Chapter 7. 3
If the output of the freqz function is not assigned, then plots are generated automatically; however, the magnitude is given in
decibels which is a logarithmic scale. For linear magnitude plots a separate call to plot is necessary.
2
1.2 Periodicity of the Frequency Response
The frequency responses of discrete-time filters are always periodic with period equal to 2⇥. Explain why
this is the case by stating a definition of the frequency response and then considering two input sinusoids
whose frequencies are ⌅ˆ and ⌅ˆ + 2⇥.
x1[n] = ej⌅ˆn versus x2[n] = ej(⌅ˆ + 2⇥)n
Consult Chapter 6 for a mathematical proof that the outputs from each of these signals will be identical
(basically because x1[n] is equal to x2[n].) The implication of periodicity is that a plot of H(ej⇥ˆ ) only
needs to extend over the interval ⇥ ⌅ ⌅ˆ ⌅ ⇥.
1.3 Frequency Response of the Four-Point Averager
In Chapter 6 we examined filters that average input samples over a certain interval. These filters are called
“running average” filters or “averagers” and they have the following form for the L-point averager:
y[n] = 1
L
L
1
k=0
x[n  k] (6)
(a) Use Euler’s formula and complex number manipulations to show that the frequency response for the
5-point running average operator is given by:
H(ej⇥ˆ ) = 1 + 2 cos(⌅ˆ) + 2 cos(2⌅ˆ)
5
ej2⇥ˆ (7)
(b) Use the geometric series
N
1
k=0
k = 1  N
1 
to show that the frequency response of the 5-point running average filter may also be expressed in the
form
H(ej⇥ˆ ) = ej2⇥ˆ sin(5⌅ˆ/2)
sin(⌅ˆ/2)
(c) Implement (7) directly in MATLAB. Specifically, use a vector that includes 400 samples between ⇥
and ⇥ for ⌅ˆ, and plot the function H(ej⇥ˆ ). Since the frequency response is a complex-valued quantity,
use abs() and angle() to extract the magnitude and phase of the frequency response for plotting.
Plotting the real and imaginary parts of H(ej⇥ˆ ) is not very informative.
(d) In this part, use freqz.m in MATLAB to compute H(ej⇥ˆ ) numerically (from the filter coefficients)
and plot its magnitude and phase versus ⌅ˆ. Write the appropriate MATLAB code to plot both the
magnitude and phase of H(ej⇥ˆ ). Follow the example in Section 1.1.1. The filter coefficient vector for
the 5-point averager is defined via:
bb = 1/5*ones(1,5);
Note: the function freqz(bb,1,ww) evaluates the frequency response for all frequencies in the
vector ww. It uses the summation in (5), not the formula in (7). The filter coefficients are defined in
the assignment to vector bb. How do your results compare with part (b)?
3
1.4 The MATLAB FIND Function
Often signal processing functions are performed in order to extract information that can be used to make
a decision. The decision process inevitably requires logical tests, which might be done with if-then
constructs in MATLAB. However, MATLAB permits vectorization of such tests, and the find function is
one way to do lots of tests at once. In the following example, find extracts all the numbers that “round” to
3:
xx = 1.4:0.33:5, jkl = find(round(xx)==3), xx(jkl)
The argument of the find function can be any logical expression. Notice that find returns a list of indices
where the logical condition is true. See help on relop for information on relational operators.
Now, suppose that you have a frequency response:
ww = -pi:(pi/500):pi; HH = freqz( 1/5*ones(1,5), 1, ww );
Use the find command to determine the indices where HH is zero, and then use those indices to display the
list of frequencies where HH is zero. Since there might be round-off error in calculating HH, the logical test
should probably be a test for those indices where the magnitude (absolute value in MATLAB) of HH is less
than some rather small number, e.g., 1 ⇥ 106. Compare your answer to the frequency response that your
plotted for the five-point averager in Section 1.3.
2 Warm-up
The first objective of this warm-up is to use a MATLAB GUI to demonstrate nulling. If you are working in
the ECE lab it is NOT necessary to install the GUI; otherwise, you must download the ZIP file and install it
into its own directory. This demo, dltidemo, can be downloaded from the web page:
http://users.ece.gatech.edu/mcclella/matlabGUIs/index.html
2.1 LTI Frequency Response Demo
The dltidemo illustrates the “sinusoid-IN gives sinusoid-OUT” property of LTI systems. In this demo,
you can change the amplitude, phase and frequency of an input sinusoid, x[n], and you can change the
digital filter that processes the signal. Then the GUI will show the output signal, y[n], which is also a
sinusoid (at the same frequency). Figure 1 shows the interface for the dltidemo GUI. It is possible to
see the formula for the output signal, if you click on the Theoretical Answer button located at the
bottom-middle part of the window. The digital filter can be changed by choosing different options in the
Filter Specifications box in the lower right-hand corner.
In the Warm-up, you should perform the following steps with the dltidemo GUI:
(a) Set the input to x[n] = 2.5 cos(0.08⇥(n  5))
(b) Set the digital filter to be a 11-point averager.
(c) Determine the formula for the output signal and write it in the form: y[n] = A cos(⌅ˆ0(n  nd)).
(d) Using nd for y[n] and the fact that the input signal had a peak at n = 5, determine the amount of delay
through the filter. In other words, how much has the peak of the cosine wave shifted?
Instructor Verification (separate page)
4
Figure 1: LTI demo interface.
(e) Now, determine the length of the averaging filter so that the output will be zero, i.e., y[n] = 0. Use
the GUI to show that you have the correct filter to zero the output. If the length is more than 15, you
will have to enter the “Filter Specifications” with the user Input option.
(f) When the output is zero, the filter acts as a Nulling Filter, because it eliminates the input at ⌅ˆ = 0.08⇥.
What other frequencies ⌅ˆ are also nulled? Find at least one.
Instructor Verification (separate page)
2.2 Cascading Two Systems
More complicated systems are often made up from simple building blocks. In Fig. 2, two FIR filters are
shown connected “in cascade.”
 FIR   Filter #1
FIR
Filter #2
x[n] w[n] y[n]
Figure 2: Cascade of two FIR filters.
Assume that the system in Fig. 2 is described by the two equations
w[n] =
M
⌥=0
⌥
x[n  ] (FIR FILTER #1)
y[n] = w[n]   w[n  1] (FIR FILTER #2)
5
(a) Use freqz() in MATLAB to get the frequency responses for each of the FIR filters for the case
where  = 0.81 and M = 15. Plot the magnitude and phase of the frequency response for Filter #1,
and also for Filter #2. Which one of these filters is a lowpass filter?
(b) Plot the magnitude and phase of the frequency response of the overall cascaded system.
(c) Explain how the individual frequency responses in part(a) are combined to get the overall frequency
response in part(b). Comment on the magnitude combinations as well as the phase combinations.
Instructor Verification (separate page)
2.3 Deconvolution
In the previous lab, the two filters from Section 2.2 were used in an image deblurring experiment. You
should now re-interpret how that experiment worked by explaining what happens in the frequency domain.
(a) If a single filter has a frequency response H(ej⇥ˆ ) = 1, how is the output of the filter y[n] related to
the input x[n]?
(b) Ideally, a “deconvolved” output should look exactly like the input prior to blurring. If Filter #1
(in Fig. 2) has a frequency response H1(ej⇥ˆ ), and Filter #2 is H2(ej⇥ˆ ), explain why the condition
H1(ej⇥ˆ )H2(ej⇥ˆ ) = 1 will guarantee “perfectly deconvolution.”
(c) The filters in Section 2.2 do not perform a perfect deconvolution (for the case  = 0.81 and M = 15).
Use the frequency response from Section 2.2(b) to explain deviations from a perfect result.
Instructor Verification (separate page)
3 Lab Exercises
3.1 Nulling Filters for Rejection
Nulling filters are filters that completely eliminate some frequency component. If the frequency to be eliminated is ⌅ˆ = 0 or ⌅ˆ = ⇥, then a two-point FIR filter will do the nulling. The simplest possible general
nulling filter can have as few as three coefficients. If ⌅ˆn is the desired nulling frequency, then the following
length-3 FIR filter
y[n] = x[n]  2 cos(⌅ˆn)x[n  1] + x[n  2] (8)
will have a zero in its frequency response at ⌅ˆ = ⌅ˆn. For example, a filter designed to completely eliminate
signals of the form Aej0.5n would have the following coefficients
b0 = 1, b1 = 2 cos(0.5⇥) = 0, b2 = 1.
because we would pick the desired nulling frequency to be ⌅ˆn = 0.5⇥.
(a) Design a filtering system that consists of the cascade of two FIR nulling filters that will eliminate the
following input frequencies: ⌅ˆ = 0.22⇥, and ⌅ˆ = 0.85⇥. For this part, derive the filter coefficients of
both nulling filters.
(b) Generate an input signal x[n] that is the sum of three sinusoids:
x[n] = 20 cos(0.22⇥n) + 25 cos(0.45⇥n  ⇥/3) + 20 cos(0.85⇥n  ⇥/4)
Make the input signal 150 samples long over the range 0 ⌅ n ⌅ 149.
6
(c) Use firfilt (or conv) to filter the sum of three sinusoids signal x[n] through the filters designed
in part (a). Show the MATLAB code that you wrote to implement the cascade of two FIR filters.
(d) Make a plot of the output signal—show the first 40 points. Determine (by hand) the exact mathematical formula (magnitude, phase and frequency) for the output signal for n ⇧ 5.
(e) Plot the mathematical formula in MATLAB to show that it matches the filter output from firfilt
over the range 5 ⌅ n ⌅ 40.
(f) Explain why the output signal is different for the first few points. How many “start-up” points are
found, and how is this number related to the lengths of the filters designed in part (a)? Hint: consider
the length of a single FIR filter that is equivalent to the cascade of two length-3 FIRs.
3.2 Simple Bandpass Filter Design
The L-point averaging filter is a lowpass filter. Its passband width is controlled by L, being inversely
proportional to L. In fact, you can use the GUI dltidemo to view the frequency response for different
averagers and measure the passband widths. It is also possible to create a filter whose passband is centered
around some frequency other than zero. One simple way to do this is to define the impulse response of an
L-point FIR as:
h[n] = 2
L cos(⌅ˆcn), 0 ⌅ n < L
where L is the filter length, and ⌅ˆc is the center frequency that defines the frequency location of the passband.
For example, we would pick ⌅ˆc = 0.45⇥ if we want the peak of the filter’s passband to be centered at 0.45⇥.
The bandwidth of the bandpass filter is controlled by L; the larger the value of L, the narrower the bandwidth.
This particular filter is also discussed in the section on useful filters in Chapter 7.
(a) Generate a bandpass filter that will pass a frequency component at ⌅ˆ = 0.45⇥. Make the filter length
(L) equal to 12. Since we are going to be filtering the signal defined in section 3.1(b), measure the
gain of the filter at the three frequencies of interest: ⌅ˆ = 0.22⇥, ⌅ˆ = 0.45⇥ and ⌅ˆ = 0.85⇥.
(b) The passband of the BPF filter is defined by the region of the frequency response where |H(ej⇥ˆ )| is
close to its maximum value. If we define the maximum to be Hmax, then the passband width is defined
as the length of the frequency region where the ratio |H(ej⇥ˆ )|/Hmax is greater than 1/
2 = 0.707.
Figure 3 shows how to define the passband and stopband. Note: you can use MATLAB’s find
function to locate those frequencies where the magnitude satisfies |H(ej⇥ˆ )| ⇧ 0.707Hmax.
The stopband of the BPF filter is defined by the region of the frequency response where |H(ej⇥ˆ )| is
close to zero. In this case, we will define the stopband as the region where |H(ej⇥ˆ )| is less than 10%
of the maximum.
Make a plot of the frequency response for the L = 12 bandpass filter from part (a), and determine the
passband width (at the 0.707 level). Repeat the plot for L = 24 and L = 48, so you can explain how
the width of the passband is related to filter length L, i.e., what happens when L is doubled or halved.
(c) Comment on the selectivity of the L = 12 bandpass filter. In other words, which frequencies are
“passed by the filter.” Use the frequency response to explain how the filter can pass one component at
⌅ˆ = 0.45⇥, while reducing or rejecting the others at ⌅ˆ = 0.22⇥ and ⌅ˆ = 0.85⇥.
(d) Generate a bandpass filter that will pass the frequency component at ⌅ˆ = 0.45⇥, but now make the
filter length (L) long enough so that it will also greatly reduce frequency components at (or near)
⌅ˆ = 0.22⇥ and ⌅ˆ = 0.85⇥. Determine the smallest value of L so that
7
0 0.5 1 1.5 2 2.5 3
0
0.2
0.4
0.6
0.8
1
Frequency (radians)
Magnitude
BANDPASS FILTER (centered at 0.4π)
PASSBAND
STOPBAND STOPBAND
Figure 3: Passband and Stopband for a typical FIR bandpass filter. In this case, the maximum value is 1,
the passband is the region where the frequency response is greater than 1/
2 = 0.707, and the stopband is
defined as the region where the frequency response is less than 25% of the maximum.
• Any frequency component satisfying |⌅ˆ| ⌅ 0.22⇥ will be reduced by a factor of 10 or more.4
• Any frequency component satisfying 0.85⇥ ⌅ |⌅ˆ| ⌅ ⇥ will be reduced by a factor of 10 or
more.
This can be done by making the passband width very small.
(e) Use the filter from the previous part to filter the “sum of 3 sinusoids” signal from Section 3.1. Make a
plot of 100 points of the input and output signals, and explain how the filter has reduced or removed
two of the three sinusoidal components.
(f) Make a plot of the frequency response (magnitude only) for the filter from part (d), and explain how
H(ej⇥ˆ ) can be used to determine the relative size of each sinusoidal component in the output signal.
In other words, connect a mathematical description of the output signal to the values that can be
obtained from the frequency response plot.
4
For example, the input amplitude of the 0.22 component is 20, so its output amplitude must be less than 2.
8
ECES-352