% m-file file name: dicbpf.m
% Function for designing digital inverse Chevyshev bandpass filter.
% Specifications are given by amax, amin, f1bpf, f2bpf, f3bpf, f4bpf.
% 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.
% f1bpf = lower passband cutoff frequency in Hz. f3bpf < f1bpf < f2bpf
% f2bpf = upper passband cutoff frequency in Hz. f1bpf < f2bpf < f4bpf
% f3bpf = lower stopband cutoff frequency in Hz. f3bpf < f1bpf 
% f4bpf = upper stopband cutoff frequency in Hz. f2bpf < f4bpf 
% 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.
% Lonnie C. Ludeman, Fundamentals of Digital Signal Processing, Wiley, 1986.
% 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,h] = dicbpf(amax,amin,f1bpf,f2bpf,f3bpf,f4bpf,fs)
Ts = 1/fs;     % Ts = sampling interval.
q1bpf = 2*pi*f1bpf*Ts;     % q1bpf = digital lower pass band cutoff frequency in radians.
q2bpf = 2*pi*f2bpf*Ts;     % q2bpf = digital upper pass band cutoff frequency in radians.
q3bpf = 2*pi*f3bpf*Ts;     % q3bpf = digital lower stop band cutoff frequency in radians.
q4bpf = 2*pi*f4bpf*Ts;     % q4bpf = digital upper stop band cutoff frequency in radians.
w1bpf = 2*fs*tan(q1bpf/2);   % w1bpf = prewarped lower pass band cutoff frequency in radians/s.
w2bpf = 2*fs*tan(q2bpf/2);   % w2bpf = prewarped upper pass band cutoff frequency in radians/s.
w3bpf = 2*fs*tan(q3bpf/2);   % w3bpf = prewarped lower stop band cutoff frequency in radians/s.
w4bpf = 2*fs*tan(q4bpf/2);   % w4bpf = prewarped upper stop band cutoff frequency in radians/s.
if w1bpf*w2bpf < w3bpf*w4bpf w4bpf = w1bpf*w2bpf/w3bpf; end
q4abpf = 2*atan(w4bpf*Ts/2);  %q4abpf = q4bpf corresponding to adjusted w4bpf 
if w1bpf*w2bpf > w3bpf*w4bpf w3bpf = w1bpf*w2bpf/w4bpf; end
q3abpf = 2*atan(w3bpf*Ts/2);   %q3abpf = q3bpf corresponding to adjusted w3bpf 
wobpf = sqrt(w1bpf*w2bpf); %wobpf = geometric mean of w1bpf and w2bpf.
N = ceil(acosh(sqrt((10^(0.1*amin)-1)/(10^(0.1*amax)-1)))/acosh((w4bpf-w3bpf)/(w2bpf-w1bpf))); % N = order of the filter.
disp('Digital Inverse Chebyshev (Chebyshev Type II) Bandpass Filter Design')
disp(' ')
disp('Lower passband cutoff frequency q1bpf =')
disp(q1bpf)
disp('Upper passband cutoff frequency q2bpf =')
disp(q2bpf)
disp('Lower stopband cutoff frequency q3bpf =')
disp(q3bpf)
disp('Upper stopband cutoff frequency q4bpf =')
disp(q4bpf)
disp('Prewarped lower passband cutoff frequency w1bpf =')
disp(w1bpf)
disp('Prewarped upper passband cutoff frequency w2bpf =')
disp(w2bpf)
disp('Prewarped lower stopband cutoff frequency w3bpf =')
disp(w3bpf)
disp('Prewarped upper stopband cutoff frequency w4bpf =')
disp(w4bpf)
disp('Center frequency = sqrt(w1bpf*w2bpf) = wobpf =')
disp(wobpf)
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
   %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
   %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
%
disp('Normalized Lowpass Pole Locations')
disp(s)
disp('Normalized Lowpass Zero Locations')
disp(sz(1))
for i = 2:N
	disp(sz(i))
end
%
%Adjustments
%
Wca = cosh((1/N)*acosh(sqrt(10^(0.1*amin)-1)));
for i = 1:N
    s(i) = Wca*s(i);
    sz(i) = Wca*sz(i);
end
disp('Normalized Lowpass Pole Locations After Adjustments')
disp(s)
disp('Normalized Lowpass Zero Locations After Adjustments')
disp(sz(1))
for i = 2:N
	disp(sz(i))
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)
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)
hold on
plot(x,by)
plot(real(s),imag(s),'bx')
if rem(N,2)==0
   for i=1:N
      szp(i)=sz(i);
   end
else
   for i=1:N-1
      szp(i)=sz(i);
   end
end   
plot(real(szp),imag(szp),'ro'); grid;
title('Normalized Lowpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
%
% Calculating the frequency transformed analog poles and zeros.
% Lowpass to bandpass transformation.
%
Wo = (w4bpf-w3bpf)/(w2bpf-w1bpf); 
w21 = (w2bpf-w1bpf)/2;
wcon = 4*w1bpf*w2bpf/(w2bpf-w1bpf)^2;
for k = 1:N
   R1(k) = real(s(k));
   I1(k) = imag(s(k));
   R2(k) = R1(k)^2;
   I2(k) = I1(k)^2;
   A(k) = wcon + I2(k) - R2(k);
   Ra(k) = sqrt((sqrt(A(k)^2 + 4*R2(k)*I2(k)) - A(k))/2);
   Rb(k) = sqrt((sqrt(A(k)^2 + 4*R2(k)*I2(k)) + A(k))/2);
end
for m = 1:2*N
   if rem(m,2) == 1
      if I1((m+1)/2) >= 0 sfs(m) = w21*(R1((m+1)/2) + Ra((m+1)/2)) + j*w21*(I1((m+1)/2) - Rb((m+1)/2)); end
      if I1((m+1)/2) < 0 sfs(m) = w21*(R1((m+1)/2) - Ra((m+1)/2)) + j*w21*(-I1((m+1)/2) + Rb((m+1)/2)); end      
      if I1((m+1)/2) == 0 sfs(m) = w21*R1((m+1)/2) + j*w21*sqrt(wcon - R1((m+1)/2)^2); end  
      if I1((m+1)/2) > 0 & R1((m+1)/2) == 0 sfs(m) = j*w21*(sqrt(wcon + I1((m+1)/2)^2) - I1((m+1)/2)); end  
      if I1((m+1)/2) < 0 & R1((m+1)/2) == 0 sfs(m) = j*w21*(sqrt(wcon + I1((m+1)/2)^2) - I1((m+1)/2)); end  
   else
      sfs(m) = conj(sfs(m-1));
   end
end
% Transformation of zeros
for k = 1:N
   Rz1(k) = real(sz(k));
   Iz1(k) = imag(sz(k));
   Rz2(k) = Rz1(k)^2;
   Iz2(k) = Iz1(k)^2;
   Az(k) = wcon + Iz2(k) - Rz2(k);
   Raz(k) = sqrt((sqrt(Az(k)^2 + 4*Rz2(k)*Iz2(k)) - Az(k))/2);
   Rbz(k) = sqrt((sqrt(Az(k)^2 + 4*Rz2(k)*Iz2(k)) + Az(k))/2);
end
for m = 1:2*N
   if rem(m,2) == 1  
      if Iz1((m+1)/2) >= 0 szfs(m) = w21*(Rz1((m+1)/2)+Raz((m+1)/2))+j*w21*(Iz1((m+1)/2)-Rbz((m+1)/2)); end
      if Iz1((m+1)/2) < 0 szfs(m) = w21*(Rz1((m+1)/2)-Raz((m+1)/2))+j*w21*(-Iz1((m+1)/2)+Rbz((m+1)/2)); end
      if Rz1((m+1)/2) == 0 & Iz1((m+1)/2) > 0 szfs(m) = j*w21*(sqrt(wcon+Iz2((m+1)/2))-Iz1((m+1)/2)); end
      if Rz1((m+1)/2) == 0 & Iz1((m+1)/2) < 0 szfs(m) = j*w21*(sqrt(wcon+Iz2((m+1)/2))-Iz1((m+1)/2)); end
   else
      szfs(m) = conj(szfs(m-1));
   end
end
if rem(N,2) == 1 szfs(2*N-1) = 0; end
if rem(N,2) == 1 szfs(2*N) = 10^50; end
%
% Plotting the frequency transformed bandpass poles and zeros.
%
wo = sqrt(w1bpf*w2bpf)/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)
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx'); grid;
title('Frequency Transformed Bandpass Poles')
xlabel('Real Part')
ylabel('Imaginary Part')
figure
plot(xa,tya)
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx')
hold on
if rem(N,2)==0
   for i=1:2*N
      szfsp(i)=szfs(i);
   end
else
   for i=1:2*N-1
      szfsp(i)=szfs(i);
   end
end 
plot(real(szfsp),imag(szfsp),'ro'); grid;
title('Frequency Transformed Bandpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Frequency Transformed Bandpass Pole Locations')
disp(sfs)
disp('Frequency Transformed Bandpass Zero Locations')

for i = 1:2*N-1
	disp(szfs(i));
end
disp(szfs(2*N));
%
% 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)
hold on
plot(x,by)
plot(real(dp),imag(dp),'bx')
plot(real(dz),imag(dz),'ro'); grid;
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
if rem(N,2)==1 n2(N,1)=1; n2(N,2)=0; n2(N,3)=-1; 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)
b01=1/max;
disp('1/max')
disp(b01);
%
%Finding b0
%
q0 = 2*atan(wobpf*Ts/2);
if N==1 
    H0 = (exp(j*q0) - dz(1))*(exp(j*q0) - dz(2))/((exp(j*q0) - dp(1))*(exp(j*q0) - dp(2)));
else
    H0 = (exp(j*q0) - dz(1))/(exp(j*q0) - dp(1)); 
    for i = 2:2*N
        H0 = H0*(exp(j*q0) - dz(i))/(exp(j*q0) - dp(i));
    end
end
if rem(N,2) == 1
    b0 = 1/real(H0);
else
    b0 = 10^(-0.05*amax)/real(H0);
end
disp('Constant Term b0 =')
disp(b0);
%
%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 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')
%
%Verifying
%
    Hq1 =abs(1 - dz(1)*exp(-j*q1bpf))/abs(1 - dp(1)*exp(-j*q1bpf));
    Hq2 =abs(1 - dz(1)*exp(-j*q2bpf))/abs(1 - dp(1)*exp(-j*q2bpf));
    Hq3 =abs(1 - dz(1)*exp(-j*q3bpf))/abs(1 - dp(1)*exp(-j*q3bpf));
    Hq3a =abs(1 - dz(1)*exp(-j*q3abpf))/abs(1 - dp(1)*exp(-j*q3abpf));
    Hq4 =abs(1 - dz(1)*exp(-j*q4bpf))/abs(1 - dp(1)*exp(-j*q4bpf));  
    Hq4a =abs(1 - dz(1)*exp(-j*q4abpf))/abs(1 - dp(1)*exp(-j*q4abpf));
    for p = 2:2*N 
        Hq1 =Hq1*abs(1 - dz(p)*exp(-j*q1bpf))/abs(1 - dp(p)*exp(-j*q1bpf));
        Hq2 =Hq2*abs(1 - dz(p)*exp(-j*q2bpf))/abs(1 - dp(p)*exp(-j*q2bpf));
        Hq3 =Hq3*abs(1 - dz(p)*exp(-j*q3bpf))/abs(1 - dp(p)*exp(-j*q3bpf));
        Hq3a =Hq3*abs(1 - dz(p)*exp(-j*q3abpf))/abs(1 - dp(p)*exp(-j*q3abpf));
        Hq4 =Hq4*abs(1 - dz(p)*exp(-j*q4bpf))/abs(1 - dp(p)*exp(-j*q4bpf));
        Hq4a =Hq4a*abs(1 - dz(p)*exp(-j*q4abpf))/abs(1 - dp(p)*exp(-j*q4abpf));
    end
    Hq1 = -20*log10(abs(b0*Hq1));
    Hq2 = -20*log10(abs(b0*Hq2));
    Hq3 = -20*log10(abs(b0*Hq3));
    Hq3a = -20*log10(abs(b0*Hq3a));
    Hq4 = -20*log10(abs(b0*Hq4));
    Hq4a = -20*log10(abs(b0*Hq4a));
disp('Loss in dB at the lower passband cutoff frequency (q1bpf)')
disp(Hq1)
disp('Loss in dB at the upper passband cutoff frequency (q2bpf)')
disp(Hq2)
disp('Loss in dB at the lower stopband cutoff frequency (q3bpf)')
disp(Hq3)
disp('Loss in dB at the upper stopband cutoff frequency (q4bpf)')
disp(Hq4)
if q3abpf ~= q3bpf
   disp('Loss in dB at the lower stopband cutoff frequency adjusted (q3abpf)')
   disp(Hq3a) 
end 
if q4abpf ~= q4bpf
   disp('Loss in dB at the upper stopband cutoff frequency adjusted (q4abpf)')
   disp(Hq4a) 
end