% m-file file name: dchbpf.m
% Function for designing digital 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.
% 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] = dchbpf(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
if w1bpf*w2bpf > w3bpf*w4bpf w3bpf = w1bpf*w2bpf/w4bpf; end
w0=sqrt(w1bpf*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 Chebyshev (Type I) 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('Resonant frequency w0 = sqrt(w1bpf*w2bpf) =')
disp(w0)
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(1/sqrt(10^(0.1*amax)-1));
if rem(N,2) == 1 s(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) = -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) = -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
%
% Plotting the normalized analog lowpass poles.
%
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 Analog Lowpass Pole Locations')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Normalized Analog Lowpass Pole Locations')
disp(s)
%
% Calculating the normalized analog lowpass zero locations.  
% Vector szb is the normalized analog zeros.
%
for k = 1:N
   sz(k) = 10^50;
end
%
% Calculating the frequency transformed analog poles and zeros.
% Lowpass to bandpass transformation.
%
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
for m = 1:2*N
   if rem(m,2) == 1 szfs(m) = 0; else szfs(m) = 10^50; end
end
%
% Plotting the frequency transformed bandpass poles.
%
wo = sqrt(w1bpf*w2bpf);
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 Pole Locations')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Frequency Transformed Bandpass Pole Locations')
disp(sfs)
disp('Frequency Transformed Bandpass 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)
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)=0; n2(i,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 = 512;  % 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 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')