FIR Filter Design and Testing Using the DFT

From Class Wiki
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
% This program shows the design of an N element FIR low pass filter using the DFT.
% Rob Frohne
close all;
N=512;
fmax=4000;
T=1/(2*fmax);
H =[ones(1,N/32+1),zeros(1,N-N/16-1),ones(1,N/32)];% make the low pass filter function.
h = ifft(H); % Get the impulse response.
check = max(abs(imag(h)./real(h))) % Check to make sure is really is real.
figure(1)
f=-1/2/T:1/N/T:1/2/T-1/N/T;
Hshift=fftshift(H);
fc=-1/2/T:1:1/2/T;
hshift=fftshift(real(h));
hshiftk=kaiser(N,10)'.*hshift;  %If you want the kaiser window.
t=-N*T/2:T:(N/2-1)*T;
plot(t,real(15*hshift),'b-',t,15*hshiftk,'g-',t,kaiser(N,10),'c-') % The 15 makes the scales better.
title('The filter coefficients')
legend('h(t)','Kaiser window applied to h(t)','Kaiser Window');
figure(2);
Hc=zeros(size(fc));
Hck=zeros(size(fc));
for n=-N/2:1:N/2-1,
	Hc=Hc+hshift(n+N/2+1).*exp(-i*2*pi*fc*(n-N/2)*T);
	Hck=Hck+hshiftk(n+N/2+1).*exp(-i*2*pi*fc*(n-N/2)*T);
end
plot(f,Hshift,'r-',fc,abs(Hc),'b-',fc,abs(Hck),'g-');
legend('Desired Filter Response','Actual Response Without Kaiser', 'Actual Response with Kaiser');
title('Filter Response')
xlabel('Frequency (Hz)')
ylabel('Response as a Ratio')
figure(3);
plot(fc,20*log10(abs(Hc)),'b-',fc,20*log10(abs(Hck)),'g-');
title('Filter Response')
xlabel('Frequency (Hz)')
ylabel('Response in dB')
legend('Actual Response Without Kaiser (dB)', 'Actual Response with Kaiser(dB)');
% Try the filter out.
Totaltime = 5;
t = 0:T:Totaltime;
nt = sin(2*pi*(1/140/T.*t).*t);  %Swept signal generator.
Fs = 1/T;
system('espeak "Here is the unfiltered signal."') % These are using the linux espeak command.
soundsc(nt,Fs);
Y=filter(hshift,1,nt); % Filter without Kaiser window.
Yk=filter(hshiftk,1,nt); % Filter with Kaiser window.
system('espeak "Here is the filtered signal without windowed coefficients."')
soundsc(Y,Fs);
system('espeak "Here is the filtered signal with windowed coefficients."')
soundsc(Yk,Fs);
figure(4);
subplot(3,1,1), plot(nt,'y-');
title('Unfiltered Signal');
subplot(3,1,2), plot(Y,'b-');
title('Filtered Unwindowed Signal');
subplot(3,1,3), plot(Yk,'g-');
title('Filtered Windowed Signal');
xlabel('Sample number');