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

Order your essay today and save **15%** with the discount code: VACCINE