%
% Digital Butterworth Highpass
%
% Oversatisfy stopband magnitude spec
%
%  T = pi/100
%
% Analog specs
wpha=150*pi; wsha=100*pi; Gp=-2; Gs=-10;
T = 1/400;
%
% Step 1: prewarp stop and pass frequencies
%
wph=(2/T)*tan(wpha*T*0.5);
wsh=(2/T)*tan(wsha*T*0.5);
%
% Step 2: determine n of prototype lowpass filter
%
wp=1; ws=wph/wsh;
P1=-Gs/10; P2=-Gp/10; Wsp=ws/wp;
nc=log((10^P1 - 1)/(10^P2 - 1))/(2*log(Wsp));
n= ceil(nc);
fprintf('Filter order is %i\n', n)
%
% Step 3 determine two wc's, one that exactly satisfies passband,
%  and the other satisfies stopband.
Wcp = wp/(10^P2-1)^(1/(2*n));
Wcs = ws/(10^P1-1)^(1/(2*n));
%Wcp, Wcs
%
% Determine two transfer functions, discrete also
[num1, den1] = butter(n, Wcp, 's');
[num2, den2] = butter(n, Wcs, 's');
[numd1, dend1] = bilinear(num1, den1, 1/T);
[numd2, dend2] = bilinear(num2, den2, 1/T);
%
%  Display prototype filter transfer function
%  Cutoff freq is the one that meets passband spec.
disp('The prototype lowpass transfer functions are')
format short e
disp('This satisfies pass band and oversatisfies stop band')
num1, den1
disp('This satisfies stop band and oversatisfies pass band')
num2, den2

%
%  Now design highpass
% first adjust high pass Wp based on cutoff freq of Butterworth
Whcp = wph / Wcp;
Whcs = wph / Wcs;
%
% Determine two highpass transfer functions and display magnitudes
[num1, den1] = butter(n, Whcp, 'high', 's');
[num2, den2] = butter(n, Whcs, 'high', 's');
[numd1, dend1] = bilinear(num1, den1, 1/T);
[numd2, dend2] = bilinear(num2, den2, 1/T);
disp('This satisfies pass band and oversatisfies stop band')
numd1, dend1
disp('This satisfies stop band and oversatisfies pass band')
numd2, dend2
%
%
% Check magnitude at pass and stop freq
% must transform frequencies to discrete
disp('Checking magnitude at pass and stop frequencies:')
% transform discrete frequency to radians around unit circle for freqz()
zs=wsha*T; zp=wpha*T;
mag1=abs(freqz(numd1, dend1, [zp zs]));
mag1
mag2=abs(freqz(numd2, dend2, [zp zs]));
%mag2
%
%  Do comparison between two digital filters
%
figure(1);
%
% 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));
H=freqz(numd2, dend2, z);
H2mag=abs(H);
H2phs=180/pi*unwrap(angle(H));
subplot(211);
plot(wd, H1mag,  'c');
title('Example 12.7, Butterworth design, red meets stop spec');
ylabel('Mag');
xlabel('Discrete freq (rad/sec)');
hold on;
plot(wd, H2mag, 'r');
grid;
hold off;
zoom on;
subplot(212);
plot(wd, H1phs,  'c');
ylabel('Phase (deg)');
hold on;
plot(wd, H2phs, 'r');
grid;
hold off;

%
% SHOULD NOT USE BUTTORD FOR THIS PROBLEM SINCE BUTTORD OVERSATISFIES
% PASSBAND
%
%
%  Now repeat using buttord and butter 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] = buttord( Wpha, Wsha, -Gp, -Gs);
[b,a]=butter(n,Wn,'high');
%
% Display this result
disp('The transfer function using buttord and butter is')
b,a