%FIR_rls:FIR RLS algorithm demo
clear all
close all
hold off
%Number of system points
N=2000;
inp=randn(N,1);
n=randn(N,1);
[b,a]=butter(2,0.25);
Gz=tf(b,a,-1);
%y(n)=b1*u(n)+b2*u(n-1)+b3*u(n-2)-a1*y(n-1)-a2*y(n-2)
%h=[b-a(2:length(a))];
h=[0.0976 0.28732 0.3360 0.2210 0.0976];
%channel system order fixed as we have 5 elements(3 in a and 2 in b)
inporder=3;
outorder=2;
sysorder=inporder+outorder;
y=lsim(Gz,inp);
%add some noise
n=n*std(y)/(15*std(n));
d=y+n;
totallength=size(d,1);
%Take only 50 points for training(N-inporder 47=50-3)
N=50;
%begin of the algorithm
%forgetting factor
lamda=0.999;
%initial P matrix
delta=1e2;
P=delta*eye(sysorder);
w=zeros(sysorder,1);
for n=inporder:N
%u(n),u(n-1),u(n-2)
u=inp(n:-1:n-inporder+1);
%d(n-1),d(n-2)
outp=d(n-1:-1:n-outorder);
u=[u;outp];
phi=u'*P;
k=phi'/(lamda+phi*u);
y(n)=w'*u;
e(n)=d(n)-y(n);
w=w+k*e(n);
P=(P-k*phi)/lamda;
%Just for plotting
Recordedw(1:sysorder,n)=w;
end
%check of results
for n=N+1:totallength
%u(n),u(n-1),u(n-2)
u=inp(n:-1:n-inporder+1);
%d(n-1),d(n-2)
outp=d(n-1:-1:n-outorder);
u=[u;outp];
y(n)=w'*u;
e(n)=d(n)-y(n);
end
hold on
plot(d)
plot(y,'r');
title('系統(tǒng)輸出');
xlabel('抽樣')
ylabel('實際和估計輸出')
figure
semilogy((abs(e)));
title('誤差曲線');
xlabel('抽樣');
ylabel('誤差值');
figure
plot(h,'r+')
hold on
plot(w,'.')
legend('濾波器權(quán)值','估計濾波器權(quán)值',4);
title('濾波器權(quán)值與估計權(quán)值的對比');
figure
plot(Recordedw(1:sysorder,sysorder:N)');
title('估計權(quán)值收斂');
xlabel('抽樣');
ylabel('權(quán)值');
axis([1 N-sysorder min(min(Recordedw(1:sysorder,sysorder:N)'))max(max(Recordedw(1:sysorder,sysorder:N)'))]);