Note: it is strongly recommended that you print this document in Firefox, not IE. Thanks!
Filter method one simply uses the convolution of the filter and the input to obtain the output. This could not be used for an IIR filter becasue convolution does not look at past values of the output.
Filter method two uses the filter function, which constructs a discrete time filter using the numerator and denominator polynomials and applies it to the input data sequence. It can be used for an IIR filter.
Filter method three generates a discrete time transfer function using the Control System Toolbox. The transfer function is then applied to the input. This method can also be used for an IIR filter.
Filter method four is very similar to filter method three, except that the transfer function is in terms of z-1 instead of z. This format is more commonly used in the field of DSP.
Shown below are the graphs of each frequency sinusoid with its 100-point running average. Each graph exhibited irregular behavior for the first 100 datapoints before settling on a steady-state response. For example, in Figure 1a, the first 100 points are a ramp from 0 to 1, and proceeds to stay at 1 for the remaining points. This is because the function is zero-padded for negative n values on the input; the filter sees zeros everywhere except where the input is explicitly defined. At point 50, the output in the first graph is 0.5, since the average of 50 * 0 + 50 * 1 datapoints is 0.5.
Other graphs show more interesting behavior in the first 100 datapoints. Figure 1b shows 1/2 cycle of a sine in the first 100 datapoints, demarcated by the sharp change in the output slope at point 100. Figure 1c shows a full cycle, 1d shows 1 1/2 cycles, etc. This is a result of the fact that 100 datapoints is an exact multiple of half the period in each case: 40Hz -> 200 datapoints per cycle, 80Hz -> 100 datapoints, 120Hz -> 50 datapoints, etc.
 Figure 1a: 0Hz Sine and 100-point Average |
 Figure 1b: 40Hz Sine and 100-point Average |
 Figure 1c: 80Hz Sine and 100-point Average |
 Figure 1d: 120Hz Sine and 100-point Average |
 Figure 1e: 160Hz Sine and 100-point Average |
 Figure 1f: 200Hz Sine and 100-point Average |
 Figure 1g: 240Hz Sine and 100-point Average |
 Figure 1h: 280Hz Sine and 100-point Average |
Below is the code for plotting the eight sampled sines and their running averages.
Fs=8000;
b=ones(1,100)./100;
n=0:399;
t=n/Fs;
x=zeros(8,400);
y=zeros(8,499);
clf;
for i=1:8
x(i,:)=cos(2*pi*(i-1)*40*t);
y(i,:)=conv(x(i,:),b);
figure();
hold on;
plot(t, x(i,:), '--r');
plot(t, y(i,1:400), 'b');
legend(sprintf('%dHz input',(i-1)*40),'100 point average');
xlabel('Time (s)');
ylabel('Amplitude');
hold off;
end
Figure 2 below shows a graph of steady-state magnitude against frequency. Note that at each multiple of 80Hz, the magnitude is zero. This is another effect of the fact that 100 datapoints, the length of the average window, is an exact multiple of the entire period of the sinusoid in each case. The average of the cosine function over an entire period, or some integer multiple of a period, is zero.
 Figure 2: Magnitudes of Steady State Averages Vs. Frequency |
Since each point of the 100-element filter must be multiplied with the specific datapoint lining up with it, and all of these products must be added together, a single output value requires 100 multiplys and 99 additions.
Here is the code for the frequency graph.
Fs=8000;
b=ones(1,100)./100;
n=0:399;
t=n/Fs;
x=zeros(8,400);
y=zeros(8,499);
clf;
for i=1:8
x(i,:)=cos(2*pi*(i-1)*40*t);
y(i,:)=conv(x(i,:),b);
if i==1
mags(i)=max(y(i,1:400));
else
mags(i)=max(y(i,400-Fs/((i-1)*40):400));
end
end
stem(0:40:280,mags);
axis([0,280,0,1.2]);
set(gca,'xtick',0:40:280);
xlabel('Frequency (Hz)');
ylabel('Magnitude');
We observed that the running average could be calculated with the following equation:
.
In words, take the last value of the running average, add the weighted influence of the current input value, and subtract off the weighted influence of the input value 100 datapoints ago. This technique requires only two additions and two multiplications per data point, but requires us to have the previous output value to calculate the next. From the above equation, we found a transfer function:
.
We entered our transfer function into MATLAB using Method 4 from Task 1. This work is shown in the code below.
Fs=8000;
b=ones(1,100)./100;
n=0:399;
t=n/Fs;
x=zeros(8,400);
y=zeros(8,400);
clf;
numerator = zeros(1,101);
numerator(1) = .01;
numerator(101) = -.01;
denominator = [1, -1];
h = tf(numerator,denominator,1/Fs, 'variable', 'z^-1')
for i=1:8
x(i,:)=cos(2*pi*(i-1)*40*t);
y(i,:)=lsim(h,x(i,:))';
y_fir(i,:)=conv(x(i,:),b);
end
[Y,T]=impulse(h);
for i=1:length(Y)
if Y(i) < .00001
break;
end
end
i-1
clf;
subplot(2,1,1);
hold on;
plot(t,y(2,1:400));
plot(t,y_fir(2,1:400),':g');
legend('IIR Filter','FIR Filter');
xlabel('Time (s)');
ylabel('Amplitude');
title('Comparison of Techniques');
subplot(2,1,2);
plot(t,y(2,1:400)-y_fir(2,1:400));
xlabel('Time (s)');
ylabel('Magnitude');
title('Difference Between Techniques');
Figure 3 below compares the outputs of the two techniques. The top subplot shows the two outputs plotted on the same axes, while the bottom subplot shows the difference between the outputs of the two techniques. Note that the magnitudes of difference are on the order of 10-16, which is comparable to rounding error for floating point numbers. This demonstrates the correctness of our code. Additionally, lines 21-27 include code to measure the number of datapoints for which the impulse response is nonzero. Since 100 datapoints are averaged, this number should be 100, which was confirmed by execution of the code.
 Figure 3: Comparison of FIR and IIR Techniques |