% Matlab code for simulation of x(k+1) = ad x(k) + bd u(k)
% Prabir Barooah

clear all;

% parameters of the continuous time model
a = -0.2;
b = 10;

T = 1e-2; %sampling period, in second;

forcing_freq = 8; %in Hz

% parameters of the discrete time model
ad = exp(a*T);
bd = b*(ad-1)/a;

%initial condition
xinit = 0.5;


N = round(1/T); %number of data points collected in a single experiment
time_vec = [0:1:N-1]'*T;

std_dev_u = 0.01;
std_dev_x = 0.02;

u_vec = sin(2*pi*forcing_freq*time_vec); %input - unit amplitude sinusoid


%z = zeros(N,1);
%H = zeros(N,2);

%intialize
x_prev = xinit;
for ii=1:N

    u_actual = u_vec(ii) + std_dev_u*randn(1,1);
    x_next = ad * x_prev + bd*u_actual;

    % update
    x_prev = x_next;

end;
% once you construct H and z,
% use H\z in Matlab to compute the least squares solution of Hx = z
% theta_hat = H\z;