%Simple matlab script to design an IIR low pass filter using
%butterworth in either impulse inv. or bilinear method.

%
% EE 420, CSU Fullerton
%
% by Nasser M. Abbasi
% 5/5/201

% Filter SPECIFICATIONS
clear all;
close all;
Fs=10000;
fp=1000;
fs=2000;
dbp=-3;
dbs=-10;
BILINEAR=1;  %set this to 0 to do impulse inv.

%%%%%%%%%%%%%%%%%%%%

fprintf('***** STARTING DESIGN *******\n');
fprintf('Sampling frequency=%f Hz\n',Fs);
fprintf('freq at passband=%f Hz\n',fp);
fprintf('freq at stopband=%f Hz\n',fs);
fprintf('db at passband=%f \n',dbp);
fprintf('db at stopband=%f \n',dbs);
if BILINEAR
   fprintf('Doing Bilinear method\n');
else
   fprintf('Doing impulse invariance method\n');
end



if BILINEAR
    T=1/Fs;
else
    T=1;
end

fprintf('T=%f\n',T);

wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
alphap=1/(10^(dbp/10));
alphas=1/(10^(dbs/10));

if BILINEAR
    gammap=2/T * tan(wp/2);
else
    gammap=wp/T;
end

if BILINEAR
    gammas=2/T * tan(ws/2);
else
    gammas=ws/T;
end

oldn=.5*(log10(alphap-1)-log10(alphas-1))/(log10(gammap)-log10(gammas));
n=ceil(oldn);
fprintf('n=%f, rounded to %d\n',oldn,n);

if BILINEAR
    gammaC=gammas/(10^(1/(2*n)*log10(alphas-1)));
else
    gammaC=gammap/(10^(1/(2*n)*log10(alphap-1)));
end

fprintf('Gamma C =%f\n',gammaC);

poles_of_hs=zeros(n,1);
for i=0:n-1
    poles_of_hs(i+1)=gammaC*exp(sqrt(-1)*(pi*(1+2*i+n)/(2*n)));
end

fprintf('POLES Of H(s)\n');
poles_of_hs

k=prod(-poles_of_hs);
den=poly(poles_of_hs);

fprintf('k=%d\n',k);

fprintf('H(s)=\n',k);


if BILINEAR
    hs=tf(k,den)
else
    hs=tf(T*k,den)    
end


[r,p,k]=residue(k,den);
hzp=zeros(n,1);
for i=1:n
    hzp(i)=exp(p(i)*T);
end

fprintf('POLES Of H(z)\n');
hzp

[B,A]=residue(r,hzp);
fprintf('H(z)=\n',k);
hz=tf(T*B',A')







