% m-file file name: papbsf.m
% FUnction for designing FIR bandstop filter using Papoulis window.
% Written by Dr. James S. Kang, Professor
% Department of Electrical and Computer Engineering
% California State Polytechnic University, Pomona
% N = length (order) of the filter (impulse response).
% f1bsf = lower cutoff frequency.
% f2bsf = upper cutoff frequency.
% 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.
function[h,H] = papbsf(N,f1bsf,f2bsf,fs)
Ts = 1/fs;  % Ts = sampling interval.

q1bsf = 2*pi*f1bsf*Ts; % q1bsf = digital lower passband cutoff frequency.
q2bsf = 2*pi*f2bsf*Ts; % q2bsf = digital upper passband cutoff frequency.
%
% Calculating the desired shifted response.
%
if rem(N,2) == 0 N = N+1; end;
disp('The length (order) of the filter is')
disp(N)
for k = 1:N
   hds(k) = -(q2bsf/pi)*sinc(q2bsf*(k-1-(N-1)/2)/pi)+(q1bsf/pi)*sinc(q1bsf*(k-1-(N-1)/2)/pi);
   if k == (N-1)/2 + 1 hds(k) = hds(k) + 1; end
end
disp('The desired shifted response is')
disp(hds)
%
% Calculating and plotting the window function.
%
for k = 1:N
   w(k) = (1/pi)*abs(sin(2*pi*(k-1-(N-1)/2)/(N-1)))+(1-2*abs(k-1-(N-1)/2)/(N-1))*cos(2*pi*(k-1-(N-1)/2)/(N-1));
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);grid;
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
%
%Normalization
%
sum1=sum(h);
for k = 1:N
   h(k) = h(k)/sum1;
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); grid
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')
%
% Calculating and plotting the frequency response.
%
M = 400;
for n = 1:M+1
   f(n) = (n-1)/M;
   q(n) = 2*pi*f(n);
   H(n) = h((N-1)/2+1)+2*h(1)*cos(q(n)*(-(N-1)/2));
   for k = 2:(N-1)/2
      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)); 
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')
%
% Calculating and plotting poles and zeros.
%
for j = 1:N-2
   h1(j) = h(j+1);
end
dz = ROOTS(h1);
maxz = 0;
for j = 1:N-3
   if abs(real(dz(j))) > maxz & abs(real(dz(j))) < 10^2 maxz = abs(real(dz(j))); end
end
for j = 1:N-3
   if abs(dz(j)) > 10^15 dz(j) = 1.2*maxz; end
end
dz(N-2) = 1.2*maxz;
for k = 1:N-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,'-');grid;
axis equal;
hold on
plot(x,by,'-')
plot(real(dp),imag(dp),'x')
for i = 1:N-2
   if dz(i) == 1.2*maxz 
      plot(real(dz(i)),imag(dz(i)),'o')
   else
      plot(real(dz(i)),imag(dz(i)),'o') 
   end     
end
title('Digital Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Digital Pole Locations')
disp(dp)
disp('Digital Zero Locations')
for i = 1:N-2
   if dz(i) == 1.2*maxz
      disp('zero is at infinity')
   else
      disp(dz(i))
   end
end