% m-file file name: delhpf.m
% Function for designing digital elliptic highpass filter.
% 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.
% fphpf = pass band cutoff frequency in Hz.
% fshpf = stop band cutoff frequency in Hz.  Note that fshpf < fphpf
% 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] = delhpf(amax,amin,fphpf,fshpf,fs)
Ts = 1/fs;     % Ts = sampling interval.
qphpf = 2*pi*fphpf*Ts;     % qphpf = digital pass band cutoff frequency in radians.
qshpf = 2*pi*fshpf*Ts;     % qphpf = digital stop band cutoff frequency in radians.
wphpf = 2*fs*tan(qphpf/2);   % wphpf = prewarped pass band cutoff frequency in radians/s.
wshpf = 2*fs*tan(qshpf/2);   % wphpf = prewarped stop band cutoff frequency in radians/s.
%
% Calculating the order of the normalized lowpass filter.
%
k = wshpf/wphpf;
Ws = sqrt(k);
Wp = 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 Highpass Filter Design')
disp(' ')
disp('digital passband cutoff frequency qphpf =')
disp(qphpf)
disp('digital stopband cutoff frequency qshpf =')
disp(qshpf)
disp('prewarped passband cutoff frequency wphpf =')
disp(wphpf)
disp('prewarped stopband cutoff frequency wshpf =')
disp(wshpf)
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
%
sqrtps = sqrt(wphpf/wshpf);
for i = 1:N
    s(i) = sqrtps*s(i);
    sz(i) = sqrtps*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.
%
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 on
title('Normalized Lowpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
%
% Calculating the frequency transformed poles and zeros.
% Lowpass to highpass transformation.
%
wohpf = sqrt(wphpf*wshpf); 
for k = 1:N
   if s(k) == 0 sfs(k) = 0; else sfs(k) = wphpf/s(k); end
   if sz(k) == 0 szfs(k) = 0; else szfs(k) = wphpf/sz(k); end
end
%
% Plotting the frequency transformed poles and zeros.
%
for na = 1:201
   xa(na) = wohpf*(na-101)/100;
   tya(na) = sqrt(wohpf^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')
plot(real(szfs),imag(szfs),'ro')
grid on
title('Frequency Transformed Highpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Frequency Transformed Highpass Pole Locations')
disp(sfs)
disp('Frequency Transformed Highpass Zero Locations')
disp(szfs)
%
% 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 on
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 is ')
disp(max)
b01=1/max;
disp('Constant Term b0 =')
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))/pi;
end
figure
plot(f,mag)
grid on
title('Magnitude Response')
xlabel('Frequency')
ylabel('Magnitude')
figure
plot(f,dB)
grid on
title('Magnitude Response in dB')
xlabel('Frequency')
ylabel('Magnitude in dB')
figure
plot(f,pha)
grid on
title('Phase Response')
xlabel('Frequency')
ylabel('Phase in Radians/pi')
%
% 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 on
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')
%
%Verifying
%
if N == 1
    Hqp = -20*log10(abs(b0)*abs(1 - dz(1)*exp(-j*qphpf))/abs(1 - dp(1)*exp(-j*qphpf)));
    Hqs = -20*log10(abs(b0)*abs(1 - dz(1)*exp(-j*qshpf))/abs(1 - dp(1)*exp(-j*qshpf)));    
else
    Hqp =abs(1 - dz(1)*exp(-j*qphpf))/abs(1 - dp(1)*exp(-j*qphpf));
    Hqs =abs(1 - dz(1)*exp(-j*qshpf))/abs(1 - dp(1)*exp(-j*qshpf));
    for p = 2:N
        Hqp =Hqp*abs(1 - dz(p)*exp(-j*qphpf))/abs(1 - dp(p)*exp(-j*qphpf));
        Hqs =Hqs*abs(1 - dz(p)*exp(-j*qshpf))/abs(1 - dp(p)*exp(-j*qshpf));   
    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)