% m-file file name: dicbsf.m
% Function for designing digital inverse Chebyshev bandstop filter.
% Specifications are given in terms of amax, amin, f1bsf, f2bsf, f3bsf, f4bsf.
% Written by Dr. James S. Kang, Professor 
% Department of Electrical and Computer Engineering
% California State Polytechnic University, Pomona 
% amax = maximum loss in dB in the pass band.
% amin = minimum loss in dB in the stop band.
% f1bsf = lower passband cutoff frequency in Hz. f1bsf < f3bsf 
% f2bsf = upper passband cutoff frequency in Hz. f4bsf < f2bsf
% f3bsf = lower stopband cutoff frequency in Hz. f1bsf < f3bsf < f4bsf
% f4bsf = upper stopband cutoff frequency in Hz. f3bsf < f4bsf < f2bsf 
% f1bsf < f3bsf < f4bsf < f2bsf
% fs = sampling rate in Hz.
% dp = location of digital poles.
% dz = location of digital zeros.
% H = Discrete-Time Fourier Transform (DTFT)
% h = impulse response.
% References:
% M. E. Van Valkenburg, Analog Filter Design, Holt, Rinehart and Winston, Inc., 1982.
% Richard W. Daniels, Approximation Methods for Electronic Filter Design, McGraw-Hill, 1974.
% Andreas Antoniou, Digital Filters: Analysis, Design, and Applications, 2nd ed., McGraw-Hill, 1993.
% 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.
% 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[dp,dz,h] = dicbsf(amax,amin,f1bsf,f2bsf,f3bsf,f4bsf,fs)
Ts = 1/fs;     % Ts = sampling interval.
q1bsf = 2*pi*f1bsf*Ts;     % q1bsf = digital lower pass band cutoff frequency in radians.
q2bsf = 2*pi*f2bsf*Ts;     % q2bsf = digital upper pass band cutoff frequency in radians.
q3bsf = 2*pi*f3bsf*Ts;     % q3bsf = digital lower stop band cutoff frequency in radians.
q4bsf = 2*pi*f4bsf*Ts;     % q4bsf = digital upper stop band cutoff frequency in radians.
w1bsf = 2*fs*tan(q1bsf/2);   % w1bsf = prewarped lower pass band cutoff frequency in radians/s.
w2bsf = 2*fs*tan(q2bsf/2);   % w2bsf = prewarped upper pass band cutoff frequency in radians/s.
w3bsf = 2*fs*tan(q3bsf/2);   % w3bsf = prewarped lower stop band cutoff frequency in radians/s.
w4bsf = 2*fs*tan(q4bsf/2);   % w4bsf = prewarped upper stop band cutoff frequency in radians/s.
if w1bsf*w2bsf > w3bsf*w4bsf w4bsf = w1bsf*w2bsf/w3bsf; end
if w1bsf*w2bsf < w3bsf*w4bsf w3bsf = w1bsf*w2bsf/w4bsf; end
N = ceil(acosh(sqrt((10^(0.1*amin)-1)/(10^(0.1*amax)-1)))/acosh((w2bsf-w1bsf)/(w4bsf-w3bsf))); % N = order of the filter.
disp('Digital Inverse Chebyshev (Chebyshev Type II) Bandstop Filter Design')
disp(' ')
disp('Lower passband cutoff frequency q1bsf =')
disp(q1bsf)
disp('Upper passband cutoff frequency q2bsf =')
disp(q2bsf)
disp('Lower stopband cutoff frequency q3bsf =')
disp(q3bsf)
disp('Upper stopband cutoff frequency q4bsf =')
disp(q4bsf)
disp('Prewarped lower passband cutoff frequency w1bsf =')
disp(w1bsf)
disp('Prewarped upper passband cutoff frequency w2bsf =')
disp(w2bsf)
disp('Prewarped lower stopband cutoff frequency w3bsf =')
disp(w3bsf)
disp('Prewarped upper stopband cutoff frequency w4bsf =')
disp(w4bsf)
disp('The order of the filter is')
disp(N)
%
% Calculating the normalized lowpass pole locations.  
% Vector s is the normalized poles.
%
a = (1/N)*asinh(sqrt(10^(0.1*amin)-1));
if rem(N,2) == 1 s(1) = -1/sinh(a); end
if N == 1 max = 0; 
   elseif rem(N,2) == 1 max = (N-3)/2;
      else max = (N-2)/2;

end
if N > 1
 for m = 0:max
   if rem(N,2) == 1
      s(2*m+2) = -1/[sinh(a)*sin((2*m+1)*pi/(2*N))+j*cosh(a)*cos((2*m+1)*pi/(2*N))];
      s(2*m+3) = conj(s(2*m+2));
   else
      s(2*m+1) = -1/[sinh(a)*sin((2*m+1)*pi/(2*N))+j*cosh(a)*cos((2*m+1)*pi/(2*N))];
      s(2*m+2) = conj(s(2*m+1));      
   end
 end
end
%
% Calculating the normalized lowpass zero locations.  
% Vector sz is the normalized zeros.
%
if N == 1 r = 1;
   elseif rem(N,2) == 1 
      r = (N-1)/2;
   else
      r = N/2;

end
if N == 1 sz(1) = 10^50;
else
      for p = 1:r
         sz(2*p-1) = j*sec(pi*(2*p-1)/(2*N));
         sz(2*p) = conj(sz(2*p-1));
      end
end
if rem(N,2) == 1 sz(N) = 10^50; end
%
% Plotting the normalized analog lowpass poles and zeros.
%
for n = 1:201
   x(n) = (n-101)/100;
   ty(n) = sqrt(1-x(n).*x(n));
   by(n) = - ty(n);
end
figure
plot(x,ty)
axis equal;
hold on
plot(x,by)
plot(real(s),imag(s),'bx');grid;
title('Normalized Lowpass Poles')
xlabel('Real Part')
ylabel('Imaginary Part')
figure
plot(x,ty)
axis equal;
hold on
plot(x,by)
plot(real(s),imag(s),'bx')
plot(real(sz),imag(sz),'ro');grid;
title('Normalized Lowpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Normalized Lowpass Pole Locations')
disp(s)
disp('Normalized Lowpass Zero Locations')
for i=1:N-1
disp(sz(i));
end
disp(sz(N));
%
% Calculating the frequency transformed analog poles and zeros.
% Lowpass to bandstop transformation.
%
Wo=(w2bsf-w1bsf)/(w4bsf-w3bsf);
w21=(w2bsf-w1bsf)/2;
w2m1=w2bsf-w1bsf;
w12t4=4*w1bsf*w2bsf;
for k = 1:N
   s(k) = Wo*s(k);
   R1(k) = real(s(k));
   I1(k) = imag(s(k));
   R2(k) = R1(k)^2;
   I2(k) = I1(k)^2;
   mag(k) = R2(k) + I2(k);
   A(k) = (I2(k)-R2(k))*w2m1^2/mag(k)^2 + w12t4;
   B(k) = 2*abs(R1(k))*abs(I1(k))*w2m1^2/mag(k)^2;
   RA(k) = (1/2)*sqrt((sqrt(A(k)^2+B(k)^2)-A(k))/2);
   RB(k) = j*(1/2)*(sqrt((sqrt(A(k)^2+B(k)^2)+A(k))/2)-I1(k)*w2m1/mag(k));
end
for m = 1:2*N
   if rem(m,2) == 1
      if I1((m+1)/2) > 0 sfs(m) = w21*R1((m+1)/2)/mag((m+1)/2) + RA((m+1)/2) + RB((m+1)/2); end
      if I1((m+1)/2) < 0 sfs(m) = w21*R1((m+1)/2)/mag((m+1)/2) - RA((m+1)/2) + RB((m+1)/2); end    
      R1a = (w2m1/R1((m+1)/2))^2;
      if I1((m+1)/2) == 0 & R1a >= w12t4 sfs(m) = w21/R1((m+1)/2) + (1/2)*sqrt(R1a-w12t4); end  
      if I1((m+1)/2) == 0 & R1a < w12t4 sfs(m) = w21/R1((m+1)/2) + j*(1/2)*sqrt(-R1a+w12t4); end     
   else
      sfs(m) = conj(sfs(m-1));
      if I1(m/2) == 0 & (w2m1/R1(m/2))^2 >= w12t4  sfs(m) = w21/R1(m/2) - (1/2)*sqrt(w2m1^2/R1(m/2)^2-w12t4); end  
      if I1(m/2) == 0 & (w2m1/R1(m/2))^2 < w12t4  sfs(m) = w21/R1(m/2) - j*(1/2)*sqrt(-w2m1^2/R1(m/2)^2+w12t4); end  
   end
end
% Transformation of zeros
for k = 1:N
   sz(k) = Wo*sz(k);
   Rz1(k) = real(sz(k));
   Iz1(k) = imag(sz(k));
   Rz2(k) = Rz1(k)^2;
   Iz2(k) = Iz1(k)^2;
   mag(k) = Rz2(k) + Iz2(k);
   Az(k) = (Iz2(k)-Rz2(k))*w2m1^2/mag(k)^2 + w12t4;
   Bz(k) = 2*abs(Rz1(k))*abs(Iz1(k))*w2m1^2/mag(k)^2;
   RAz(k) = (1/2)*sqrt((sqrt(Az(k)^2+Bz(k)^2)-Az(k))/2);
   RBz(k) = j*(1/2)*(sqrt((sqrt(Az(k)^2+Bz(k)^2)+Az(k))/2)-Iz1(k)*w2m1/mag(k));
   if Iz1(k) ~= 0 AAz(k) = ((w2bsf-w1bsf)/Iz1(k))^2 + w12t4; end
end
for m = 1:2*N
   if rem(m,2) == 1
      if Iz1((m+1)/2) >= 0 szfs(m) = w21*Rz1((m+1)/2)/mag((m+1)/2) + RAz((m+1)/2) + RBz((m+1)/2); end
      if Iz1((m+1)/2) < 0 szfs(m) = w21*Rz1((m+1)/2)/mag((m+1)/2) - RAz((m+1)/2) + RBz((m+1)/2); end    
      if Iz1((m+1)/2) > 0 & Rz1((m+1)/2) == 0 szfs(m) = j*(1/2)*(sqrt(AAz((m+1)/2))-w2m1/Iz1((m+1)/2)); end  
      if Iz1((m+1)/2) < 0 & Rz1((m+1)/2) == 0 szfs(m) = j*(1/2)*(sqrt(AAz((m+1)/2))-w2m1/Iz1((m+1)/2)); end  
   else
      szfs(m) = conj(szfs(m-1));
   end
end
%
% Plotting the frequency transformed bandpass poles and zeros.
%
wo = sqrt(w1bsf*w2bsf)/Wo;
for na = 1:201
   xa(na) = wo*(na-101)/100;
   tya(na) = sqrt(wo^2-xa(na)^2);
   bya(na) = - tya(na);
end
figure
plot(xa,tya)
axis equal;
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx');grid;
plot(real(szfs),imag(szfs),'ro')
title('Frequency Transformed Bandstop Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Frequency Transformed Bandstop Pole Locations')
disp(sfs)
disp('Frequency Transformed Bandstop Zero Locations')
disp(szfs)
%
% Transforming analog poles and zeros to digital poles and zeros using bilinear transformation.
%
for m = 1:2*N
   dp(m) = (1 + (Ts/2)*sfs(m))/(1 - (Ts/2)*sfs(m));
   dz(m) = (1 + (Ts/2)*szfs(m))/(1 - (Ts/2)*szfs(m));
end
%
% Plotting digital poles and zeros on the complex z-plane
%
figure
plot(x,ty)
axis equal;
hold on
plot(x,by)
plot(real(dp),imag(dp),'bx');grid;
plot(real(dz),imag(dz),'ro')
title('Digital Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Digital Pole Locations')
disp(dp)
disp('Digital Zero Locations')
disp(dz)
%
for i = 1:N
	d2(i,1)=1; d2(i,2)=-2*real(dp(2*i-1));d2(i,3) = (abs(dp(2*i-1)))^2;
	n2(i,1)=1; n2(i,2)=-2*real(dz(2*i-1));n2(i,3) = (abs(dz(2*i-1)))^2;
end
disp('Numerator Second Order Terms')
disp(n2)
disp('Denominator Second Order Terms')
disp(d2)
%
% Plotting the magnitude response and the phase response with respect to the normalized digital frequency.
%
M = 128;  % M is the number of points the frequency response is calculated.
step = 2*pi/M;
for m = 1:M
   H(m) = (1-dz(1)*exp(j*(-1)*(m-1)*step))*(1-dz(2)*exp(j*(-1)*(m-1)*step))/[(1-dp(1)*exp(j*(-1)*(m-1)*step))*(1-dp(2)*exp(j*(-1)*(m-1)*step))];
end
if N > 1
   for n = 3:2*N
     for m = 1:M
       H(m) = H(m).*[(1-dz(n)*exp(j*(-1)*(m-1)*step))/(1-dp(n)*exp(j*(-1)*(m-1)*step))];
     end
   end
end
max = abs(H(1));
for m = 2:M
   if abs(H(m)) > max max = abs(H(m)); end
end
disp('maximum gain is')
disp(max)
b0=1/max;
disp('Constant Term b0 =')
disp(b0);
b0n=b0^(1/N);
disp('Constant Term per 2nd-order term b0n =')
disp(b0n);
%
%Polynomial Form
%
nump(1) = n2(1,1); nump(2) = n2(1,2); nump(3) = n2(1,3); 
denp(1) = d2(1,1); denp(2) = d2(1,2); denp(3) = d2(1,3);
if N > 1
	for i=2:N nump = conv(nump,n2(i,:)); denp = conv(denp,d2(i,:)); end
end
disp('Numerator Polynomial')
disp(b0*nump)
disp('Denominator Polynomial')
disp(denp)
%
%
for i = 1:N

	n2(i,1)=b0n*n2(i,1); n2(i,2)=b0n*n2(i,2);n2(i,3) = b0n*n2(i,3);
end
disp('Numerator Second Order Terms Multiplied by b0n')
disp(n2)
disp('Denominator Second Order Terms')
disp(d2)
%
for m = 1:M
   f(m) = (m-1)/M; 
   H(m) = H(m)/max;
   mag(m) = abs(H(m));
   if H(m) == 0 dB(m) = -200; else dB(m) = 20*log10(mag(m)); end
   if dB(m) < -200 dB(m) = -200; end
   pha(m) = angle(H(m));
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 the Impulse Response of the Filter.
%
h = ifft(H);
for m = 1:M
   k(m) = m-1; 
   z(m) = 0;
   if real(h(m)) >= 0 Lo(m) = 0; else Lo(m) = -real(h(m)); end
   if real(h(m)) >= 0 Hi(m) = real(h(m)); else Hi(m) = 0; end
end
figure
%errorbar(k,z,Lo,Hi,'y')
stem(k,h);grid;
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')