% EML 6934, Fall 09, Prabir Barooah
% define a continuous-time plant, discretize it, 
% convert it to ARX form and simulate 


clear all

% ---- continuous time to discrete time -----
Hc = tf(100,[1 0.2 100]);
sampling_period = 0.01;
Hd = c2d(Hc,sampling_period,'tustin');
[alpha_vec,beta_vec] = tfdata(Hd,'v');
alpha_rowvec = [alpha_vec(:)]';
tmp = beta_vec(2:end);
beta_rowvec = [tmp(:)]';

%%-------- choose simulation parameters -----
N = 1000;
time_vec = [0:1:N-1]'*sampling_period;

%choose input signal
u_vec = 0.25*sin(2*pi*6*time_vec);%6 Hz sinusoid

%intialize data
y = zeros(N,1);

% choose initial conditions (I choose 0 but anything else is also OK)
y(1) = 0;
y(2) = 0;

% run simulation
for ii=3:N
    y_rhs = [y(ii-1) y(ii-2)]';
    u_rhs = [u_vec(ii) u_vec(ii-1) u_vec(ii-2)]';
    y(ii,1)= -beta_rowvec*y_rhs + alpha_rowvec*u_rhs;    
end;


figure;
plot(time_vec,y,'r.-',time_vec,u_vec,'b.-');

%compare with MATLAB simulation with 0 initial condition
y_m = filter(alpha_vec,beta_vec,u_vec);
figure
plot(time_vec,y,'r.-',time_vec,y_m,'b.');
%I can;t explain the difference between the two.
%whoever can figure out the reason gets a prize.