%
% In this demo, we'll use the strategy of inverse Fourier transforming the
% desired frequency response, and then truncating and windowing the resulting
% infinite sequence to get a FIR filter.  
%
%
% First, read in a 20 second audio clip.  (Bruce Cockburn)
%
[Y,FS,NBITS]=wavread('clip.wav');
%
% First, we'll play with a simple 5 point running average filter.
%
disp('Here are the filter coefficients for a 17 point running average filter');
M=17;
w=ones(M,1)/M
%
% Plot the spectrum of the filter.
%
figure(1);
df=1/1024;
f=-0.5:df:0.5;
s=zeros(size(f));
for k=1:length(f)
  for l=1:M,
    s(k) = s(k) + (w(l)*exp(sqrt(-1)*2*pi*(l-1)*f(k)));
  end;
end;
plot(f,20*log10(abs(s)));
xlabel('Relative frequency (f/fs)');
ylabel('frequency response (db)');
print -deps runningaverage.ps
%
%
%
YFILTERED=conv(Y,w);
%
% The convolution will increase the length of the signal by M-1 samples. This
% happens because the filter has a response that lasts beyond the end of the
% the signal. We'll take these extra samples out.
%
YFILTERED=YFILTERED(1:882000);
%
% Save the filtered signal.
%
wavwrite(YFILTERED,FS,NBITS,'clipg.wav');
disp('Figure 1 shows the frequency response of the 17 point running mean filter.');
disp('Play clipg.wav');
disp('press enter to continue');
pause;
%
% Next, construct an M point filter sequence using (10) from the notes.  Feel
% free to adjust M to any odd number that you might want.  We'll start
% with M=15.
%
% In our case, we want a cutoff at a frequency of 2Khz, which is roughly
% 44100/22=2004.5 Hz. So, alpha=22.
%
% The only trick here is that MATLAB wants the sequence to be from 1 to M
% instead of -(M-1)/2 to (M-1)/2, so we shift the index of summation.
%
alpha=22;
M=15;
w=zeros(M,1);
for i=-(M-1)/2:(M-1)/2,
  w(i+(M-1)/2+1)=(2/alpha)*sinc(2*i/alpha);
end;
%
%
%
disp('Here are the truncated filter coefficients.');
w
%
% Plot the spectrum of the filter.
%
figure(2);
df=1/1024;
f=-0.5:df:0.5;
s=zeros(size(f));
for k=1:length(f)
  for l=1:M,
    s(k) = s(k) + (w(l)*exp(sqrt(-1)*2*pi*(l-1)*f(k)));
  end;
end;
plot(f,20*log10(abs(s)));
xlabel('Relative frequency (f/fs)');
ylabel('frequency response (db)');
print -deps unwindowedfreq.ps
%
% Now, apply the filter to our signal.
%
YFILTERED=conv(Y,w);
%
% The convolution will increase the length of the signal by M-1 samples. This
% happens because the filter has a response that lasts beyond the end of the
% the signal. We'll take these extra samples out.
%
YFILTERED=YFILTERED(1:882000);
%
% Save the filtered signal.
%
wavwrite(YFILTERED,FS,NBITS,'cliph.wav');
disp('Figure 2 shows the frequency response of the filter');
disp('To hear the clip, play cliph.wav');
disp('press enter to continue');
pause;
%
% There's noticable ripple in the frequency response, although it may not
% be audible.  
%
% Next, we apply a Bartlett window 
%
window=bartlett(M);
w2=w.*window;
df=1/1024;
f=-0.5:df:0.5;
s=zeros(size(f));
for k=1:length(f)
  for l=1:length(w2),
    s(k) = s(k) + (w2(l)*exp(sqrt(-1)*2*pi*(l-1)*f(k)));
  end;
end;
figure(3);
plot(f,20*log10(abs(s)));
xlabel('frequency, relative to sampling frequency');
ylabel('frequency response (db)');
print -deps windowedfreq.ps
%
% Now, apply the filter to our signal.  This time, we'll use the filter
% command.
%
YFILTERED=filter(w2,[1],Y);
%
% Save the filtered signal.
%
wavwrite(YFILTERED,FS,NBITS,'clipi.wav');
disp('Figure 3 shows the frequency response of the windowed filter');
disp('To hear the clip, play clipi.wav');
disp('press enter to continue');
pause;
%
% Finally, design a 2Khz low pass filter using the frequency sampling
% technique.
%
N=64;
fstep=FS/N;
freqs=(-22050:fstep:22050-fstep)';
w3t=ones(N,1);
w3t=w3t.*(abs(freqs)<2000);
w3=ifft(ifftshift(w3t));
%
% Plot the frequency response.
%
df=1/4096;
f=-0.5:df:0.5;
s=zeros(size(f));
for k=1:length(f)
  for l=1:length(w3),
    s(k) = s(k) + (w3(l)*exp(sqrt(-1)*2*pi*(l-1)*f(k)));
  end;
end;
figure(4);
plot(f,abs(s));
xlabel('frequency, relative to sampling frequency');
ylabel('|W(f)|');
print -deps freqsamp.ps
%
% Now, apply the filter to our signal.  This time, we'll use the filter
% command.
%
YFILTERED=filter(w3,[1],Y);
%
% Save the filtered signal.
%
wavwrite(YFILTERED,FS,NBITS,'clipj.wav');
%
% 
%
disp('Figure 4 shows the frequency response of the frequency sampling filter');
disp('This plot is not in db units, because the cutoff is down to 0 (where');
disp('db is not defined) at some frequencies.');
disp('To hear the clip, play clipj.wav');
disp('Notice that there is still a lot of high frequency content left.');
%
% All of the following clips try to low pass the original clip.
%
% clipg       17 point running mean filter
% cliph       Truncated inverse FFT filter, no window.
% clipi       Truncated inverse FFT filter, Bartlett window.
% clipj       64 point frequency sampling filter 
%
% The filters in clipg and cliph are fairly effective low pass filters at
% about 2Khz.  In clip i, the use of the Barlett window has smoothed
% out the ripple, but at the expense of broadening the frequency response
% to about 4Khz- this sounds much closer to the original clip than 
% clipg or cliph.  In clip j, there's lots of energy at some but not
% all higher frequencies, so this sounds much like the original clip, but
% "weird"
%