% m-file file name: delbpf.m
% Function for designing digital elliptic bandpass filter.
% Written by Dr. James S. Kang, Professor 
% Department of Electrical and Computer Engineering
% California State Polytechnic University, Pomona 
% 3801 W. Temple Ave, Pomona, CA 91768
% jskang@csupomona.edu
% 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:
% Andreas Antoniou, Digital Filters: Analysis, Design, and Applications, 2nd ed., McGraw-Hill, 1993.
% 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.
% 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] = delbpf(amax,amin,f1bpf,f2bpf,f3bpf,f4bpf,fs)
Ts = 1/fs;     % Ts = sampling interval.
q1bpf = 2*pi*f1bpf*Ts;     % q1bpf = digital lower passband cutoff frequency in radians.
q2bpf = 2*pi*f2bpf*Ts;     % q2bpf = digital upper passband cutoff frequency in radians.
q3bpf = 2*pi*f3bpf*Ts;     % q3bpf = digital lower stopband cutoff frequency in radians.
q4bpf = 2*pi*f4bpf*Ts;     % q4bpf = digital upper stopband cutoff frequency in radians.
w1bpf = 2*fs*tan(q1bpf/2);   % w1bpf = prewarped lower passband cutoff frequency in radians/s.
w2bpf = 2*fs*tan(q2bpf/2);   % w2bpf = prewarped upper passband cutoff frequency in radians/s.
w3bpf = 2*fs*tan(q3bpf/2);   % w3bpf = prewarped lower stopband cutoff frequency in radians/s.
w4bpf = 2*fs*tan(q4bpf/2);   % w4bpf = prewarped upper stopband cutoff frequency in radians/s.
Ws1 = (w3bpf^2 - w1bpf*w2bpf)/(w3bpf*(w2bpf-w1bpf)); %backward transformation.
Ws2 = (w4bpf^2 - w1bpf*w2bpf)/(w4bpf*(w2bpf-w1bpf)); %backward transformation.
Wslpf = min(abs(Ws1),abs(Ws2)); %Wslpf = stopband cutoff frequency of the normalized lowpass fiter.
q3abpf = q3bpf;
q4abpf = q4bpf;
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.
%
% Finding the order of the filter.
%
k = (w2bpf-w1bpf)/(w4bpf-w3bpf);
Wp = sqrt(k);
Ws = 1/sqrt(k);
kp = sqrt(1-k^2);
qo = (1/2)*(1-sqrt(kp))/(1+sqrt(kp));
q = qo + 2*qo^5 + 15*qo^9 + 150*qo^13;
D = (10^(0.1*amin)-1)/(10^(0.1*amax)-1);
N = ceil(log10(16*D)/log10(1/q)); % N = order of the filter.
disp('Digital Elliptic 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.
%
L = (1/(2*N))*log((10^(0.05*amax)+1)/(10^(0.05*amax)-1));
sumn = sinh(L);
for m = 1:50
   sumn = sumn + (-1)^m*q^(m*(m+1))*sinh((2*m+1)*L);
end
sumd = -q*cosh(2*L);
for m = 2:50
   sumd = sumd + (-1)^m*q^m^2*cosh(2*m*L);
end
so = abs(2*q^0.25*sumn/(1+2*sumd));
W = sqrt((1+k*so^2)*(1+so^2/k));
if N == 1 r = 1;
   elseif rem(N,2) == 1
      r = (N-1)/2;
   else 
      r = N/2;
      %   end
end
if rem(N,2) == 1 muminus = 0; else muminus = 0.5; end
for p = 1:r
   sumna = sin(pi*(p-muminus)/N);
   for m = 1:50
      sumna = sumna + (-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*(p-muminus)/N);            
   end
   sumda = -q*cos(2*pi*(p-muminus)/N);
   for m = 2:50
      sumda = sumda + (-1)^m*q^m^2*cos(2*m*pi*(p-muminus)/N);
   end
   WG(p) = 2*q^0.25*sumna/(1+2*sumda);
   V(p) = sqrt((1-k*WG(p)^2)*(1-WG(p)^2/k));
   if N == 1 WG(1) = 10^(-50); end
   Ao(p) = 1/WG(p)^2;
   Bo(p) = [(so*V(p))^2 + (WG(p)*W)^2]/[1+so^2*WG(p)^2]^2;
   B1(p) = 2*so*V(p)/[1+so^2*WG(p)^2];
end
if rem(N,2) == 1 Ho = so*Bo(1)/Ao(1); end
if rem(N,2) == 0 Ho = 10^(-0.05*amax)*Bo(p)/Ao(p); end
if N > 1
   for p = 2:r
      Ho = Ho*Bo(p)/Ao(p);
   end
end
if rem(N,2) == 1 s(1) = -so; end
if N > 1
   if rem(N,2) == 1
      for m = 2:2*r+1
         if rem(m,2) == 0
            if B1(m/2)^2 >= 4*Bo(m/2) s(m) = -B1(m/2)/2+0.5*sqrt(B1(m/2)^2-4*Bo(m/2));
            else
            s(m) = -B1(m/2)/2+j*0.5*sqrt(-B1(m/2)^2+4*Bo(m/2));  
            end 
         end
         if rem(m,2) == 1
            if B1((m-1)/2)^2 >= 4*Bo((m-1)/2) s(m) = -B1((m-1)/2)/2-0.5*sqrt(B1((m-1)/2)^2-4*Bo((m-1)/2));
            else
            s(m) = conj(s(m-1));
            end
         end            
      end
   else
      for m = 1:2*r
         if rem(m,2) == 1
            if B1((m+1)/2)^2 >= 4*Bo((m+1)/2) s(m) = -B1((m+1)/2)/2+0.5*sqrt(B1((m+1)/2)^2-4*Bo((m+1)/2));
            else
            s(m) = -B1((m+1)/2)/2+j*0.5*sqrt(-B1((m+1)/2)^2+4*Bo((m+1)/2));   
            end
         end
         if rem(m,2) == 0
            if B1(m/2)^2 >= 4*Bo(m/2) s(m) = -B1(m/2)/2-0.5*sqrt(B1(m/2)^2-4*Bo(m/2));
            else
            s(m) = conj(s(m-1));
            end
         end            
      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 rem(N,2) == 1 sz(1) = 10^50; end
if rem(N,2) == 1 & N > 1
      for m = 1:2*r
         if rem(m,2) == 1 sz(m+1) = j*sqrt(Ao((m+1)/2));
         else
         sz(m+1) = conj(sz(m));
         end
      end
end
if rem(N,2) == 0
      for m = 1:2*r
         if rem(m,2) == 1 sz(m) = j*sqrt(Ao((m+1)/2));
         else
         sz(m) = conj(sz(m-1));
         end
      end
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
%
sqrtws = sqrt(Wslpf);
for i = 1:N
    s(i) = sqrtws*s(i);
    sz(i) = sqrtws*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 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;
if rem(N,2) == 1 & N > 1
	for i = 2:N
   		plot(real(sz(i)),imag(sz(i)),'ro');
	end
end
if rem(N,2) == 0
	plot(real(sz),imag(sz),'ro');
end
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 = sqrt((w4bpf-w3bpf)/(w2bpf-w1bpf)); 
w21 = (w2bpf-w1bpf)/2;
wcon = 4*w1bpf*w2bpf/(w2bpf-w1bpf)^2;
for k = 1:N
   s(k) = s(k);
   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
   sz(k) = sz(k);
   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(1) = 10^50; end
if rem(N,2) == 1 szfs(2) = 0; end
%
% Plotting the frequency transformed bandpass poles and zeros.
%
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 Poles')
xlabel('Real Part')
ylabel('Imaginary Part')
figure
plot(xa,tya)
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx')
hold on
grid;
if rem(N,2) == 1 & N > 1
	for i = 2:2*N
    		plot(real(szfs(i)),imag(szfs(i)),'ro');
	end
end
if rem(N,2) == 0
    		plot(real(szfs),imag(szfs),'ro');
end
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')
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'); 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;
   if ((dz(2*i-1)==1)&(dz(2*i)==-1)) | ((dz(2*i-1)==-1)&(dz(2*i)==1))
	n2(i,1)=1; n2(i,2)=0; n2(i,3) = -1;
   else
	n2(i,1)=1; n2(i,2)=-2*real(dz(2*i-1));n2(i,3) = (abs(dz(2*i-1)))^2;
   end
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('Constant Term b0 =')
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