%% Kittiphan Techakittiroj (27 Jan 1996)

%% Input
figure(1);
k = 0:1:100;
x = [ zeros(2,1); [1:50]'; 50*ones(25,1); zeros(24,1)];
w = [ zeros(2,1); 100*ones(50,1); zeros(25,1); zeros(24,1)];
xw = 10e-3*x.*w;
xx = 10e-3*x.*x;
plot(k, xw, k, x, k, xx);
title('an example input');
ylabel('Input (unit)');  xlabel('time (unit)');  grid;
input = [ xw x xx];

%% Model
figure(2);
B0 = [   0   0   2 ];   A0 = [    1 ];  %% delay 0
B1 = [   5   3   0 ];   A1 = [ 0.98 ];  %% delay 1
AA = [ A0 A1];
BB = [ B0 B1];
Model = arx2th( AA, BB, 1, 3);
y = idsim(input, Model);
plot(k, y);
title('an example output from the simulation');
ylabel('Output (unit)');  xlabel('time (unit)');  grid;


%% Add noise to the output
figure(3);
y_n = y + max(y)*0.1*(randn(size(y)) - 0.5);
plot(k, y_n);
title('an example output from the simulation + noise');
ylabel('Output (unit)');  xlabel('time (unit)');  grid;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ARX Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Identification
ident_data = [y input];
na = [1];
nb = [1 1 1];
nc = [1 1 0];
[I_A I_B] = th2arx(arx(ident_data, [na nb nc]));

%% Compare
I_A
I_B

%% Identification
ident_data = [y_n input];
na = [1];
nb = [1 1 1];
nc = [1 1 0];
[I_A I_B] = th2arx(arx(ident_data, [na nb nc]));

%% Compare
I_A
I_B


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RARX Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lam = 0.9;
%% Shift the input because arx doesn't except nc = 0
input = [ xw x [xx(2:length(xx));0]];

%% Identification
ident_data = [y input];
na = [1];
nb = [1 1 1];
nc = [1 1 1];
I_Model = rarx(ident_data, [na nb nc], 'ff', lam);

%% Compare
I_Model(length(ident_data),:)

%% Identification
ident_data = [y_n input];
na = [1];
nb = [1 1 1];
nc = [1 1 1];
I_Model = rarx(ident_data, [na nb nc], 'ff', lam);

%% Compare
I_Model(length(ident_data),:)
