% signal generation
% y = 1.2 f + e
close all
clear all
N=input('number of data points?');
f=randn(N,1);y=1.2*f+.1*randn(size(f));
R=0.01;
% initializing the RLS scheme
x(1)=0; P(1) = .001;
for i=2:N,
    K(i)=P(i-1)*f(i-1)*inv(f(i-1)*P(i-1)*f(i-1)' + R);
    x(i)=x(i-1)+K(i)*(y(i-1)-f(i-1)*x(i-1));
    P(i)=P(i-1)-P(i-1)*f(i-1)'*K(i)';
end
figure(1),plot(P),title('Covariance par estimate'),
figure(2),plot(x),title('True parameter is 1.2');
pause,
f=randn(N,1);y=-1.2*f+.1*randn(size(f));
P(1)=P(N);x(1)=x(N);
% rescue P(1)=P(N)+0.001;x(1)=x(N);
xf = x';
Pf = P';
for i=2:N,
    K(i)=P(i-1)*f(i-1)*inv(f(i-1)*P(i-1)*f(i-1)' + R);
    x(i)=x(i-1)+K(i)*(y(i-1)-f(i-1)*x(i-1));
    P(i)=P(i-1)-P(i-1)*f(i-1)'*K(i)';
end
figure(1),plot(P),title('Covariance par estimate');
figure(2),plot(x),title('True parameter is -1.2');
xf = [xf;x'];
Pf = [Pf;P'];
pause,
figure(1),plot(Pf),title('Covariance evolution over full time window');
figure(2),plot(xf),title('True parameter switch from 1.2 to -1.2');
