Getting Started with Filtering
To teach filtering features of MATLAB.
Students will be able to carry out filtering operations within the MATLAB interface and understand the basic functions of convolution and filter.
About two hours
Students should read the Digital Filter notes provided on the class website. They should have a good understanding of the concepts of difference equations and their representation through a delay operator. Also understand the z-plane and the meaning of transfer function of a filter (phase and amplitude response).
Before you start, enable the diary function. When you are done, save it as [yourname]A2LabDiary and edit out all unnecessary steps.
A. FIR Filtering
1) FIR Filtering by Buffering
Filtering operation can be described as a sliding dot product (vector dot product is done in MATLAB by *) between the filter coefficients, or coefficients of the tap delays, and a corresponding segment of the audio vector. Since the audio vector is longer then the number of coefficients, an appropriate segment has to be extracted from the sound vector, multiplied by the filter coefficients and then the result is written at the output.
In this exercise we will create a moving average of three samples. To do so we define a filter whose coefficients are a0 = 1, a1 = 1/2, a2 = 1/3. Note that the older samples are averaged with amplitude ½ for first and 1/3 for two delayed samples.
Create a filter vector a = [1 1/2 1/3]
Create a row of 100 ones followed by three zeros: b=[ones(1,100) 0 0 0]
Use the buffer function (type "help buffer" into Matlab to view the help file) to create a matrix from vector y of length 3 and overlap 2: y=buffer(b,3,2)
Use vector dot product to create the output signal: x=a*y
Notice the first two and last three columns. Why are they like this?
2) FIR Filtering by Convolution
Convolution is the operation equivalent of a sliding dot product. It is performed by MATLAB conv command.
Create a row of ten ones: a=ones(1,10)
Create another row of ten ones: b=ones(1,10)
Convolve the two: conv(a,b)
Explain this result (hint: try drawing a picture of the intermediate steps)
3) Impulse response
Create an input signal called imp that comprises of a single sample of amplitude 1 followed by ten zeros: x=[1 zeros(1,10)]
Create B filter coefficients: b=[0.8 0.6 0.7 0.5 0.3 0 0.9]
Why do these coefficients make this filter FIR?
Perform a convolution of the filter coefficients a with the impulse: y=conv(b,x)
Plot the result and describe what you see: plot(y)
How long is the output signal? How does it relate to the B coefficients? How would this be different if x were [1 -1 zeros(1,10)]?
B. IIR Filtering
1) Filter function.
Learn about the filter function:
FILTER One-dimensional digital filter.
Y = FILTER(B,A,X) filters the data in vector X with the
filter described by vectors A and B to create the filtered
data Y. The filter is a "Direct Form II Transposed"
implementation of the standard difference equation:
a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na)
Perform a low pass filter that executes the following difference equation:
y(n) = y(n-1) + 0.1[x(n) - y(n-1)]
a) write down the numerator and denominator of the filter (A and B vectors)
b) create 1000 samples of white noise by typing y=rand(1,1000)*2-1
c) execute the following commands
h = spectrum.welch; % Create a Welch spectral estimator.
psd(h,y,'Fs',Fs); % Calculate and plot the PSD.
Save the graphs and describe what you see.
d) use the filter command on an impulse: create an input signal imp that comprises of a single sample of amplitude1 followed by many zeros (several hundreds). Use the filter command to perform filtering of imp. Plot the result and save the graph as [yourname]A2FiltImp.jpg
C. Frequency response
Plotting z-plane in MATLAB
Copy and paste the following code into Matlab:
% Set up vector for zeros
z = [j ; -j];
% Set up vector for poles
p = [-1 ; .5+.5j ; .5-.5j];
title('Pole/Zero Plot for Complex Pole/Zero Plot Example');
Why do only some of the elements of p contain complex numbers? (hint: it's a trick question)