## Description

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