%
% Chebyshev bandpass
%
% Step 1; determine n of prototype lowpass filter

wp1=1000; wp2=2000; ws1=450; ws2=4000; Gp=-1; Gs=-20;
r=-Gp;
wsh = [ (  wp1*wp2 - ws1*ws1) / ( ws1 * ( wp2-wp1 ) );
        ( -wp1*wp2 + ws2*ws2) / ( ws2 * ( wp2-wp1 ) ) ];

wp=1; ws=min(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 2 determine transfer function of prototype lowpass

% Determine two transfer functions and display magnitudes
[num1, den1] = cheby1(n, r, 1, 's');

w=logspace(-.3, .3);
[mag1, phase1, w] = bode( num1, den1, w);
semilogx(w, mag1);
title('Prototype lowpass filter');
hold off
%  Display prototype filter transfer function

disp('The prototype lowpass transfer function is')
format short e
num1, den1
pause

%  Now design bandpass
Whcp = [wp1 wp2];
[num1, den1] = cheby1(n, r, Whcp, 's');

w=logspace(1.0, 3.0);
[mag1, phase1, w] = bode( num1, den1, w);
semilogx(w, mag1);
title('Final Bandpass filter');
hold off

%  Now display final filter transfer function

disp('The final Bandpass transfer function is')
format short e
num1, den1
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
