%truncated FIR low pass filter demonstration
%
%
clear
clf
%set alpha = 4 (frequency cutoff equal to r/4,
%where r is the sampling rate).
alpha = 4.;
%
%  Use 128 point filters.
%
N = 128;
%assign normalized frequency range, f = k/N
%where f=1 is sampling rate
df=1/1024;
f=0:df:1/2;

%
%  First, use a rectangular window.
%
figure(1)

for ii=1:N
j=ii-N/2;
if j==0
y(ii)=2/alpha;
else
y(ii) = (2/alpha)*sin(pi*2*j/alpha)/(pi*2*j/alpha);
end
end

for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end

spec=abs(s);
%
%  Plot the windowed filter.
%
subplot(2,1,1)
plot(1:N,y)
title('Rectangular Window; N=128')
axis([0 N-1 -2/alpha 2/alpha]);
%
% Plot the frequency response.
%
subplot(2,1,2)
semilogy(f,spec);
xlabel('f/r')
ylabel('|W_k|');
axis([0 1/2 1.0e-5 10]);
%
% Plot the phase response.
% 
%phase=unwrap(angle(s),3)*180/pi;
%subplot(3,1,3);
%plot(f,phase);
%xlabel('f/r');
%ylabel('\theta (degrees)');

disp('Aster and Borchers Figure 5.4: FIR weights and response for a 128-point rectangular window FIR realization of a low pass filter with a desired cutoff frequency of f=r/4.')
disp('press enter to continue')
pause;

%
%  Now, the Bartlett window.
%

ystore = y;

y=y.*(bartlett(length(y)))';

for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end

spec=abs(s);

figure(2)
subplot(2,1,1)
plot(1:N,y)
title('Bartlett Window; N=128')
axis([0 N-1 -2/alpha 2/alpha]);
subplot(2,1,2)
semilogy(f,spec);
xlabel('f/r')
ylabel('|W_k|');
axis([0 1/2 1.0e-5 10]);
%
% Plot the phase response.
% 
%phase=unwrap(angle(s),3)*180/pi;
%subplot(3,1,3);
%plot(f,phase);
%xlabel('f/r');
%ylabel('\theta (degrees)');
%
disp('Aster and Borchers Figure 5.5: FIR weights and response for a 128-point Bartlett window FIR realization of a low pass filter with a desired cutoff frequency of f=r/4.')
disp('Note the reduction in ripple near the transition band relative to the simple truncation from Figure 5.4')
disp('The Bartlett window is the same as a triangle except there are leading and trailing 0s.')
disp('press enter to continue')
pause;

y=y.*(hamming(length(y)))';

for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end

spec=abs(s);

%
% Now, the Hamming window.
%
figure(3)
subplot(2,1,1)
plot(1:N,y)
title('Hamming Window; N=128')
axis([0 N-1 -2/alpha 2/alpha]);
subplot(2,1,2)
semilogy(f,spec);
xlabel('f/r');
ylabel('|W_k|');
axis([0 1/2 1.0e-5 10]);
%
% Plot the phase response.
% 
%phase=unwrap(angle(s),3)*180/pi;
%subplot(3,1,3);
%plot(f,phase);
%xlabel('f/r');
%ylabel('\theta (degrees)');

disp('Aster and Borchers Figure 5.6: FIR weights and response for a 128-point Hamming window FIR realization of a low pass filter with a desired cutoff frequency of f=r/4. ')
disp('Note the reduction in ripple near the transition band relative to the simple truncation from Figure 5.4 and the Bartlett window from Figure 5.5.')
disp('Recall the Hamming window is a raised cosine with leading and trailing 0s')
disp('The tradeoff for smoother response and better rejection outside of the desired passband is to have a more gradual transition.')
