% m-file file name: dchhpf.m
% Function for designing digital Chevyshev highpass filter.
% Specifications are givem in terms of amax, amin, fphpf, fshpf.
% 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] = dchhpf(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 Chebyshev (Type I) 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(1/sqrt(10^(0.1*amax)-1));
disp('The a value is:')
disp(a)
if rem(N,2) == 1 s(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) = -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) = -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
%
% 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');grid;
title('Normalized Lowpass Pole Locations')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Normalized Lowpass Pole Locations')
disp(s)
%
% Calculating the normalized lowpass zero locations.  
% Vector szb is the normalized zeros.
%
for k = 1:N
   sz(k) = 10^50;
end
disp('Normalized Lowpass Zero Locations')
disp(sz)
%
% Calculating the frequency transformed poles and zeros.
% Lowpass to highpass transformation.
%
wohpf = wphpf;  % wohpf is the frequency transformed half-power 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(1));
for m = 2:M
   if abs(H(m)) > max max = abs(H(m)); end
end
disp('maximum gain is')
disp(max)
b0=1/max;
disp('constant factor 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;
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/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;
title('Impulse Response')
xlabel('Time')
ylabel('Amplitude')