%
% Butterworth Low pass
%
% Step 1; determine n of prototype lowpass filter

wp=10; ws=20; Gp=-2; Gs=-20;

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 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);
title('Prototype lowpass filter');
hold off

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