FIR Filter Design and Testing Using the DFT

% 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');