LMS算法（在线性预测器权重估计上的应用）

LMS算法原理

最陡下降算法，可用迭代的方法求解维纳-霍夫方程，但是得知道信号相关矩阵R和互相关向量p，当两者确定后，其实迭代的过程和结果都是确定的。LMS算法就是利用输入信号和期望信号都是平稳过程，用有限观测样本的时间平均来逼近统计平均：

                                                                          

                                                                        

然后同最陡下降算法的迭代公式，化简即可得到权重更新公式：   

• 首先产生信号样本序列

clc,clear all,close all
%% 产生样本序列
%由题目得差分方程为 u(n)=0.975u(n-1)-0.95u(n-2)+v(n)
% 写出差分方程初值条件以及 num 和 den
u0=[0 0];
den=[1 -0.975 0.95]; %差分左边
num=[1]; %差分右边
vsigma=0.0731;number=512;
v=normrnd(0,0.0731,number,1);
% v=sqrt(vsigma)*rand(number,1);
% 生成初始条件
Zi=filtic(num,den,u0); %v(n)每一次已经固定了，所以初始条件没有
u=filter(num,den,v,Zi);%这里的 v 的信息需包括样本序列个数
% 绘制样本序列
plot_n=1:number;
figure(1),plot(plot_n,v,'R-',plot_n,u,'b--');
xlabel('n');ylabel('y(n)')
legend('加性白噪声 v ','样本序列 u(n)');

                        

• 利用LMS和上述产生的512个样本序列进行二阶线性预测器权重估计：（用了两种步长μ）

%% 两种步长下 LMS 估计 w1 和 w2
clc,clear all,close all
u0=[0 0];
den=[1 -0.975 0.95]; %差分左边
num=[1]; %差分右边
vsigma=0.0731;number=512;
v=normrnd(0,0.0731,number,1);
% v=sqrt(vsigma)*rand(number,1);
% 生成初始条件
Zi=filtic(num,den,u0); %v(n)每一次已经固定了，所以初始条件没有
u=filter(num,den,v,Zi);
mu1=0.05;mu2=0.005;
w1=zeros(2,number);% mu1 情况下 2 阶 LP 的权系数
w2=zeros(2,number);% mu2 情况下
e1=zeros(number,1);% mu1 情况下的估计误差
e2=zeros(number,1);
d1=zeros(number,1); %mu1 情况下的估计输出
d2=zeros(number,1);
% 为了简便，取 2 阶抽头都有值的时候开始迭代
for n=3:number-1
w1(:,n+1)=w1(:,n)+mu1*u(n-1:-1:n-2,:)*conj(e1(n));
w2(:,n+1)=w2(:,n)+mu2*u(n-1:-1:n-2,:)*conj(e2(n));
d1(n+1)=w1(:,n+1)'*u(n:-1:n-1,:);
d2(n+1)=w2(:,n+1)'*u(n:-1:n-1,:);
e1(n+1)=u(n+1,:)-d1(n+1);%线性预测期 LP 的输出期望值就为 u(n)
e2(n+1)=u(n+1,:)-d2(n+1);
end
%绘制步长为 mu1 下的权系数迭代
figure(2),plot(plot_n,w1(1,:),'R-',plot_n,w1(2,:),'b--');
xlabel('n');ylabel('w(n)')
legend('w1','w2');
%绘制步长为 mu2 下的权系数迭代
figure(3),plot(plot_n,w2(1,:),'R-',plot_n,w2(2,:),'b--');
xlabel('n');ylabel('w(n)')
legend('w1','w2');

                       

• LMS会有一个剩余均方误差和失调参数来描述其均方误差的稳态统计性能：

（两种步长情况的权重收敛曲线和学习曲线比较）

vnew=normrnd(0,0.0731,number,1,100); %产生 512*1*100 的噪声
Zinew=filtic(num,den,u0);
unew=filter(num,den,vnew,Zinew);
neww1=zeros(2,number,100);% mu1 情况下 2 阶 LP 的权系数
neww2=zeros(2,number,100);% mu2 情况下
newe1=zeros(number,1,100);% mu1 情况下的估计误差
newe2=zeros(number,1,100);
newd1=zeros(number,1,100); %mu1 情况下的估计输出
newd2=zeros(number,1,100);
% J1=zeros(number,100);
% J2=zeros(number,100);
for n=3:number-1
for m=1:100
neww1(:,n+1,m)=neww1(:,n,m)+mu1*unew(n-1:-1:n-2,:,m)*conj(newe1(n,m));
neww2(:,n+1,m)=neww2(:,n,m)+mu2*unew(n-1:-1:n-2,:,m)*conj(newe2(n,m));
newd1(n+1,m)=neww1(:,n+1,m)'*unew(n:-1:n-1,:,m);
newd2(n+1,m)=neww2(:,n+1,m)'*unew(n:-1:n-1,:,m);
newe1(n+1,m)=unew(n+1,:,m)-newd1(n+1,m);%线性预测期 LP 的输出期望值就为
u(n)
newe2(n+1,m)=unew(n+1,:,m)-newd2(n+1,m);
rm=xcorr(unew(:,:,m),'biased');
R=[rm(512),rm(513);rm(511),rm(512)];
p=[rm(511);rm(510)];end
end
newe1=newe1.^2;
newe2=newe2.^2;
for i=1:2
for n=1:512
averagee1(n,1)=mean(newe1(n,1,:));
averagee2(n,1)=mean(newe2(n,1,:));
averagew1(i,n)=mean(neww1(i,n,:));
averagew2(i,n)=mean(neww2(i,n,:));
end
end
figure(4),plot(plot_n,w1(1,:),'b-',plot_n,w1(2,:),'b--
',plot_n,averagew1(1,:),'r--',plot_n,averagew1(2,:),'r--');
xlabel('n');ylabel('w(n)')
figure(5),plot(plot_n,w2(1,:),'b-',plot_n,w2(2,:),'b--
',plot_n,averagew2(1,:),'r--',plot_n,averagew2(2,:),'r--');
xlabel('n');ylabel('w(n)')
figure(6),plot(plot_n,averagee1,'b-',plot_n,averagee2,'r--');
xlabel('n');ylabel('MSE')
legend('u1','u2');
%% 计算剩余误差和失调参数
wopt=zeros(2,100); %理论求解最佳权值
Jmin=zeros(1,100);
sum_eig=zeros(100,1);
for m=1:100
%计算每次实验的自相关函数
rm=xcorr(unew(:,:,m),'biased');
%LP(2)，则自相关矩阵为[r(0),r(-1);r(1),r(0)]
R=[rm(512),rm(513);rm(511),rm(512)];
p=[rm(511);rm(510)];
wopt(:,m)=R\p; %维纳-霍夫方程求 R-1p
[v,d]=eig(R);
Jmin(m)=rm(512)-p'*wopt(:,m);%LP 模型
sum_eig(m)=d(1,1)+d(2,2); %特征值之和
end
%平均误差
meanJmin=sum(Jmin)/100;
sum_eig=sum(sum_eig)/100;
Jex1=mu1*meanJmin*sum_eig/(2-mu1*sum_eig); %计算剩余均方误差
Jex2=mu2*meanJmin*sum_eig/(2-mu1*sum_eig);
M1=Jex1/meanJmin;
M2=Jex2/meanJmin;