%
% Digital Chebyshev highpass design
%
%
% Analog specs
wpha=150*pi; wsha=100*pi; Gp=-2; Gs=-10;
T = 1/400;
%
% Step 1: prewarp pass and stop frequencies
%
wph=(2/T)*tan(wpha*T*0.5);
wsh=(2/T)*tan(wsha*T*0.5);
%
% Step 2: determine n
%
r=-Gp;
wp=1; ws=wph/wsh;
P1=-Gs/10; P2=r/10; Wsp=ws/wp;
nc=acosh( sqrt( (10^P1-1)/(10^P2-1) ) )/(acosh(Wsp));
n= ceil(nc);
fprintf('Filter order is %i\n', n)
%
% Step 3: determine transfer function of analog prototype low pass
%
[num1, den1] = cheby1(n, r, wp, 's');
%
%w=logspace(1,2);
%[mag1, phase1, w] = bode( num1, den1, w);
%semilogx(w, mag1);
%hold off
%
% Display prototype low pass transfer function
disp('The prototype lowpass transfer function is')
format short e
num1, den1
%
% Step 4 design analog high pass and transform to digital
%
[num1, den1] = cheby1(n, r, wph, 'high', 's');
[numd1, dend1] = bilinear(num1, den1, 1/T);
%
%  Now display final filter transfer function
%
disp('The final transfer function is')
numd1, dend1
format
%
%
% Check magnitude at pass and stop freq
% must transform frequencies to discrete
disp('Checking magnitude at stop and pass frequencies:')
% transform discrete frequency to radians around unit circle for freqz()
zs=wsha*T; zp=wpha*T;
mag1=abs(freqz(numd1, dend1, [zs zp]));
mag1
%
%
% Plot frequency response as a check
% z is the radian frequency to be passed to freqz(), and is the
% product of wd, the discrete frequency and T, the sample period.
% wd, the discrete frequency, is z/T
%
%
z=0:pi/200:pi/2;
wd=z/T;
H=freqz(numd1, dend1, z);
H1mag=abs(H);
H1phs=180/pi*unwrap(angle(H));
subplot(211);
plot(wd, H1mag,  'c');
title('Chebyshev design for highpass');
ylabel('Mag');
xlabel('Discrete freq (rad/sec)');
grid;
zoom on;
subplot(212);
plot(wd, H1phs,  'c');
ylabel('Phase (deg)');
grid;
%
%  Now repeat using cheby1ord and cheby1 for discrete
%  To use these functions, remember that the radian sampling
%  frequency is 2. Must normalize all pass and stop freqs by
%  multiplying by T/pi
Wpha=wpha*(T/pi); Wsha=wsha*(T/pi);
[n, Wn] = cheb1ord( Wpha, Wsha, -Gp, -Gs);
[b,a]=cheby1(n,-Gp, Wn, 'high');
%
% Display this result
disp('The transfer function using cheby1ord and cheby1 is')
b,a