% m-file file name: dellpf.m
% Function for designing digital elliptic lowpass 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
% amax = maximum loss in dB in the pass band.
% amin = minimum loss in dB in the stop band.
% fplpf = pass band cutoff frequency in Hz.
% fslpf = stop band cutoff frequency in Hz.
% 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]= dellpf(amax,amin,fplpf,fslpf,fs)
Ts = 1/fs;     % Ts = sampling interval.
qplpf = 2*pi*fplpf*Ts;     % qplpf = digital pass band cutoff frequency in radians.
qslpf = 2*pi*fslpf*Ts;     % qplpf = digital stop band cutoff frequency in radians.
wplpf = 2*fs*tan(qplpf/2);   % wplpf = prewarped pass band cutoff frequency in radians/s.
wslpf = 2*fs*tan(qslpf/2);   % wplpf = prewarped stop band cutoff frequency in radians/s.
%
% Calculating the order of the normalized lowpass filter.
%
wo2 = sqrt(wplpf*wslpf);
k = wplpf/wslpf;
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 Lowpass Filter Design')
disp(' ')
disp('Digital pass band cutoff frequency qplpf =')
disp(qplpf)
disp('Digital stop band cutoff frequency qplpf =')
disp(qslpf)
disp('Prewarped pass band cutoff frequency wplpf =')
disp(wplpf)
disp('Prewarped stop band cutoff frequency wplpf =')
disp(wslpf)
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)); % L = Lamda in Antoniou
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
%
sqrtsp = sqrt(wslpf/wplpf);
for i = 1:N
    s(i) = sqrtsp*s(i);
    sz(i) = sqrtsp*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')
for i = 1:N
 if real(sz(i)) < 10^3
   plot(real(sz(i)),imag(sz(i)),'ro')
 end
end
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')
disp(sz(1))
for i = 2:N
	disp(sz(i))
end
%
% Calculating the frequency transformed passband cutoff frequency.
% Lowpass to lowpass transformation.
%
wolpf = sqrt(wplpf*wslpf);  
for k = 1:N
   sfs(k) = wplpf*s(k);
   szfs(k) = wplpf*sz(k);
end
%
% Plotting the frequency transformed lowpass poles.
%
for na = 1:201
   xa(na) = wplpf*(na-101)/100;
   tya(na) = sqrt(wplpf^2-xa(na).*xa(na));
   bya(na) = - tya(na);
end
figure
plot(xa,tya)
hold on
plot(xa,bya)
plot(real(sfs),imag(sfs),'bx'); grid;
title('Frequency Transformed Lowpass Pole Locations')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Frequency Transformed Lowpass Pole Locations')
disp(sfs)
disp('Frequency Transformed Lowpass Zero Locations')
disp(szfs(1))
for i = 2:N
	disp(szfs(i))
end
%
% Transforming analog poles and zeros to digital poles and zeros using bilinear transformation.
%
for k = 1:N
   dp(k) = (1 + (Ts/2)*sfs(k))/(1 - (Ts/2)*sfs(k));
   dz(k) = (1 + (Ts/2)*szfs(k))/(1 - (Ts/2)*szfs(k));
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)
%
if rem(N,2) == 1 nd2 = (N+1)/2; else nd2 = N/2; end
if rem(N,2) == 1 & N == 1 d2(1,1)=1; d2(1,2)=-real(dp(1)); d2(1,3)=0; end
if rem(N,2) == 1 & N == 1 n2(1,1)=1; n2(1,2)=1; n2(1,3)=0; end
if rem(N,2) == 1 & N > 1 
	d2(1,1) = 1; d2(1,2) = -dp(1); d2(1,3)=0;
	n2(1,1) = 1; n2(1,2) = -dz(1); n2(1,3)=0;
	for i = 2:nd2
	   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
if rem(N,2) == 0
	for i = 1:nd2
	   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-dp(1)*exp(j*(-1)*(m-1)*step));
end
if N > 1
   for n = 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 at f = 0 is')
disp(max)
b01 = 1/max;
disp('1/max')
disp(b01);
%
%Finding b0
%
if N==1 
    H0 = (1 - dz(1))/(1 - dp(1));
else
    H0 = (1 - dz(1))/(1 - dp(1)); 
    for i = 2:N
        H0 = H0*(1 - dz(i))/(1 - 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:nd2 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
%
if N == 1
    Hqp = -20*log10(abs(b0)*abs(1 - dz(1)*exp(-j*qplpf))/abs(1 - dp(1)*exp(-j*qplpf)));
    Hqs = -20*log10(abs(b0)*abs(1 - dz(1)*exp(-j*qslpf))/abs(1 - dp(1)*exp(-j*qslpf)));    
else
    Hqp =abs(1 - dz(1)*exp(-j*qplpf))/abs(1 - dp(1)*exp(-j*qplpf));
    Hqs =abs(1 - dz(1)*exp(-j*qslpf))/abs(1 - dp(1)*exp(-j*qslpf));
    for p = 2:N
        Hqp =Hqp*abs(1 - dz(p)*exp(-j*qplpf))/abs(1 - dp(p)*exp(-j*qplpf));
        Hqs =Hqs*abs(1 - dz(p)*exp(-j*qslpf))/abs(1 - dp(p)*exp(-j*qslpf));   
    end
    Hqp = -20*log10(abs(b0*Hqp));
    Hqs = -20*log10(abs(b0*Hqs));
end
disp('Loss in dB at the passband cutoff frequency (qplpf)')
disp(Hqp)
disp('Loss in dB at the stopband cutoff frequency (qslpf)')
disp(Hqs)