% m-file file name: dichpf.m
% Function for designing digital inverse Chevyshev highpass filter.
% Specifications are given in terms of amax, amin, fplpf, fslpf.
% 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 passband.
% amin = minimum loss in dB in the stopband.
% fphpf = passband cutoff frequency in Hz.
% fshpf = stopband 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:
% 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] = dichpf(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.
N = ceil(acosh(sqrt((10^(0.1*amin)-1)/(10^(0.1*amax)-1)))/acosh(wphpf/wshpf)); % N = order of the filter.
disp('Digital Inverse Chebyshev (Chebyshev Type II) 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.
%
a = (1/N)*asinh(sqrt(10^(0.1*amin)-1));
if rem(N,2) == 1 s(1) = -1/sinh(a); end
if N == 1 max = 0; 
   elseif rem(N,2) == 1 max = (N-3)/2;
      else max = (N-2)/2;
   end
end
if N > 1
 for m = 0:max
   if rem(N,2) == 1
      s(2*m+2) = -1/[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) = -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
%
% 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 p = 1:r
         sz(2*p) = j*sec(pi*(2*p-1)/(2*N));
         sz(2*p+1) = conj(sz(2*p));
      end
end
if rem(N,2) == 0
      for p = 1:r
         sz(2*p-1) = j*sec(pi*(2*p-1)/(2*N));
         sz(2*p) = conj(sz(2*p-1));
      end
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')
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
grid on
title('Normalized Lowpass Poles and Zeros')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Normalized Lowpass Pole Locations')
disp(s)
disp('Normalized Lowpass Zero Locations')
for i = 1:N
	disp(sz(i))
end
%
% Calculating the frequency transformed poles and zeros.
% Lowpass to highpass transformation.
%
wohpf = wshpf;  % wohpf is the frequency transformed stopband cutoff frequecy.
for k = 1:N
   if s(k) == 0 sfs(k) = 0; else sfs(k) = wohpf/s(k); end
   if sz(k) == 0 szfs(k) = 0; else szfs(k) = wohpf/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;
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;
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(M/2 + 1)); b0 = 1/max;
disp('maximum gain at f = 1/2 is')
disp(max)
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')