%	First Order Lowpass Filter

%	Download and Run

Fs=10e3;        %   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 

hs = 1./(1+j*w/w1);     %   Continuous filter transfer function
hsdB=20*log10(abs(hs));

hz = (1-exp(-w1*T))./(1-exp(-w1*T)*z.^(-1));    %   Discrete filter transfer function
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

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;

for n=1:N-1
    y(n+1) = (1-exp(-w1*T))*x(n+1) + exp(-w1*T)*y(n);   %   Software filter
end

L=N/2;
lengthy=length(y);

yf=2*fft(y)/length(y);      %   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
ydBtest = ydB(testbins+1)
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')


