%   Butterworth filter Download and run

%   H(s) = 1/(s/w1 + (1-j)/sqrt(2))(s/w1 + (1+j)/sqrt(2)) = w1^2/(s-s1)(s-s2)
%   H(z) = (1-z1)(1-z2)/(z-z1)(z-z2)
%   where z=exp(jwT) and z1 = exp(s1T) , z2 = exp(s2T), s1 = (1-j)*w1/sqrt(2), s2 = (1+j)*w1/sqrt(2)

Fs=1e4	%	10 kHz sampling frequency
T=1/Fs;
N=256;          %   N samples
w1 = 2*pi*1e3;  %   Filter corner frequency = 1 kHz

freq = 0:Fs/N:Fs;       %   frequencies are multiples of the fundamental frequency Ff = Fs/N
w = 2*pi*freq;
z=exp(j*w*T);           %   Z-Transform variable 

s1 = -w1*(1-j)/sqrt(2);   %   First pole
s2 = -w1*(1+j)/sqrt(2);   %   Second pole

hs = w1^2./((j*w-s1).*(j*w-s2));     %   Continuous filter transfer function
hsdB=20*log10(abs(hs));

z1=exp(s1*T);
z2=exp(s2*T);
%   hz = (1-z1)*(1-z2)/((z-z1).*(z-z2));    %   Discrete filter transfer function
%   Since z1 = conj(z2)   (1-z1)*(1-z2) = 1 -2rept(z1) - z1*conj(z1)
%   and (z-z1)*(z-z2) = z^2 -(z1+z2)*z + z1*z2 = z^2 - 2rept(z1)*z + z1*z2
a=(1-z1)*(1-z2);
b= z1+z2;
c=z1*z2;
hz = a./(z.^2-b*z + c);
hzdB=20*log10(abs(hz));

%   Set up a multitone signal to test the filter.
testbins =[2 4 10 20 22 24 25 26 27 28 30 40 60 120]; %   14 tones with 8 near w1
x=zeros(1,N);         %   Initilize x at zero
y=x;
for k=1:length(testbins)
    for n=1:N
    x(n) = x(n) + cos(pi/2 + 2*pi*testbins(k)*(n-1)/N); %   Tones have 90 degree phase so x(1)=0
    end
end


y(1)=0;
y(2)=0;

for n=1:N-2
    y(n+2) = b*y(n+1)-c*y(n)+a*x(n);   %   Software filter
end
L=N/2;


yf=2*fft(y)/N;      %   Take the Fourier transform of y(n)
yf(1)=yf(1);                %   Correct the DC term
ym = abs(yf);               %   Get amplitude
ydB=20*log10(ym);           %   Convert to dB

h=semilogx(freq(1:L),hsdB(1:L), '-',freq(1:L),hzdB(1:L),':',freq(testbins+1),ydB(testbins+1),'o')
ylabel('Filter gain (dB)'), xlabel('frequency (Hz)')
legend(h,'Continuous Filter', 'Discrete Filter','Software Filter')


