% m-file file name: kailpf2.m
% Function for designing FIR lowpass filter using Kaiser window.
% Specifications are given in terms of ap, aa.
% Written by Dr. James S. Kang, Professor
% Department of Electrical and Computer Engineering
% California State Polytechnic University, Pomona
% ap = passband ripple in dB.  ap = 20*log10((1+d1)/(1-d1)).
% aa = stopband attenuation in dB.  aa = -20*log10(d2).
% fplpf = passband cutoff frequency in Hz.
% fslpf = stopband cutoff frequency in Hz.  fs > fp.
% fs = sampling rate in Hz.
% References:
% Leland B. Jackson, Digital Filters and Signal Processing, 3rd ed., Kluwer Academic Publishers, 1996.
% T. W. Parks and C. S. Burrus, Digital Filter Design, John Wiley and Sons, Inc., 1987.
% Andreas Antoniou, Digital Filters: Analysis, Design, and Applications, 2nd ed., McGraw-Hill, 1993.
% J. G. Proakis and D. G. Manolakis, Digital Signal Processing, 3rd ed., Prentice-Hall, 1996.
% A. V. Oppenheim and R. W. Shafer, Discrete-Time Signal Processing, Prentice-Hall, 1989.
% Try >>kailpf2(2,55,2000,2500,8000)
function[h,H] = kailpf2(ap,aa,fplpf,fslpf,fs)
Ts = 1/fs;  % Ts = sampling interval.
fclpf = (fplpf+fslpf)/2;
qplpf = 2*pi*fplpf*Ts; % qplpf = digital passband cutoff frequency.
qslpf = 2*pi*fslpf*Ts; % qslpf = digital stopband cutoff frequency.
qclpf = 2*pi*fclpf*Ts; % qclpf = digital cutoff frequency.
d1 = (10^(0.05*ap)-1)/(10^(0.05*ap)+1);  % d1 = fractional ripple above and below 1 in the passband.
d2 = 10^(-0.05*aa);  % d2 = ripple in the stopband.
if d1 <= d2 dm = d1; else dm = d2; end 
aa = -20*log10(dm); % aa = stopband attenuation in dB.
if aa <= 21 b = 0;
   elseif aa > 21 & aa <= 50 b = 0.5842*(aa-21)^0.4 + 0.07886*(aa-21);
      else
         b = 0.1102*(aa-8.7);
      end
end
if aa <= 21 D = 0.9222; else D = (aa-7.95)/14.36; end
N = ceil(fs*D/(fslpf-fplpf)+1);
%if rem(N,2) == 0 N = N+1; end;
disp('Delta 1 = ')
disp(d1)
disp('Delta 2 = ')
disp(d2)
disp('Delta = ')
disp(dm)
disp('Stopband attenuation in dB = aa = ')
disp(aa)
disp('Cutoff frequency in Hz = ')
disp(fclpf)
disp('Cutoff frequency in radians = ')
disp(qclpf)
bt = fslpf - fplpf;
disp('Width of transition band in Hz = bt = ')
disp(bt)
disp('beta = ')
disp(b)
disp('D = ')
disp(D)
disp('The length (order) of the filter is')
disp(N)
%
% Calculating the desired shifted response.
%
for k = 1:N
   if (k-1) ~= (N-1)/2
      hds(k) = (qclpf/pi)*sin((qclpf*((k-1)-(N-1)/2)))/(qclpf*((k-1)-(N-1)/2));
   else
      hds(k) = qclpf/pi;
   end
end
for k = 1:N
   t(k)=k-1;
end
figure
stem(t,hds);grid;
title('Desired Shifted Response')
xlabel('Time')
ylabel('Amplitude')
disp('The desired shifted response is')
disp(hds)
%
% Calculating and plotting the window function.
%
for k = 1:N
   w(k) = bessel(0,j*b*sqrt(1-(1-2*(k-1)/(N-1))^2))/bessel(0,j*b);
end
for k = 1:N
   t(k) = k-1; 
   z(k) = 0; 
   if w(k) >= 0 Hi(k) = w(k); Lo(k) = 0; 
      else Hi(k) = 0; Lo(k) = -w(k); 
   end
end
figure
stem(t,w)
grid
%errorbar(t,z,Lo,Hi,'w')
title('Window Function')
xlabel('Time')
ylabel('Amplitude')
disp('The window function is')
disp(w)
%
% Calculating and plotting the impulse response of the FIR filter.
%
for k = 1:N
   h(k) = hds(k)*w(k);
end
for k = 1:N
   t(k) = k-1; 
   z(k) = 0; 
   if h(k) >= 0 Hi(k) = h(k); Lo(k) = 0; 
      else Hi(k) = 0; Lo(k) = -h(k); 
   end
end
figure
stem(t,h)
grid on
%errorbar(t,z,Lo,Hi,'h')
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')
disp('The impulse response is')
disp(h)
%
% Calculating and plotting the frequency response.
%
M = 400;
for n = 1:M+1
   f(n) = (n-1)/M;
   q(n) = 2*pi*f(n);
   if rem(N,2) == 1
        H(n) = h((N-1)/2+1)+2*h(1)*cos(q(n)*(-(N-1)/2));
   else
        H(n) = 2*h(1)*cos(q(n)*(-(N-1)/2));
   end
   if rem(N,2) == 1 endk = (N-1)/2; else endk = N/2; end
   for k = 2:endk
      H(n) = H(n) + 2*h(k)*cos(q(n)*(k-1-(N-1)/2));
   end
   H(n) = exp(-j*q(n)*(N-1)/2)*H(n);
   mag(n) = abs(H(n));
   dB(n) = 20*log10(abs(H(n)));
   if dB(n) < -200 dB(n) = -200; end
   pha(n) = angle(H(n))/pi; 
end
figure
plot(f,mag)
grid on
title('Magnitude Response')
xlabel('Frequency')
ylabel('Magnitude')
figure
plot(f,dB)
grid on
title('Magnitude Response in dB')
xlabel('Frequency')
ylabel('Magnitude in dB')
figure
plot(f,pha)
grid on
title('Phase Response')
xlabel('Frequency')
ylabel('Phase in Radians/pi')
%
% Calculating and plotting poles and zeros.
%
nozeros = 0;
k = 1;
while abs(h(k)) < 10^(-10) & k == nozeros + 1 & k < (N-1)/2
   nozeros = nozeros + 1;
   k = k+1;
end
stop = N - 2*nozeros;
for k = 1:stop
  h1(k) = h(k+nozeros);
end
dz=ROOTS(h1);
for k = 1:stop-1
  dp(k) = 0;
end
for m = 1:201
   x(m) = (m-101)/100;
   ty(m) = sqrt(1-x(m).*x(m));
   by(m) = - ty(m);
end
figure
plot(x,ty)
axis equal;
grid on
hold on
plot(x,by) 
plot(real(dp),imag(dp),'x')
plot(real(dz),imag(dz),'o')
title('Digital Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Digital Pole Locations')
disp(dp)
disp('Digital Zero Locations')
disp(dz)