RLS_Noise.m

%% RLS_Noise : Recursive Least Square Identification with Noise
% This program objective is to identified parameters of third order transfer
% function. In real-life application usually there's some noise add in
% recorded signal. In this example we use the same transfer funtion, and
% still the same input step, only this time we add some noise.
% Author : Otniel Adrian Rachmat
% Contact information : 
% Department of Electrical Engineering,
% Atma Jaya Catholic University
% Jendral Sudirman 51,
% Jakarta 12930 ,
% Indonesia .
% email :[e-mail]
% email : [e-mail]
% Date    : 09-06-2009
% Version : 1.0
% Reference : 
% (1) System Modeling and Identification, Prentice Hall, Rolf Johansson,1993
% (2) Neuro-Fuzzy and Soft Computing, J.-S.R.Jang, C.-T.Sun, E.Mizutani,1997
% (3) FIR_RLS.m-file, Tamer Abdelazim Mellik, version 1.1.0, 
% Note : The author doesn't take any responsibility for any harm caused by
% the use of this file

%% Generate Transfer Function
% In this example we try to identified some parameters in third order transfer function.
% First we generate some transfer function, using an identified parameters manually.
% Hence we use butter to generate third order lowpass digital Butterworth filter with determine cutoff frequency.

[num,den] = butter(3,0.75)
r_param = [num,den(2:end)];
Gz = tf(num,den,0.01)   
ltiview('step',Gz);

%% Simulate Input-Output Data
% Now that we got the transfer function, we must specified some input data
% to simulate identification process. Normally in identification we record
% a sequence of input-output data, in this term of case we use step
% input.

n_input = 10;       % Number of input data each step
n_step = 10;        % Number of step signal
sample = 4*n_step*n_input;
zero = zeros(n_input,1);
one = ones(n_input,1);
onen = one.*-1;
data = [zero;one;zero;onen];
in_data = [];
for n=1:n_step
    in_data = [in_data;data];
end
out_data = lsim(Gz,in_data);

figure
hold on;
plot(in_data,'b');
plot(out_data,'r');
title('Input-Output Data Identification') ;
xlabel('Samples')
ylabel('Input-Output Signal')

%% Noise addition
% To get a better understanding of the real system we add some noise in the
% output signal we got.

out_sample = lsim(Gz,in_data);
noise = randn(sample,1);
noise = noise * std(out_sample)/(15*std(noise));
out_data = out_sample+noise;

figure
hold on;
plot(out_sample,'b');
plot(out_data,'r');
title('Output Pure-Noise Data Identification') ;
xlabel('Samples')
ylabel('Pure-Noise Signal')

%% RLS Algorithm
% After we obtain input-output data pairs, we try to identified parameters
% of the original transfer function. This can be done by applied Recursive
% Least-Square with input-output data pairs. Vk as Input data at k-times.
% Wk as Output data at k-times.
% W(z)   b0*Z^3 + b1*Z^2 + b2*Z + b3
% ---- = ------------------------------
% V(z)   a0*Z^3 + a1*Z^2 + a2*Z + a3
% Base on the function above, we must identified the parameters such as b0,
% b1, b2, b3, a1, a2, a3.

in_order = 4;
out_order = 3;
sys_order = in_order + out_order;

delta = 1e2;        
Po = delta*eye(sys_order);
Qo = zeros(sys_order,1);
lamda = 0.999;        

for n=4:sample
    for i=1:3
        W(i) = out_data(n-i);
    end
    for i=1:4
        V(i) = in_data(n-i+1);
    end
    reg = [V';W'];
    sample_out = out_data(n);
    est_out(n)= reg' * Qo;
    error(n) = out_data(n) - est_out(n); 
    Po = (Po - ( (Po*reg*reg'*Po)/( lamda+(reg'*Po*reg) ) ))/lamda;
    Qo = Qo + (Po*reg* (sample_out - (reg'*Qo) ) );
    Recordedw(1:sys_order,n) = Qo;
end
f_param = [Qo(1:4);-Qo(5:7)];
numz = [Qo(1:4,1)];
denz = [1;-Qo(5:7,1)];
Fz = tf(numz',denz',0.01);

figure
hold on
plot(r_param,'b*');
plot(f_param,'R+');
title('Parameters system') ;
xlabel('Order')
ylabel('True and estimated output')

figure
plot(Recordedw(1:sys_order,sys_order:sample)');
title('Estimated parameters convergence') ;
xlabel('Samples');
ylabel('Parameter value');
axis([1 (sample)-sys_order min(min(Recordedw(1:sys_order,sys_order:(sample))')) max(max(Recordedw(1:sys_order,sys_order:(sample))')) ]);

figure
semilogy((abs(error))) ;
title('Error curve') ;
xlabel('Samples');
ylabel('Error value');

figure
hold on
plot(out_data);
plot(est_out,'r');
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')

%% Comparison between real and identified transfer function
% In this comparison we found the identified parameters has a distinct
% value to the real one.

Gz
Fz
ltiview('step',Gz,Fz);

%% Conclusion
% Using RLS we may identified parameters of transfer function by means of
% its input-output data. However, when there's noise in the same signal, we
% won't be able to get the exact desire value. But at least we got the
% closest system. Then, is there any way to got the desire ones?