%
% Butterworth bandstop
%
% Step 1; determine n of prototype lowpass filter
wp1=60; wp2=260; ws1=100; ws2=150; Gp=-2.2; Gs=-20;

wsh = [ ( ws1 * ( wp2-wp1 ) / (  wp1*wp2 - ws1*ws1) );
        ( ws2 * ( wp2-wp1 ) / ( -wp1*wp2 + ws2*ws2) ) ];
wp=1; ws=min(wsh)
P1=-Gs/10; P2=-Gp/10; Wsp=ws;
nc=log((10^P1 - 1)/(10^P2 - 1))/(2*log(Wsp));
n= ceil(nc);
fprintf('Filter order is %i\n', n)

% Step 2 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 and display magnitudes
[num1, den1] = butter(n, Wcp, 's');
[num2, den2] = butter(n, Wcs, 's');

w=logspace(-.3, .3);
[mag1, phase1, w] = bode( num1, den1, w);
[mag2, phase2, w] = bode( num2, den2, w);
semilogx(w, mag1);
hold on
semilogx(w, mag2, 'g');
title('Prototype lowpass filters');
hold off

%  Display prototype filter transfer function
%  Cutoff freq is the one that meets stopband spec.
disp('The prototype lowpass transfer function is')
format short e
num2, den2
pause

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

% Determine two transfer functions and display magnitudes
[num1, den1] = butter(n, Whcp, 'stop', 's');
[num2, den2] = butter(n, Whcs, 'stop', 's');
w=logspace(1, 3);
[mag1, phase1, w] = bode( num1, den1, w);
[mag2, phase2, w] = bode( num2, den2, w);
semilogx(w, mag1);
hold on
semilogx(w, mag2, 'g');
title('Final Bandstop filters');
hold off

%  Now display final filter transfer function
%  Cutoff freq is the one that meets stopband spec.
disp('The final bandpass transfer function is')
format short e
num2, den2
format

% Figure response at stop and pass
[mag1s, phase1, x] = bode(num1, den1, [ws1 ws2]);
[mag1p, phase1, x] = bode(num1, den1, [wp1 wp2]);
mag1s, mag1p
[mag2s, phase2, x] = bode(num2, den2, [ws1 ws2]);
[mag2p, phase2, x] = bode(num2, den2, [wp1 wp2]);
mag2s, mag2p