%
% Digital Butterworth bandpass
% satisfy pass band spec
%
% Analog specs
wp1a=120; wp2a=300; ws1a=45; ws2a=450; Gp=-2; Gs=-12;
T = 1/1000;

% Step 1: prewarp stop and pass frequencies
wp1=(2/T)*tan(wp1a*T*0.5);
wp2=(2/T)*tan(wp2a*T*0.5);
ws1=(2/T)*tan(ws1a*T*0.5);
ws2=(2/T)*tan(ws2a*T*0.5);

% Step 2: determine n of prototype lowpass filter
wsh = [ (  wp1*wp2 - ws1*ws1) / ( ws1 * ( wp2-wp1 ) );
        ( -wp1*wp2 + ws2*ws2) / ( ws2 * ( wp2-wp1 ) ) ];
wp=1; ws=min(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 prototype LP transfer functions
[num1, den1] = butter(n, Wcp, 's');
[num2, den2] = butter(n, Wcs, 's');

%  Display prototype filter transfer function
%  Cutoff freq is the one that meets passband spec.
disp('The prototype analog lowpass transfer function is')
format short e
num1, den1
pause

%  Now design bandpass
% first adjust Wp's based on cutoff freq of Butterworth
%  Conversion from buttord() function
W0 = [-Wcp Wcp];
Whcp = -W0*(wp2-wp1)/2 + sqrt( W0.^2/4*(wp2-wp1)^2 + wp1*wp2);
Whcp = sort(abs(Whcp));
W0 = [-Wcs Wcs];
Whcs = -W0*(wp2-wp1)/2 + sqrt( W0.^2/4*(wp2-wp1)^2 + wp1*wp2);
Whcs = sort(abs(Whcs));

% Determine two transfer functions and display magnitudes
[num1, den1] = butter(n, Whcp, 's');
[num2, den2] = butter(n, Whcs, '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
%

