Lab  #7:  FIR    Filtering    of    Images    

ECES-352

Lab    #7:        FIR    Filtering    of    Images    
(Lab Report Du e at beginning of next lab)
This    is    the    official    Lab    #7    description.    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.    
Formal    Lab    Report:        You    must    write    a    formal    lab    report    that    describes    your    
approach    to    music    synthesis    (Section    4).    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    learn    how    to    implement    FIR    filters    in    MATLAB,    and    then    
study    the    response    of    FIR    filters    to    various    signals,    including    images.    As    a    result,    you    
should    learn    how    filters    can    create    effects    such    as    blurring    and    ghosts.    In    addition,    
we    will    use    FIR    filters    to    study    the    convolution    operation    and    properties    such    as    
linearity    and    time-invariance.
In    the    experiments    of    this    lab,    you    will    use    firfilt(    ),    or    conv(),    to    implement    1-D    
filters    and    conv2()    to    implement    two-dimensional    (2-D)    filters.    In    this    lab    the    2-D    
filtering    operation    actually    consists    of    1-D    filters    applied    to    all    the    rows    of    the    image    
and    then    to    all    of    the    columns.    
1.1 Two GUIs
This    lab    involves    the    use    of    two    MATLAB    GUIs:    one    for    sampling    and    aliasing    and    
one    for    convolution.
1. con2dis: GUI for sampling and aliasing. An input sinusoid and its spectrum is tracked
through A/D and D/A converters.
2. dconvdemo: GUI for discrete-time convolution. This is exactly the same as the MATLAB
functions conv() and firfilt() used to implement FIR filters.
Both of these demos can be downloaded from WebCT for ECE2025. From the Homepage, click ”Extra
M-files for Labs”. Then click ”Educational Matlab GUIs”. Then go to the two demos (continuous-discrete
sampling demo and discrete convolution demo) and click on the download statement.
1.2 Overview of Filtering
For this lab, we will define an FIR filter as a discrete-time system that converts an input signal x[n] into an
output signal y[n] by means of the weighted summation:
y[n] = ⇥
M
k=0
bk x[n  k] (1)
Equation (1) gives a rule for computing the nth value of the output sequence from certain values of the input
sequence. The filter coefficients {bk} are constants that define the filter’s behavior. As an example, consider
the system for which the output values are given by
y[n] = 1
3x[n] + 1
3x[n  1] + 1
3x[n  2] (2)
= 1
3 {x[n] + x[n  1] + x[n  2]}
This equation states that the nth value of the output sequence is the average of the nth value of the input
sequence x[n] and the two preceding values, x[n  1] and x[n  2]. For this example the bk’s are b0 = 1
3 ,
b1 = 1
3 , and b2 = 1
3 .
MATLAB has built-in functions, conv( ) and filter( ), for implementing this operation in (1),
but we have also supplied another M-file firfilt( ) for the special case of FIR filtering. The function
⇤⇥⌅
CD-ROM
firfilt.m filter implements a wider class of filters than just the FIR case. Technically speaking, the both the conv
and the firfilt function implement the operation called convolution. The following MATLAB statements
implement the three-point averaging system of (2):
nn = 0:99; %<--Time indices
xx = cos( 0.08*pi*nn ); %<--Input signal
bb = [1/3 1/3 1/3]; %<--Filter coefficients
yy = firfilt(bb, xx); %<--Compute the output
In this case, the input signal xx is a vector containing a cosine function. In general, the vector bb contains the filter coefficients {bk} needed in (1). These are loaded into the bb vector in the following way:
bb = [b0, b1, b2, ... , bM].
In MATLAB, all sequences have finite length because they are stored in vectors. If the input signal has,
for example, L samples, we would normally only store the L samples in a vector, and would assume that
x[n] = 0 for n outside the interval of L samples; i.e., we do not have to store any zero samples unless it suits
our purposes. If we process a finite-length signal through (1), then the output sequence y[n] will be longer
than x[n] by M samples. Whenever firfilt( ) implements (1), we will find that
length(yy) = length(xx)+length(bb)-1
In the experiments of this lab, you will use firfilt( ) to implement FIR filters and begin to understand
how the filter coefficients define a digital filtering algorithm. In addition, this lab will introduce examples to
show how a filter reacts to different frequency components in the input.
2
1.3 Pre-Lab: Run the GUIs
The first objective of this lab is to demonstrate that you understand the concepts behind the two GUIs. First
of all, you must download the ZIP files for each and install them. Each ZIP file installs as a directory
containing a number of files. You can put the GUIs on the matlabpath, or you can run the GUIs from
their home directories.
1.4 Sampling and Aliasing Demo
In this demo, you can change the frequency of an input signal that is a sinusoid, and you can change the
sampling frequency. The GUI will show the sampled signal, x[n], its spectrum, and also the reconstructed
output signal, y(t) with its spectrum. Figure 1 shows the interface for the con2dis GUI. In order to see
the entire GUI, you must select Show All Plots under the Plot Options menu.
Figure 1: Sampling/reconstruction demo interface.
In the pre-Lab, you should perform the following steps with the con2dis GUI:
(a) Set the input to x(t) = cos(40⇥t)
(b) Set the sampling rate to fs = 24 samples/sec.
(c) Determine the locations of the spectrum lines for the discrete-time signal, x[n], found in the middle
panels. Click the Radian button to change the axis to from ˆf to ⇤ˆ.
(d) Determine the formula for the output signal, y(t) shown in the rightmost panels. What is the output
frequency in Hz?
3
1.5 Discrete-Time Convolution Demo
In this demo, you can select an input signal x[n], as well as the impulse response of the filter h[n]. Then
the demo shows the “flipping and shifting” used when a convolution is computed. This corresponds to the
sliding window of the FIR filter. Figure 2 shows the interface for the dconvdemo GUI.
Figure 2: Interface for discrete-time convolution GUI.
In the pre-lab, you should perform the following steps with the dconvdemo GUI.
(a) Click on the Get x[n] button and set the input to a finite-length pulse: x[n] = (u[n]  u[n  10]).
(b) Set the filter to a three-point averager by using the Get h[n] button to create the correct impulse
response for the three-point averager. Remember that the impulse response is identical to the bk’s for
an FIR filter. Also, the GUI allows you to modify the length and values of the pulse.
(c) Use the GUI to produce the output signal.
(d) When you move the mouse pointer over the index “n” below the signal plot and do a click-hold, you
will get a hand tool that allows you to move the “n”-pointer. By moving the pointer horizontally you
can observe the sliding window action of convolution. You can even move the index beyond the limits
of the window and the plot will scroll over to align with “n.”
1.6 Filtering via Convolution
You can perform the same convolution as done by the dconvdemo GUI by using the MATLAB function
firfilt, or conv. For ECE-2025, the preferred function is firfilt.
4
(a) For the Pre-Lab, you should do the filtering with a 3-point averager. The filter coefficient vector for
the 3-point averager is defined via:
bb = 1/3*ones(1,3);
Use firfilt to process an input signal that is a length-10 pulse:
x[n] =
 1 for n = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
0 elsewhere
NOTE: in MATLAB indexing can be confusing. Our mathematical signal definitions start at n =
0, but MATLAB starts its indexing at “1”. Nevertheless, we can ignore the difference and pretend
that MATLAB is indexing from zero, as long as we don’t try to write x[0] in MATLAB. For this
experiment, generate the length-10 pulse and put it inside of a longer vector with the statment xx
= [ones(1,10),zeros(1,5)]. This produces a vector of length 15, which has 5 extra zero
samples appended.
(b) To illustrate the filtering action of the 3-point averager, it is informative to make a plot of the input
signal and output signals together. Since x[n] and y[n] are discrete-time signals, a stem plot is
needed. One way to put the plots together is to use subplot(2,1,*) to make a two-panel display:
nn = first:last; %--- use first=1 and last=length(xx)
subplot(2,1,1);
stem(nn-1,xx(nn))
subplot(2,1,2);
stem(nn-1,yy(nn),’filled’) %--Make black dots
xlabel(’Time Index (n)’)
This code assumes that the output from firfilt is called yy. Try the plot with first equal to
the beginning index of the input signal, and last chosen to be the last index of the input. In other
words, the plotting range for both signals will be equal to the length of the input signal, even though
the output signal is longer. Notice that using nn-1 in the call to stem( ) causes the x-axis to start
at zero in the plot.
(c) Explain the filtering action of the 3-point averager by comparing the plots in the previous part. This
filter might be called a “smoothing” filter. Note how the transitions in x[n] from zero to one, and from
one back to zero, have been “smoothed.”
2 Warm-up
2.1 Sampling and Aliasing
Use the con2dis GUI to do the following problem:
(a) Input frequency is 12 Hz.
(b) Sampling frequency is 15 Hz.
(c) Determine the frequency of the reconstructed output signal
(d) Determine the locations in ⇤ˆ of the lines in the spectrum of the discrete-time signal. Give numerical
values.
Instructor Verification (separate page)
(e) Change the sampling frequency to 12 Hz, and explain the appearance of the output signal.
5
2.2 Discrete-Time Convolution
In this section, you will generate filtering results needed in a later section. Use the discrete-time convolution
GUI, dconvdemo, to do the following:
(a) Set the input signal to be x[n] = (0.92)n (u[n]  u[n  10]). Use the “Exponential” signal type
within Get x[n].
(b) Set the impulse response to be h[n] = [n]  0.92[n  1]. Once again, use the “Exponential” signal
type within Get h[n].
(c) Illustrate the output signal y[n] and explain why it is zero for almost all points. Compute the numerical
value of the last point in y[n], i.e., the one that is negative and non-zero.
Instructor Verification (separate page)
2.3 Loading Data
In order to exercise the basic filtering function firfilt, we will use some “real” data. In MATLAB you
can load data from a file called lab6dat.mat file by using the load command as follows:
load lab6dat
The data file lab6dat.mat contains two filters and three signals, stored as separate MATLAB variables:
⇤⇥⌅
CD-ROM
lab6dat.mat
x1: a stair-step signal such as one might find in one sampled scan line from a TV test pattern
image.
xtv: an actual scan line from a digital image.
x2: a speech waveform (“oak is strong”) sampled at fs = 8000 samples/second.
h1: the coefficients for a FIR discrete-time filter of the form of (1).
h2: coefficients for a second FIR filter.
For this lab we will only use the signal x2. After loading the data, use the whos function to verify that all
five vectors are in your MATLAB workspace.
2.4 Filtering Images: 2-D Convolution
One-dimensional FIR filters, such as running averagers and first-difference filters, can be applied to onedimensional signals such as speech or music. These same filters can be applied to images if we regard each
row (or column) of the image as a one-dimensional signal. For example, the 50th row of an image is the
N-point sequence xx[50,n] for 1 ⌅ n ⌅ N, so we can filter this sequence with a 1-D filter using the
conv or firfilt operator.
One objective of this lab is to show how simple 2-D filtering can be accomplished with 1-D row and
column filters. It might be tempting to use a for loop to write an M-file that would filter all the rows. For a
first-difference filter, this would create a new image made up of the filtered rows:
y1[m, n] = x[m, n]  x[m, n  1].
However, this image y1[m, n] would only be filtered in the horizontal direction. Filtering the columns would
require another for loop, and finally you would have the completely filtered image:
y2[m, n] = y1[m, n]  y1[m  1, n]
6
In this case, the image y2[m, n] has been filtered in both directions by a first-difference filter.
These filtering operations involve a lot of conv calculations, so the process can be slow. Fortunately,
MATLAB has a built-in function conv2( ) that will do this with a single call. It performs a more general
filtering operation than row/column filtering, but since it can do these simple 1-D operations it will be very
helpful in this lab.
(a) Load in the image echart.mat with the load command (it will create the variable echart whose
size is 257 ⇥ 256). We can filter all the rows of the image at once with the conv2( ) function. To
filter the image in the horizontal direction using a second-difference filter, we form a row vector of
filter coefficients and use the following MATLAB statements:
bdiffh = [0.25, -0.5, 0.25];
yy1 = conv2(echart, bdiffh);
In other words, the filter coefficients bdiffh for the first-difference filter are stored in a row vector
and will cause conv2( ) to filter all rows in the horizontal direction. Display the input image
echart and the output image yy1 on the screen at the same time. Compare the two images and give
a qualitative description of what you see.
(b) Now filter the “eye-chart” image echart in the vertical direction with a second-difference filter
to produce the image yy2. This is done by calling yy2 = conv2(echart,bdiffh’) with a
column vector of filter coefficients. Display the image yy2 on the screen and describe in words how
the output image compares to the input.
Instructor Verification (separate page)
3 Lab: FIR Filters
In the following sections we will study how a filter can produce the following special effects:
1. Echo: FIR filters can produce echoes and reverberations because the filtering formula (1) contains
delay terms. In an image, such phenomena would be called “ghosts.”
2. Deconvolution: One FIR filter can (approximately) undo the effects of another—we will investigate a
cascade of two FIR filters that distort and then restore an image. This process is called deconvolution.
3.1 Deconvolution Experiment for 1-D Filters
Use the function firfilt( ) to implement the following FIR filter
w[n] = x[n]  0.92x[n  1] (3)
on the input signal x[n] defined via the MATLAB statement: xx = 256*(rem(0:100,30)10); In
MATLAB you must define the vector of filter coefficients bb needed in firfilt.
(a) Plot both the input and output waveforms x[n] and w[n] on the same figure, using subplot. Make
the discrete-time signal plots with MATLAB’s stem function, but restrict the horizontal axis to the
range 0 ⌅ n ⌅ 90. Explain why the output appears the way it does by figuring out (mathematically)
the effect of the filter coefficients in (3).
(b) Note that w[n] and x[n] are not the same length. Determine the length of the filtered signal w[n], and
explain how its length is related to the length of x[n] and the length of the FIR filter. (If you need a
hint refer to Section 1.2.)
7
3.1.1 Restoration Filter
The following FIR filter
y[n] = ⇥
M
⇥=0
r⇥
w[n  ⌦] (FIR FILTER-2)
can be use to undo the effects of the FIR filter in the previous section (see the block diagram in Fig. 3). It
performs restoration, but it only does this approximately. Use the following steps to show how well it works
when r = 0.92 and M = 15.
(a) Process the signal w[n] from (3) with FILTER-2 to obtain the output signal y[n].
(b) Make stem plots of w[n] and y[n] using a time-index axis n that is the same for both signals. Put the
stem plots in the same window for comparison—using a two-panel subplot.
(c) Since the objective of the restoration filter is to produce a y[n] that is almost identical to x[n], make a
plot of the error between x[n] and y[n] over the range 0 ⌅ n ⌅ 90.
3.1.2 Worst-Case Error
(a) Evaluate the worst-case error by doing the following: find the maximum of the difference between
y[n] and x[n] in the range 0 ⌅ n ⌅ 90.
(b) What does the error plot and worst case error tell you about the quality of the restoration of x[n]?
What parameter of FILTER-2 would you change to improve the quality of the restoration? How small
do you think the worst case error has to be so that it cannot be seen on a plot?
3.1.3 An Echo Filter
The following FIR filter can be interpreted as an echo filter.
y1[n] = x1[n] + r x1[n  P] (4)
Explain why this is a valid interpretation by working out the following:
(a) You have an audio signal sampled at fs = 8000 Hz and you would like to add a delayed version of
the signal to simulate an echo. The time delay of the echo should be 0.2 seconds, and the strength of
the echo should be 90% percent of the original. Determine the values of r and P in (4); make P an
integer.
(b) Describe the filter coefficients of this FIR filter, and determine its length.
(c) Implement the echo filter in (4) with the values of r and P determined in part (a). Use the speech
signal in the vector x2 found in the file lab6dat.mat. Listen to the result to verify that you have
produced an audible echo.
8
3.2 Cascading Two Systems
More complicated systems are often made up from simple building blocks. In the system of Fig. 3 two FIR
filters are connected “in cascade.” For this section, assume that the the filters in Fig. 3 are described by the
two equations:
w[n] = x[n]  q x[n  1] (FIR FILTER-1)
y[n] = ⇥
M
⇥=0
r⇥
w[n  ⌦] (FIR FILTER-2)
 FIR   FILTER-1
FIR
FILTER-2
x[n] w[n] y[n]
Figure 3: Cascading two FIR filters: the second filter attempts to “deconvolve” the distortion introduced by
the first.
3.2.1 Overall Impulse Response
(a) Implement the system in Fig. 3 using MATLAB to get the impulse response of the overall cascaded
system for the case where q = 0.92, r = 0.92 and M = 15. Use two calls to firfilt(). Plot the
impulse response of the overall cascaded system.
(b) Work out the impulse response h[n] of the cascaded system by hand to verify that your MATLAB
result in part (a) is correct. (Hint: consult Problem 6.7 of Problem Set #6.)
(c) In a deconvolution application, the second system (FIR FILTER-2) tries to undo the convolutional effect of the first. Perfect deconvolution would require that the cascade combination of the two systems
be equivalent to the identity system: y[n] = x[n]. If the impulse responses of the two systems are
h1[n] and h2[n], state the condition on h1[n] ⇤ h2[n] to achieve perfect deconvolution.1
3.2.2 Distorting and Restoring Images
If we pick q to be a little less than 1.0, then the first system (FIR FILTER-1) will cause distortion when
applied to the rows and columns of an image. The objective in this section is to show that we can use the
second system (FIR FILTER-2) to undo this distortion (more or less). Since FIR FILTER-2 will try to undo
the convolutional effect of the first, it acts as a deconvolution operator.
(a) Load in the image echart.mat with the load command. It creates a matrix called echart.
⇤⇥⌅
CD-ROM
echart.mat
(b) Pick q = 0.92 in FILTER-1 and filter the image echart in both directions: apply FILTER-1 along the
horizontal direction and then filter the resulting image along the vertical direction also with FILTER-1.
Call the result ech92.
1
Note: the cascade of FIR FILTER-1 and FILTER-2 does not perform perfect deconvolution.
9
(c) Deconvolve ech92 with FIR FILTER-2, choosing M = 15 and r = 0.92. Describe the visual
appearance of the output, and explain its features by invoking your mathematical understanding of the
cascade filtering process. Explain why you see “ghosts” in the output image, and use some previous
calculations to determine how big the ghosts (or echoes) are, and where they are located. Evaluate the
worst-case error in order to say how big the ghosts are relative to “black-white” transitions which are
0 to 255.
3.2.3 A Second Restoration Experiment
(a) Now try to deconvolve ech92 with several different FIR filters for FILTER-2. You should set r =
0.92 and try several values for M such as 10, 15 and 30. Pick the best result and explain why it
is the best. Describe the visual appearance of the output, and explain its features by invoking your
mathematical understanding of the cascade filtering process. HINT: determine the impulse response
of the cascaded system and relate it to the visual appearance of the output image.
Hint: you can use dconvdemo to generate the impulse responses of the cascaded systems, like you
did in the Warm-up.
(b) Furthermore, when you consider that a gray-scale display has 256 levels, how large is the worst-case
error (from the previous part) in terms of number of gray levels? Do this calculation for each of the
three filters in part (a). Think about the following question: “Can your eyes perceive a gray scale
change of one level, i.e., one part in 256?”
Include all images and plots for the previous two parts to support your discussions in the lab report.
10
Lab    #7
ECES-352
Instructor    Verification    Sheet
For each verification, be prepared to explain your answer and respond to other related questions
that the lab TA’s or professors might ask. Turn this page in at the end of your lab period.
Name:________________________________ Date:________________________________
Part    2.3    (a),(b)    Process    the    input    image    echart 2-D    filter    that    filters    in    both    the    
horizontal    and    vertical    directions    with    a    first    difference    filter.    Explain    how    the    filter    
changes    the    “image    signal.”