% m-file file name: kaibpf2.m
% Function for designing FIR bandpass 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).
% f1bpf = lower passband cutoff frequency in Hz.
% f2bpf = upper passband cutoff frequency in Hz. 
% f3bpf = lower stopband cutoff frequency in Hz.
% f4bpf = upper stopband cutoff frequency in Hz. 
% f3bpf < f1bpf < f2bpf < f4bpf < fs/2.
% 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 >>kaibpf2(2,55,2000,2500,1600,3000,8000)
function[h,H] = kaibpf2(ap,aa,f1bpf,f2bpf,f3bpf,f4bpf,fs)
Ts = 1/fs;  % Ts = sampling interval.
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 f1bpf-f3bpf <= f4bpf-f2bpf bt = f1bpf-f3bpf; else bt = f4bpf-f2bpf; end
f1bpf = f1bpf-bt/2; % f1bpf = lower passband cutoff frequency.
f2bpf = f2bpf+bt/2; % f1bpf = upper passband cutoff frequency.
q1bpf = 2*pi*f1bpf*Ts; % q1bpf = digital lower passband cutoff frequency.
q2bpf = 2*pi*f2bpf*Ts; % q2bpf = digital upper passband cutoff frequency.
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/bt+1);
disp('Delta 1 = ')
disp(d1)
disp('Delta 2 = ')
disp(d2)
disp('Delta = ')
disp(dm)
disp('Stopband attenuation in dB = aa = ')
disp(aa)
disp('Lower cutoff frequency in Hz = ')
disp(f1bpf)
disp('Upper cutoff frequency in Hz = ')
disp(f2bpf)
disp('Lower cutoff frequency in radians = ')
disp(q1bpf)
disp('Upper cutoff frequency in radians = ')
disp(q2bpf)
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.
%
%if rem(N,2) == 0 N = N+1; end;
for k = 1:N
  if (k-1) ~= (N-1)/2
   hds(k) = sin(q2bpf*(k-1-(N-1)/2))/(pi*(k-1-(N-1)/2))-sin(q1bpf*(k-1-(N-1)/2))/(pi*(k-1-(N-1)/2));
  else
   hds(k) = q2bpf/pi - q1bpf/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;
%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;
title('Magnitude Response')
xlabel('Frequency')
ylabel('Magnitude')
figure
plot(f,dB);grid;
title('Magnitude Response in dB')
xlabel('Frequency')
ylabel('Magnitude in dB')
figure
plot(f,pha);grid;
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)