% m-file file name: delbsf.m
% Function for designing digital elliptic bandstop 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.
% f1bsf = lower pass band cutoff frequency in Hz. f1bsf < f3bsf 
% f2bsf = upper pass band cutoff frequency in Hz. f4bsf < f2bsf
% f3bsf = lower stop band cutoff frequency in Hz. f1bsf < f3bsf < f4bsf
% f4bsf = upper stop band 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:
% 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] = delbsf(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.
Ws1 = (w3bsf*(w2bsf-w1bsf))/(-w3bsf^2 + w1bsf*w2bsf); %backward transformation.
Ws2 = (w4bsf*(w2bsf-w1bsf))/(-w4bsf^2 + w1bsf*w2bsf); %backward transformation.
Wslpf = min(abs(Ws1),abs(Ws2)); %Wslpf = stopband cutoff frequency of the normalized lowpass fiter.
if w1bsf*w2bsf > w3bsf*w4bsf w4bsf = w1bsf*w2bsf/w3bsf; end
q4absf = 2*atan(w4bsf*Ts/2);  %q4absf = q4bsf corresponding to adjusted w4bsf
if w1bsf*w2bsf < w3bsf*w4bsf w3bsf = w1bsf*w2bsf/w4bsf; end
q3absf = 2*atan(w3bsf*Ts/2);   %q3absf = q3bsf corresponding to adjusted w3bsf 
wobsf = sqrt(w1bsf*w2bsf); %wobsf = geometric mean of w1bsf and w2bsf.
%
% Calculating the order of the normalized lowpass filter.
%
k = (w4bsf-w3bsf)/(w2bsf-w1bsf);
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 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('Center frequency = sqrt(w1bsf*w2bsf) = wobsf =')
disp(wobsf)
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;
for i = 1:N
 if real(sz(i)) < 10^3
   plot(real(sz(i)),imag(sz(i)),'ro')
 end
end
title('Normalized Lowpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
%
% Calculating the frequency transformed analog poles and zeros.
% Lowpass to bandstop transformation.
%
Wo=sqrt((w2bsf-w1bsf)/(w4bsf-w3bsf));
w21=(w2bsf-w1bsf)/2;
w2m1=w2bsf-w1bsf;
w12t4=4*w1bsf*w2bsf;
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;
   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) = 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)
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx');grid;
for i = 1:2*N
 if real(szfs(i)) < 10^30
    plot(real(szfs(i)),imag(szfs(i)),'ro')
 end
end
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)
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
    if rem(N,2) == 1 
        if i == 1
            d2(i,1) = 1; d2(i,2) = -dp(1)-dp(2); d2(i,3) = (-dp(1))*(-dp(2));
            n2(i,1) = 1; n2(i,2) = -dz(1)-dz(2); n2(i,3) = (-dz(1))*(-dz(2));
        else           
	        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
    else
	        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
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
%

    H0 = (1 - dz(1))/(1 - dp(1)); 
    for i = 2:2*N
        H0 = H0*(1 - dz(i))/(1 - dp(i));
    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*q1bsf))/abs(1 - dp(1)*exp(-j*q1bsf));
    Hq2 =abs(1 - dz(1)*exp(-j*q2bsf))/abs(1 - dp(1)*exp(-j*q2bsf));
    Hq3 =abs(1 - dz(1)*exp(-j*q3bsf))/abs(1 - dp(1)*exp(-j*q3bsf));
    Hq3a =abs(1 - dz(1)*exp(-j*q3absf))/abs(1 - dp(1)*exp(-j*q3absf));
    Hq4 =abs(1 - dz(1)*exp(-j*q4bsf))/abs(1 - dp(1)*exp(-j*q4bsf));  
    Hq4a =abs(1 - dz(1)*exp(-j*q4absf))/abs(1 - dp(1)*exp(-j*q4absf));
    for p = 2:2*N 
        Hq1 =Hq1*abs(1 - dz(p)*exp(-j*q1bsf))/abs(1 - dp(p)*exp(-j*q1bsf));
        Hq2 =Hq2*abs(1 - dz(p)*exp(-j*q2bsf))/abs(1 - dp(p)*exp(-j*q2bsf));
        Hq3 =Hq3*abs(1 - dz(p)*exp(-j*q3bsf))/abs(1 - dp(p)*exp(-j*q3bsf));
        Hq3a =Hq3*abs(1 - dz(p)*exp(-j*q3absf))/abs(1 - dp(p)*exp(-j*q3absf));
        Hq4 =Hq4*abs(1 - dz(p)*exp(-j*q4bsf))/abs(1 - dp(p)*exp(-j*q4bsf));
        Hq4a =Hq4a*abs(1 - dz(p)*exp(-j*q4absf))/abs(1 - dp(p)*exp(-j*q4absf));
    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 (q1bsf)')
disp(Hq1)
disp('Loss in dB at the upper passband cutoff frequency (q2bsf)')
disp(Hq2)
disp('Loss in dB at the lower stopband cutoff frequency (q3bsf)')
disp(Hq3)
disp('Loss in dB at the upper stopband cutoff frequency (q4bsf)')
disp(Hq4)
if q3absf ~= q3bsf
   disp('Loss in dB at the lower stopband cutoff frequency adjusted (q3absf)')
   disp(Hq3a) 
end 
if q4absf ~= q4bsf
   disp('Loss in dB at the upper stopband cutoff frequency adjusted (q4absf)')
   disp(Hq4a) 
end