%
%  Digital Butterworth lowpass, for both wc's
%
% Analog specs
wpa=750; wsa=1500; Gp=-1; Gs=-20;

T = 1/(2*pi*900);

% Step 1: prewarp stop and pass frequencies
wp=(2/T)*tan(wpa*T*0.5);
ws=(2/T)*tan(wsa*T*0.5);

% Step 2: determine n
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));

% Step 4 Determine two transfer functions and discrete versions. 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 final filter transfer functions (analog and digital)
format short e
disp('The final transfer function is (meeting passband spec.)')
num1, den1
disp('----------DIGITAL---------------');
numd1, dend1
disp('=-----------------------------------------------------=');
disp('The final transfer function is (meeting stopband spec.)')
num2, den2
disp('----------DIGITAL---------------');
numd2, dend2
%

%  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
%Wpa=wpa*(T/pi); Wsa=wsa*(T/pi);
%[n, Wn] = buttord( Wpa, Wsa, -Gp, -Gs);
%[b,a]=butter(n,Wn);
%
% Display this result
%disp('The transfer function using buttord and butter is')
%b,a