Description
EE 453: Lab 2
Introduction to Filter Design and Image Processing
March 15, 2021
Overview
This lab will be based around filter design and image processing. You will design
and find the characteristics of the Butterworth, Chebyshev, and Elliptic filters.
The final part of the lab will be a connection between the Radon Transform
used in medical imaging, and the Fourier Transform.
1
1
Simulink Filter Design
For the first part of this experiment, you will learn how to implement Butterworth, Chebyshev, and Elliptic filters and verify their properties.
1. Open up Simulink by typing simulink into the command prompt as usual,
and open a blank model
2. Setup the signal I/O: DSP System toolbox → Sources → From Multimedia File and DSP System toolbox → Sinks → Audio Device
Writer. You will be using sinusoids at three main frequencies for this part
of the lab: 500 Hz, 1000 Hz, and 3500 Hz, all sampled at Fs = 8000Hz.
Create wav files with each of these frequencies, and select one to be the
input signal (it doesn’t matter which one for now). You may just want to
confirm everything is working by listening to the tone
3. Next we are going to place the Digital Filter Design tool to start modifying
signals, and the spectrum analyzer to read the results. The Digital Filter
Design tool can be found in DSP System toolbox → Filtering →
Filter Implementations → Digital Filter Design. Connect the signal
input (from multimedia) to the input of the Digital Filter Design tool, and
the output of the Digital Filter Design tool to the signal output (device
writer).
4. The spectrum analyzer is going to be found in DSP System Toolbox
→ Sinks → Spectrum Analyzer. You are going to want two spectrum
analyzers. Connect one in parallel to the signal input (from multimedia).
This is your reference signal. Connect the second in parallel to the signal
output (device writer). The difference between the two spectrum analyzers
will show the effect of the filter. Doing it this way avoids the need to have
normalized input signals.
5. The options for the digital filter design tool are as follows:
• Design Method: Bottom Left
• Filter Order: Middle
• Filter Specifications: Left
• Ripple: Right (Apass = 6 means a 6 dB ripple in passband)
• All dB measurements are in voltage (20log10)
Butterworth LPF
6. The first filter we are going to design will be a lowpass Butterworth filter. The Butterworth filter is associated with having as flat a frequency
response as possible in the passband. This would be ideal when the frequency components you want to filter out are far away from the frequency
components of interest. For example, a speech audio file is corrupted
2
with the noise of a jackhammer outside (high frequency). By using a
Butterworth filter the speech is unaffected as much as possible, while the
jackhammer is suppressed.
The frequency response for a low pass Butterworth filter is given as:
H(ejω ) = q
1
jω 2N
1 + ( jω
)
c
• N is the filter order
• ωc is the center frequency of the filter
Answer Q1 on your worksheet
7. Set the Simulink Filter to model a 4th order Butterworth low pass filter,
with a cutoff frequency of 1 kHz.
Answer Q2 on your worksheet
Call your TA over to verify that this program is working. Your
TA will sign your lab worksheet.
Answer Q3 on your worksheet
3
Chebyshev HPF
The second filter you are going to design is a Chebyshev High Pass Filter.
The Chebyshev filters are characterized by the equiripple response in the
passband. This means that the error between the original and filtered
signal can be bounded by some maximum error. The filter response will
vary smoothly within that bound, and always go from bound to bound.
The filter is guaranteed not to exceed this error. Therefore, a Chebyshev
high pass filter can be used just like a Butterworth filter, if a small amount
of error (ripple) is acceptable in the passband. A Chebyshev filter will
have a smaller transition band between the passband and the stopband.
Therefore, lets say we have recorded an electronic piece of music, but the
60-cycle hum from the wall outlet adds unwanted noise at 60 Hz. Music
of interest is too close to the 60 Hz for a Butterworth to be used, but we
can bound the ripple of the Chebyshev filter to provide an approximation
close to the original piece of music, without the 60-cycle hum.
8. The magnitude response for a low pass Chebyshev filter is given as:
1
|H(ejω )| = q
1 + ε2 TN2 ( ωω0 )
• N is the filter order
• TN (x) ≈ cos(N cos−1 (x)) is the N th degree Chebyshev polynomial
• ε is the ripple parameter. If δ is the amount of ripple in dB that is
desired in the passband, then ε2 = 100.1δ − 1
• ω0 is the center frequency of the filter, but actually not where the -3
dB point is. That point is at ω3dB = ω0 cosh( N1 cosh−1 ( 1ε )
The transformation from a low pass to a high pass filter is:
ω
−ω0
⇒
ω0
ω
Answer Q4 on your worksheet
9. Set the Simulink Filter to model a 4th order Chebyshev Type 1 high pass
filter, with a cutoff frequency of 1 kHz, and a ripple of 6 dB.
Answer Q5 on your worksheet
Q6 on your worksheet
4
Elliptic BPF
The benefit of an elliptic filter is it is a generalization of the Butterworth
and Chebyshev filters. By changing parameters, the Elliptic filter can act
like a Butterworth. Using different parameters it has the properties of a
Chebyshev filter. It has the added advantage of adding ripple in both the
passband and stopband. While this may not seem ideal, the more ripple
you have the smaller the transition band. Therefore, of all three filters the
Elliptic filter has the smallest transition band, but largest ripple.
10. The magnitude response for a low pass elliptic filter is given as:
|H(ejω )| = q
1
1+
2 (ξ, ω )
ε2 RN
ω0
• N is the filter order
• RN (x) ≈
Πi (x−xi,z )
Πi (x−xi,p )
is the elliptic rational function
• ε and ξ are ripple parameters
• As ξ → ∞ the Elliptic filter approaches a Chebyshev filter
• As ξ → ∞, ε → 0 the Elliptic filter approaches a Butterworth filter
The transformation from a low pass to a band pass filter is:
ω
1 ω
ω0
⇒ (
−
)
ω0
∆ ω0
ω
1
Bandwidth
=
∆
ω0
Answer Q7 on your worksheet
11. Set the Simulink Filter to model a 4th order Elliptic band pass filter, with
a first pass band Fpass,1 = 1000 to Fpass,2 = 3000. Use a passband ripple
of 3 dB, and a stopband ripple of 9 dB.
Answer Q8 on your worksheet
Call your TA over to verify that this program is working. Your
TA will sign your lab worksheet.
5
Radon Transform
The Radon transform is the key component of computed tomography
(CAT) scans for X-ray imaging in hospitals, as well as a host of other
biomedical imaging modalities. The radon transform can be derived from
the attenuation of X-rays within a biologic tissue. As electromagnetic
waves, X-rays experience an exponential decay when traveling through a
lossy medium. Therefore, a differential equation describing the attenuation would be
dI
= −µ(x)I
dx
Where I is the intensity of the wave, and µ(x) is the attenuation coefficient.
R
This ordinary differential equation has the solution I1 = I0 exp(− L µ(x)dx).
R
Taking logarithms, this can be written as L µ(x)dx = log( II10 ). Therefore,
by taking a measurement of the intensity I1 at a sensor, that is equivalent
to the integral of the absorption over the line connecting sensor I0 and I1 .
The Radon transform in general is equal to the line integral of a function,
with the specific line being a parameter of the transform.
Z
R{f (x)} = R(l) = f (x)ds
l
In biomedical imaging, the data received by sensors is R(l). The goal is to
recover f (x), which holds information about the tissue, from R(l). This
is accomplished via inverting the Radon transform. In order to invert the
Radon transform, an extremely useful result called the Projection-Slice
Theorem is used. The Projection-Slice Theorem states that if you take
the Fourier transform of the Radon Transform, that will be equivalent to
the two dimensional Fourier transform of the function evaluated along the
same line that the Radon transform was taken over.
F1 R = F2 (l)
This extremely important result means that given Radon transform (sensor) data, the f (x) can be recovered by extremely efficient fast Fourier
transform algorithms, namely one forward FFT, and then a two dimensional FFT. Unfortunately, because discrete measurements are taken, and
not all possible paths are reconstructed, the reconstruction of the function
is often not ideal. Therefore, filters are often applied during the inversion
process, leading to techniques known as filtered back projection. We will
not explore filtered back projection algorithms in this lab, but they do use
windowing techniques. Windowing techniques will be the theme of Lab 3,
which we will cover filtered back projection as well as other processes that
use windowing in digital signal processing.
12. Our exploration of the radon transform will begin with a brain phantom
often used in medical imaging. It can be loaded by the MATLAB command phantom(N), where N is the size of the image (N = 256 is a good
start).
6
13. The Radon transform is taken over lines through the origin, which can
be parameterized by angles between 0 and 180 (MATLAB’s transform
expects angles in degrees). It is suggested to have many angles, so as to
recreate the image as close as possible (between 90 and 180 would be a
good start). Important: The angles must always lie between 0 and 180,
so for instance, 90 angles could be given by linspace(0,180,90).
14. The Radon transform in MATLAB is given by R = radon(I,th). The
inverse likewise is Irec = iradon(R,th). The images can be viewed using
imagesc(R) or imagesc(Irec). It is advised to use the gray colormap.
15. Add noise to the signal in the spatial domain using In = imnoise(I). Separately, add noise to the radon transformed signal, Rn = imnoise(R).
Keep in mind, in real sensing applications, all the measurement and environmental noise on the sensors adds noise to the Radon domain. This is
why signal processing is extremely important in biomedical sensing applications.
16. Compare the reconstruction when noise is added to the signal in the spatial
domain vs the radon transformed signal. Use PSNR (Peak Signal to Noise
ratio, MATLAB function psnr ) to have a quantitative measurement of
the error. Notice: The reconstructed signal may have to be truncated to
match the size of the original image for comparison.
Answer Q9 and Q10 on your worksheet
Call your TA over to verify that this program is working. Your
TA will sign your lab worksheet.
7
EE 453: Lab 2 Questions
Introduction to Filtering and Image Processing
March 15, 2021
Name:
Lab Partners:
1. What is the Magnitude Response of the Butterworth LPF
2. Calculate analytically what the expected attenuation is for each of the
following frequencies, and record what your measurements indicate when
you actually pass a sinusoid through the filter. Remember to take the
difference between two spectrum analyzers to find the attenuation.
Frequencies (Hz)
500
Analytic Attenuation (dB)
1000
3500
1
Measured Attenuation (dB)
TA Check Butterworth:
3. Qualitatively, what is the difference when the filter order is changed to
2nd order (use the display in the Digital Filter Model). How about when
it is changed to 6th order?
4. What is the Magnitude Response of the Chebyshev HPF
5. Calculate analytically what the expected attenuation is for each of the
following frequencies, and record what your measurements indicate when
you actually pass a sinusoid through the filter. Remember to take the
difference between two spectrum analyzers to find the attenuation.
Frequencies (Hz)
500
Analytic Attenuation (dB)
1000
3500
2
Measured Attenuation (dB)
6. Qualitatively, what is the difference when the filter order is changed to
2nd order (use the display in the Digital Filter Model). How about when
it is changed to 6th order?
7. What is the Magnitude Response of the Elliptic BPF
8. Qualitatively, what is the difference when the filter order is changed to
2nd order (use the display in the Digital Filter Model). How about when
it is changed to 6th order?
3
TA Check Elliptic:
9. Add different types of noise (Poisson, Gaussian, salt and pepper). Which
types of noise effect the reconstructed image the most? Which domain
(spatial or radon) is more sensitive to noise when it comes to the reconstructed image?
10. Apply Gaussian blurring (imgaussfilt), or other methods you prefer, to
denoise the the signal. Which technique is more effective (if you tried several)? Which domain saw the most improvement in reconstructed image
quality via denoising?
4
TA Check Radon:
5
From Kangrui Han to Everyone:
how to change the parameters of Filter Designer?
From Christopher D Fadden to Everyone
you double click to bring up the menu
From Kangrui Han to Everyone:
I couldn’t save them successfully
From Christopher D Fadden to Everyone:
at the bottom right
there is, I believe something like implement
if you want to share your screen I can find it for you
From Kangrui Han to Everyone:
got it, thank you
From Christopher D Fadden to Everyone
sure
From Kangrui Han to Everyone:
For the step 7, the order is 4, and fc = 1000hz, is Fs
equal to 8000?
From Justin Henry to Everyone
That’s what I did
From Christopher D Fadden to Everyone:
yes
From Kangrui Han to Everyone:
ok thanks
Purchase answer to see full
attachment