rls.m

写了著名的最小二乘算法的matlab程序

%--------------------------RLS算法--------------------------%
%-----模型：x(k)+a1*x(k-1)+a2*x(k-2)=w(k)-------------------%
%-----为自适应一步预测器，其中w(k)为方差0.01的高斯白噪声-------%
close all;
clc;clear all;
N=400;
a=0.99;
% 第一组参数
% a1=-0.1950;
% a2=0.95;
% 第二组参数
a1=-1.5955;
a2=0.95;
u1=0.001;
x=zeros(1,N);
w=1*randn(1,N);
wt=zeros(2,N);
M=zeros(2,N);
e=zeros(1,N);
e2=zeros(1,N);
x(1)=w(1);
x(2)=w(2)-a1*x(1);
%----初试值选取
R_inv1=eye(2);
R0=[0,0;0,0];
dxx=0;
pxx=0;
wt(:,3)=[4,-2];
wt1=zeros(2,N);
wt1(:,3)=[4,-2];
e1=zeros(1,N);
for n=3:N
    x(n)=w(n)-a1*x(n-1)-a2*x(n-2);
    xn=[x(n-1);x(n-2)];
        e1(:,n)=x(:,n)-xn'*wt1(:,n);
    wt1(:,n+1)=wt1(:,n)+2*u1*e1(:,n)*xn;
    ee2(:,n)=e1(:,n).^2;
    R1=xn*xn';
    R0=R1+R0;
    dx=x(n)*x(n);
    dxx=dx+dxx;
    px=x(n)*[x(n-1);x(n-2)];
    pxx=px+pxx;
    e(:,n)=x(n)-wt(:,n)'*[x(n-1);x(n-2)];
    e2(:,n)=e(:,n).^2;
    R_inv2=(eye(2)-R_inv1*[x(n-1);x(n-2)]*[x(n-1),x(n-2)]/...
      (a+[x(n-1),x(n-2)]*R_inv1*[x(n-1);x(n-2)]))*R_inv1/a;
    M(:,n)=R_inv2*[x(n-1);x(n-2)];
    wt(:,n+1)=wt(:,n)+e(:,n).*M(:,n);
    R_inv1=R_inv2;
end
wtt=-wt';
R0=R0/(N-1);
% figure(1);plot(wtt)
dxx=dxx/(N-1);
pxx=pxx/(N-1);
wopt=inv(R0)*pxx;
L=length(wtt);
woptL=zeros(N,2);
for i=1:L
    woptL(i,:)=wopt;
end
% figure(1);plot(wtt);
% hold on,plot(-woptL,'r-.');
[w1,w2]=meshgrid(wopt(1)-2:.1:wopt(1)+2,wopt(2)-2:.1:wopt(2)+2);
err=dxx+w1.^2*R0(1,1)+w1.*w2.*(R0(1,2)+R0(2,1))+w2.^2*R0(2,2)-2*w1.*pxx(1)-2*w2.*pxx(2);
figure(2);
contour(w1,w2,err,40);
hold on
plot(wt1(1,3:N),wt1(2,3:N),'k',wt(1,3:N),wt(2,3:N),'r',wopt(1),wopt(2),'m-*');
title('权值收敛曲线');
legend('红色为RLS算法；黑色为LMS算法');
figure(3),subplot(2,1,1),plot(e2);title('RLS学习曲线');
subplot(2,1,2),plot(ee2);title('LMS学习曲线');