% Program id51.m to illustrate use of first order filters
% Copyright June 1996, B. Joseph Washington University in St.Louis

%
% This program generates a first order filter and uses it to filter a 
% combination of sinewaves
clc;clf;echo on;
%
% Define filter constants
tau= 10 ; % This is filter time constant= 1/wc, wc= cutoff frequency
T=1; %Sampling time used for discretizing the filter
% 
% Define filter transfer function
%
num=[0 1]; den=[10 1];
printsys(num,den);
% Get frequency response
w=logspace(-2,1);% create frequency axis scale 
[mag,phase,w]=bode(num,den,w);
% plot freqency response
loglog(w,mag,'-')
% 
% Now generate the corresponding discrete time filter
%
[numd,dend]=c2dm(num,den,T,'tustin');
[magd,phased,w]=dbode(numd,dend,T,w); % frequency response of discrete filter
hold on; loglog(w,magd,'+');xlabel('frequency rad/sec');
ylabel('Amplitude Ratio continuous(-) and discrete(+)');
 

% see graph to compare with cont version
printsys(numd,dend,'z');
pause
%
% Generate some data to illustrate low and high pass filtering
%
t=0:1:50*pi;ns=length(t);
y=sin(.05*t)+sin(t);
hold off; clf;
% Now apply the filter to the data
% define coefficients
a0=dend(1);a1=dend(2);b0=numd(1);b1=numd(2);
% let yf= filtered data
yf=y;% initilize filtered data;
for i=2:ns; yf(i)=(b0*y(i)+b1*y(i-1)-a1*yf(i-1))/a0;end;
plot(t,y,'-',t,yf,'+');
xlabel('time, sec');ylabel('low pass filtered data')

xlabel('time, sec');ylabel('filtered(+) and unfiltered(-) data')
pause % to see filtered and unfiltered data;
% Now we can try a high pass filter
% This is obtained by merely subtracting the low pass signal
yhp=y-yf;
plot(t,yhp,'-');xlabel('time, sec');ylabel('high pass filtered data')
pause
%
% Now let us construct a second order Butterworth filter 
% and compare its performance on the same data
%
wb=.1; % Desired bandwidth of filter
butter_num=[1];butter_den=[1/wb^2 1.42/wb 1];
printsys(butter_num,butter_den);
[b,a]=c2dm( butter_num,butter_den,T,'tustin');
printsys(b,a,'z');
% 
yf=filter(b,a,y);
plot(t,yf,'-');xlabel('time, sec');ylabel(' 2nd order butterworth  filtered data')