mat lab report

First Download to your desktop the following images that you will find in the lab folder on Blackboard:

ImageDec03.jpg’,

image.jpg,

tennis-ball.jpg

tennis-ball2.png

tennis-ball3.jpg

tennis-ball4.jpg

tennis-ball5.png

blackwhitephoto.jpg

Then upload the above image files to Matlab by clicking on Home, Upload, filename from your desktop. All these files should appear now in the list on the left window “Current Folder”

Create a new script file named Lab3_part2_edge_detection

Enter the following commands in the script file:

% Detects the edges of an image by using the Matlab 2D

% convolution function that acts on all rows or on all columns of the

% image array

subplot(3,2,1);

I = imread('ImageDec03.jpg');

I = I(:,:,1); % This reduce the image dimension

imshow(I); % Display origina limage imageDec03

subplot(3,2,2)

bdiffh = [1 -1]; % Fist difference order FIR filter coefficents

yy1 = conv2(I,bdiffh); % Matlab 2D concvolution function that acts on all rows

imshow(yy1) % Display the row-filtered image

subplot(3,2,3)

yy2 = conv2(I,bdiffh'); % Matlab 2D convolution that acts on all columns

imshow(yy2) % Display the column-filtered image

Note that conv2 is used twice, to perform convolution on all rows and all columns of the image array (image Dec03) respectively and display the results.

Repeat the above steps for each of the other images by replacing in the above script file the command:

I = imread(‘ImageDec03.jpg');

successively with:

I = imread(‘image.jpg’);

I = imread(‘tennis-ball.jpg’);

I = imread('tennis-ball2.png');

I = imread('tennis-ball3.jpg');

I = imread('tennis-ball4.jpg');

I = imread('tennis-ball5.png');

I = imread(‘blackwhitephoto.jpg');

Obtain screenshots for each image simulations as you did above.

Part 3: Design of FIR filters in the frequency domain

Matlab has the following function to calculate the frequency response of a FIR filter and plot its magnitude and phase vs. normalized frequency:

Freqz(bb, aa) where bb = [b0, b1 …] is the array of filter coefficients and aa is 1 for FIR filters.

Design of a low pass filter by using an L-points running average filter:

Both in the lecture and textbook, a L-point running sum lowpass filter and a bandpass filter were discussed. In this section ,you will design this filter with Matlab and study how the order affects the frequency response and the pole-zero plots. Matlab uses the zplane () function to compute and plot the pole-zero plots. The real (x) and imaginary (y) parts can be obtained by hovering the cursor with the mouse over their symbols (x: poles, o: zeros).

Create a new script file called Lab3_Part3_1

Write the following Matlab commands:

% Design of a lowpass filer using an L-point running sum

M = 10; %--Filter order

bb = ones(1,M); %--Filter coefficients

ww = -pi:(pi/500):pi; %-- omega hat frequency vector

H = freqz(bb, 1, ww); %<--Matlab function for frequency response

subplot(2,1,1);

plot(ww, abs(H)), grid on

subplot(2,1,2);

plot(ww, angle(H)), grid on

xlabel('Normalized Radian Frequency')

zplane(bb,1); % pole-zero plots

Save (click on edit then save as) and run the script file and get a screenshot of the script file and plots. This should give the results of the 11-point running average filer discussed in class and in your textbook. Are the results as discussed?

Now practice with increasing the filter order M by repeat the above steps and changing the order to M = 15, 20, and 30. What happens on the pole-zero plot and the frequency response. Make sure you get a screenshot of all the results (script file and plots). Write the real and imaginary parts of all the zeros and poles. Do the zeros lie on the unit circle as discussed in class?

Design of a bandpass filter using an 11-points running average filter:

In class and the textbook (chapt. 7) it was shown that a bandpass filter can be designed from an L-point running sum filter.

Create a new script file called Lab3_Part3_2

Write the following Matlab commands:

% Bandpass filter obtained form 11 points running average filter

M = 10; %--Filter order

%bb = [1 -1]; %-- Filter Coefficients

kk = [1 2 3 4 5 6 7 8 9 10]; % for a general M order write a similar array up to M

w0 = pi/3; %—Center of the bandpass

bb = cos(w0*kk); % New filter coefficients that rotate the zeros by w0

ww = -pi:(pi/500):pi; %-- omega hat frequency vector

H = freqz(bb, 1, ww); %<--Matlab function for frequency response

subplot(3,1,1);

plot(ww, abs(H)), grid on

subplot(3,1,2);

plot(ww, angle(H)), grid on

xlabel('Normalized Radian Frequency')

subplot(3,1,3);

zplane(bb,1); % pole-zero plots

Save (click on edit then save as) and run the script file and get a screenshot of the script file and plots. Record the real and imaginary parts of the zeros and poles. Are the zeros on the unit circle? Where is the band pass now centered (check this both in the frequency response plot and the pole-zero plot)?

Now in the above script file, change the center frequency of the bandpass to w0 = 2*pi/3 and repeat the above steps. How does the new plots compare to the previous one? Are the plots as expected?

Design of a high-pass filter from a First Difference Filter:

Here you will design a simple high pass filter by using the first difference filter.

Repeat by changing the filter coefficients to that of a first difference filter (filter coefficients bb = [1 -1].

Create a new script file called Lab3_Part3_3

bb = [1 -1]; %-- Filter Coefficients

ww = -pi:(pi/500):pi; %-- omega hat frequency vector

H = freqz(bb, 1, ww); %<--Matlab function for frequency response

subplot(2,1,1);

plot(ww, abs(H)), grid on

subplot(2,1,2);

plot(ww, angle(H)), grid on

xlabel('Normalized Radian Frequency')

zplane(bb,1); % pole-zero plot

How many zeros are there? Which frequency is blocked? What type of filter is this (lowpass, bandpass, high pass…)? Read and record the real and imaginary parts of the zeros and poles.

Part 4: Design of FIR filters using the windowing method:

The concept of windowing is widely used in signal processing. The basic idea is to extract a finite section of a very long signal x[n] via multiplication w[n].x[n]. This approach works if the window function w[n] is zero outside of an interval. For example, consider the simplest window function which is the L-point rectangular window defined as

w[n] =1 for 0< n < L-1

w[n] 0 elsewhere

The impulse response for an ideal lowpass filter can be shown to be the sinc function

is the cutoff frequency of the ideal LPF, which separates the passband

from the stop band ( defined in Fig. 1 below).

One simple approach to designing a practical FIR filter is to truncate the impulse response of an ideal filter (which is infinite in length). This can be accomplished with a window function, so we usually say that the practical FIR filter has an impulse response that is a windowed version of the ideal impulse response. For example, we could make a length-23 FIR lowpass filter by taking the center portion of the sinc function:

where wr[n] is a 23-point rectangular window.

Procedure:

Create a new script file called Lab3_part4

Write the following commands:

% Filter design with the Window method

M = 22; %-- even Filter order

nn = 0:M; %--'n' index vector

nsh = nn - M/2; %--- Put main lobe in center of h[n]

wb = 0.32*pi; %-- 'Bandwidth parameter' of sinc

bb = sin(wb*nsh)./(pi*nsh); %-- Filter Coefficients

bb(find(nsh==0)) = wb/pi; %-- fix the divide by zero

ww = -pi:(pi/500):pi; %-- omega hat frequency vector

H = freqz(bb, 1, ww); %<--Matlab function for frequency response

subplot(3,2,1);

plot(ww, abs(H)), grid on % Plot magnitude

subplot(3,2,2);

plot(ww, angle(H)), grid on % Plot phase

xlabel('Normalized Radian Frequency')

subplot (3,2,3)

zplane(bb,1); % Pole-zero plot

Save and run the file and get a screenshot of the file and plots. Is this a lowpass, highpass or bandpass filter? Are the results as expected. Read and record the real (x) and imaginary (y) parts values of the zeros and poles by hovering with the mouse over their symbols (x: poles, o: zeros). How close are the zeros to the unit circle where frequencies are strongly attenuated?

Repeat the above by changing the filter order to M = 10, 30, 50, 80.

Each time save and run the file and get a screenshot of the file and the plots.

Is this a low pass, high pass or bandpass filter? How do the plots change and compare? How does this filter compare to the lowpass filter obtained a previous section (based on L-point running sum filter) What do you conclude?

Part 5: Design of a Nulling filter to block undesired tones:

In many applications, we wish to block some undesirable frequencies, such as the parasitic 60 Hz frequency originating from the mains. In the following, you will design a nulling filter that blocks 3 frequencies of tones that were added to a voice recording. You will first upload the recoded voice “SunshineSquare. wav”. To do this you go to Blackboard in the labs folder, and download this file to your desktop. Then go to Matlab click on Home, and upload this file. This file will be added to the list on the left “current folder” window.

As discussed in class and your textbook (Chapt. 7) a nulling filter can be considered a cascade of first order filters whose system function is of the form:

Hk(z) = (1 - 2cos(ωk)z-1 + z-2)

where the kth blocked normalized frequency is ωk . To block dc the following system function is required:

H0(z) = (1 - z-1)

Here 3 tones are to be blocked, in addition to dc. The system function of the cascade is

H(z) = H0(z)*H1(z)*H2(z)*H3(z)

= (1 - z-1)(1 - 2cos(ω1)z-1 + z-2) (1 - 2cos(ω2)z-1 + z-2) (1 - 2cos(ω3)z-1 + z-2)

The following script file will plot the spectrogram of the audio file. From this spectrogram you will be able to identify the undesirable frequencies (correspond to the yellow vertical bands).

To get the filter coefficients bk of the cascade you can multiply out the above factor form of the system function and then take the inverse z- transform. The bk coefficients will be the coefficients of the overall system function H(z). Here you will use the convolution function conv( ) to find the unit impulse response directly in the time n-domain as h = convolution of h0, h1, h2, and h3. The output yy[n] is then computed as the convolution of h[n] and the input xx[n].

This is output as an analog sound where you will hear how this filter has blocked the undesirable frequencies. You will also calculate the frequency response and plot its magnitude and phase vs. You will then output the unfiltered audio file and compare.

Create a new script file named Lab3_tone_removal.m.Write the following commands:

% This will block three tones from an audio recording,

% where the tones were added to a voice recording.

%The frequencies of the tones are identified from the spectrogram

% then a filter developed to block them as well as dc.

% The output is heard and compared to the input. The frequency response is plotted.

[xx, fs] = audioread('SunshineSquare .wav’); % input audio file

xx = xx’; % transforms the row vector into a column vector

subplot (4,1,1)

spectrogram(xx, [ ],fs) % spectrogram and sampling frequency fs. toes to

% be removed are the higher amplitude (yellow vertical stripes) frequencies)

subplot (4,1,2)

bb0 = [1 -1]; % Filter coefficients to block dc

bb1 = [1 -2*cos(0.285*pi) 1]; % Filter coefficients to block the first one

bb2 = [1 -2*cos(0.57*pi) 1]; % Filter coefficents to block the 2nd tone

bb3 = [1 -2*cos(0.86*pi) 1]; % Filter coefficients to block the 3d tone

bb =conv(bb3,conv(bb2,conv(bb1,bb0))); % convolution of the impulse response of

% the cascade of the 4 blocking filters

yy = conv(bb,xx); % filtered output

spectrogram(yy,[],fs) % spectrogram of the output

soundsc(yy,fs) % analog output as sound

subplot (4,1,3)

ww = -pi:pi/500:pi; % Normalized frequency array

H = freqz(bb,1,ww); % frequency response of the nulling filter

plot(ww, abs(H)), grid on % plot magnitude of frequency response

subplot (4,1,4)

plot (ww,angle(H)), grid on % plot phase of frequency response

xlabel('Normalized Radian Frequency')

Save this file then run it (click on editor, save, run).

Listen to the filtered recorded voice.

Now to play the original unfiltered audio file in the Matlab Command window (bottom) enter:

>> soundsc(xx, fs)

How do the filtered and unfiltered audio files compare? Are the undesirable tones noticeably suppressed? Note that if you still hear some of the tones this will be because the actual frequencies read from the spectrogram are approximate and need to be more fine tuned.

Look at the spectrograms of the input and output audio files? are the undesirable tones (yellow vertical stripes) removed?

Now look at the frequency response magnitude plot. At what normalized frequencies is it = 0? Compare these frequencies to the undesirable tine frequencies that you intended to suppress. Are they the same as indicated on the spectrogram of the unfiltered audio signal? What do you conclude?

Take a screenshot of the script file and plots.

Part 6: Design of IIR filters

An IIR Filter of order N is defined by the recursive difference equation:

y[n] = a1 y[n-1] + a2 y[n-2] + ….. + aN y[n-N] + b0x[n] b1x[n-1] +…. + bM x[n-M]

and alternatively by the system function

H(z) = ( b0 + b1 z-1 ….. + bM z-M ) / (1 + a1 z-1 + a2 z-2 …… + aN z-N )

The filter is specified by 2 arrays of coefficients:

aa = [1, a1,…, aN]

bb = [b0, b1,…, bM]

In the z-domain it is defined by its system function H(z). The zeros of H are the values of z for which H(z) = 0. The poles are the values for which H(z) is infinite. The frequency response is given by the Matlab function freqz(bb,aa). The plot of poles and zeros on the complex or z plane, including the unit circle (on which the system function becomes the frequency response) is given by zplane(bb,aa). To get the real part (x) and imaginary part (y) of the poles and zeros, hover with the mouse over their symbols ( x: poles, o: zeros).

IIR filters are most common as analog filters. They provide more versatility than FIR filters in that one can move both the poles and zeros at wish close to the desired region near the unit circle in order to better shape the frequency response. With FIR filters all poles are fixed at the origin.

However, with IIR filters, one has to make always sure the poles are within the unit circle, for otherwise it will be unstable and therefore unrealizable which is not an issue with FIR filters.

Create a script file called Lab3_part6.

Write the following commands:

% Computes the frequency response for an IIR filter given its coefficients.

% Displays magnitudes and phase plots and pole-zero plot of the corresponding

% system function

bb = [64 300 50]; %-- Numerator coefficients of system function H(z)

aa = [1 0.68 0.64] % —Denominator coefficients of system function H(z)

ww = -pi:(pi/500):pi; %-- omega hat radian frequency vector

H = freqz(bb, aa, ww); %— frequency response

subplot(3,1,1);

plot(ww, abs(H)), grid on % Plot magnitude of H

ylabel('Magnitude')

subplot(3,1,2);

plot(ww, angle(H)), grid on % Plot phase of H

ylabel('Phase')

xlabel('Normalized Radian Frequency')

subplot(3,1,3)

zplane(bb,aa)% computes and plots the poles and zeros of system function H(z)

Save and run the file and get a screenshot. What type of filter is this (lowpass, bandpass, highpass, etc.)? In what range of frequencies is the output strongly attenuated? Which frequency ranges are passed (near peaks). Look at the pole-zero plot. Is there a zero or pole near these frequencies. Hover with the cursor of the mouse over the poles and zeros and record their real (x) and imaginary (y) parts. Note that frequencies (normalized) are the angles on the unit circle. Repeat for regions where there is a peak in the frequency response and output.

Now repeat the above steps after changing the IIR filter coefficients in the script file as follows:

bb = [64 180 100];

aa = [1 0.68 4.64];

Repeat again for the following filter coefficients:

bb = [64 300 50];

aa = [1 0.68 0.64];

Practice as above by using different filter coefficients aa and bb. You do not need to submit these

latter results.

IV- Conclusion:

Write a conclusion summarizing your results. Were your objectives met? Where your results as expected? If not why? What are your thoughts about this lab project and suggestions if any to possibly improve it.

Image files to upload to Matlab

4