adsp华科实验报告材料Lms算法.docx
- 文档编号:25459544
- 上传时间:2023-06-09
- 格式:DOCX
- 页数:19
- 大小:129.61KB
adsp华科实验报告材料Lms算法.docx
《adsp华科实验报告材料Lms算法.docx》由会员分享,可在线阅读,更多相关《adsp华科实验报告材料Lms算法.docx(19页珍藏版)》请在冰豆网上搜索。
adsp华科实验报告材料Lms算法
现代数字信号处理实验报告
院系:
班级:
姓名/学号:
题目一
题目:
按照课本第三章63页的要求,仿真实现LMS算法和RLS算法,比较两种算法的权值收敛速度,并对比不同λ值对RLS算法的影响。
解答:
1.1数据模型
(1)均值为0、方差为1的白噪声信号:
用randn函数产生均值为0、方差为1的标准正态分布随机矩阵来实现。
根据a1=-1.6,a2=0.8,x(n)+a1*x(n-1)+a2*x(n-2)=e(n);来生成正真的x(n)序列,其中e(n)表示的是噪声信号。
(2)信号点数这里取为1000,用1000个信号来估计滤波器系数。
分别用LMS,RLS算法来实现,并通过信号点来反映自适应权系数的过度过程,以此来比较。
(3)分别取4个不同的λ值来分析RLS算法对收敛效果的影响。
其中lambda=[1,0.98,0.96,0.94]。
1.2算法模型
1.2.1LMS算法模型
fori=3:
1:
N%权系数迭代N次
W=W+2*u*X*e(i-1);%LMS算法的权系数迭代公式
X=[x(i-1)x(i-2)]';%LMS算法中输入信号矢量的递推
e(i)=x(i)-W'*X;
a1L(i)=-W
(1);%LMS算法中权系数a1的提取
end
LMS算法最核心的思想是用平方误差代替均方误差,每次迭代中使用很粗略的梯度估计值来代替精确的梯度值,因而权系数的调整过程是有噪声的,w不再是确定性函数而变成了随机变量,当迭代过程收敛后,权矢量将在最佳权矢量附近随机起伏。
1.2.2RLS算法模型
fori=3:
1:
N
X=[x(i-1),x(i-2)]';
Rxx=lbda*Rxx+X*X';%迭代公式中自相关矩阵的计算
e=x(i)-W'*X;%输出信号误差e(n\n-1)
W=W+inv(Rxx)*X*e;%RLS算法的权系数迭代公式
a1R(i)=-W
(1);%LMS算法中权系数a1的提取
end;
RLS算法的关键是用二乘方的时间平均的最小化准则取代最小均方准则,并按时间进行迭代计算。
对于非平稳随机输入信号引入一个加权因子对i时刻的误差值进行修正,由于加权因子取值范围从0~1,因此时间越久远的数据对误差的贡献就越小。
其最大的特点是收敛速度快,但每次的迭代计算量大。
1.3仿真过程简介
仿真过程按照如下过程进行
(1)信号产生:
首先产生高斯白噪声序列Noise(n),Noise=randn(1,n);然后通过一个简单的二阶自回归滤波器生成信号
,该滤波器a1=-1.6,a2=0.8。
头两个可以表达出来x
(1)=Noice
(1);x
(2)=Noise
(2)-a2*x
(1);由x(n)=Noise(n)-a1*x(n-1)-a2*x(n-2);从而生成序列。
(2)将步骤一生成的信号通过LMS和RLS自适应滤波器进行处理
(3)通过改变λ值对
收敛速度的影响来分析RLS算法的性能。
(4)绘制各种图形曲线
(5)源代码如下:
按照步骤一,生成信号的程序如下所示,生成的信号包含n个点
(1)%输入信号序列的产生
clearall;
clc;
%初始化参数
a1=-1.6;%生成信号x(n)的参数
a2=0.8;
N=1000;%信号点数
%信号及白噪声信号序列的初始化
x=zeros(1,N)';%信号的初始化
Noise=randn(1,N)';%白噪声的初始化,均值为0,方差为1
x
(1)=Noise
(1);%信号前两点的初始赋值
x
(2)=Noise
(2)-a1*x
(1);
%信号序列的产生
forn=3:
N
x(n)=Noise(n)-a1*x(n-1)-a2*x(n-2);
end
生成信号之后,将信号输入LMS和RLS滤波器进行处理,其matlab程序如下:
(2)%LMS和RLS算法下的参数a1的收敛曲线
%LMS滤波
u=0.002;%LMS算法下自适应增益常数初始化
Rxx=[x
(1)*x
(1)0;00];
T=0;
e
(1)=0;
W=[0;1];
X=[x
(2);x
(1)];
e
(2)=x
(2)-W'*X;
fori=3:
1:
N%权系数迭代N次
Rxx=[x(i-1)*x(i-1)x(i-1)*x(i-2);x(i-1)*x(i-2)x(i-2)*x(i-2)];%列出自相关矩阵
T=1/(T+trace(Rxx));%求出迹的值,为后续u的判断做准备
if(u>T)
error('uislargerthan1/t[R]');%判断u的值是否小于迹的倒数
end
W=W+2*u*X*e(i-1);%LMS算法的权系数迭代公式
X=[x(i-1)x(i-2)]';%LMS算法中输入信号矢量的递推
e(i)=x(i)-W'*X;
a1L(i)=-W
(1);%LMS算法中权系数a1的提取
end
%RLS滤波
lbda=1;%RLS算法下lambda取值
T=eye(2,2)*10;%%RLS算法下T参数的初始化,T初始值为10
Rxx=inv(T);
w1
(1)=0;w2
(1)=0;
w1
(2)=0;w2
(2)=0;
W=[w1
(2);w2
(2)];%初始权值
P=[0;0];
fori=3:
1:
N
X=[x(i-1),x(i-2)]';
Rxx=lbda*Rxx+X*X';%迭代公式中自相关矩阵的计算
e=x(i)-W'*X;%输出信号误差e(n\n-1)
W=W+inv(Rxx)*X*e;%RLS算法的权系数迭代公式
a1R(i)=-W
(1);%LMS算法中权系数a1的提取
end;
%绘制LMS与RLS算法下a1、a2收敛曲线
figure
(1)
plot(1:
N,a1L,'r-',1:
N,a1R,'b-',1:
N,a1,'k-');
legend('LMS-a1变化','RLS-a1变化','a1收敛值',0);%图例
title('自适应权系数a1(n)的过渡过程(RLS和LMS算法比较)');
xlabel('信号点数N');
ylabel('对应a1的值');
(3)%RLS算法下不同lambda值的参数收敛曲线
w1
(1)=0;w2
(1)=0;
w1
(2)=0;w2
(2)=0;
W=[w1
(2);w2
(2)];
lambda=[1,0.98,0.96,0.94];
forj=1:
4
fori=3:
N
X=[x(i-1),x(i-2)]';
Rxx=lambda(j)*Rxx+X*X';%迭代公式中自相关矩阵的计算
e=x(i)-W'*X;%输出信号误差e(n\n-1)
W=W+inv(Rxx)*X*e;%RLS算法的权系数迭代公式
a1R1(i)=-W
(1);
a1R2(i)=-W
(1);
a1R3(i)=-W
(1);
a1R4(i)=-W
(1);
end
end
figure
(2)
plot(1:
N,a1R1,'b-',1:
N,a1R2,'r:
',1:
N,a1R3,'c-.',1:
N,a1R4,'g-.',1:
N,a1,'k')%画图显示不同lamda值下RLS算法性能差别
legend('lambda=1','lambda=0.98','lambda=0.96','lambda=0.94','a1收敛值',0)%图例
title('RLS算法下遗忘因子对权系数a1收敛速度影响')
xlabel('信号点数N')
ylabel('对应a1的值')
1.4仿真结果及分析
LMS和RLS算法下的参数a1的收敛曲线如下图:
图1自适应权系数a1(n)的过渡过程(RLS和LMS算法比较)
上图是两种算法下a1的收敛曲线,其中黑色的线表示a1原值,红色的线表示LMS算法关于a1的收敛曲线,蓝色的表示RLS的收敛曲线。
由可以看到在u和
值相同的情况下(u=0.002,
=1),RLS比LMS具有更快的收敛速度,但是RLS权系数收敛后出现了较大的噪声。
由于LMS算法是以梯度增量逐渐逼近最佳值,因此较为平缓的变化趋势;而RLS算法曲线起步变化幅度比较大,但其能够迅速找到最佳值,权系数收敛后出现了较大的噪声因为RLS中有效记忆长度只有49,在计算最加权系数时没有充分利用能够获取的全部采样数据。
图2RLS算法下遗忘因子对权系数a1(n)的影响
图中分别画出来关于四种不同
情况下,RLS算法下遗忘因子对权系数a1(n)的影响可知
越大,RLS算法的收敛速度越慢,算法的有效记长度
=
/(1-
),
越小对应的
越小,意味着对信号的非平稳性的跟踪越好,而
太小
会小于信号每个平稳段的有效时间,这时不能充分利用所有能够获取的取样数据,结果算出的权矢量w(n)将会受到噪声的影响。
是遗忘因子,其值越接近1表明算法越趋近最小均方准则,那么其曲线的变化就越平稳(上面已经分析了)。
对于平稳模型
的最佳值为1。
但权系收敛后噪声越小,当
为1时看出收敛后是几乎没有噪声的,因为有效记忆长度此时为无穷大。
题目二
题目:
WriteasmallMATLABprogramthatimplementsthepthorderLevinson-Durbin(L-D).Run/TesttheprogramusingaAR
(2)process(b0=1,a1=0,a2=0.81)andanMA
(2)(bn=1,1,1)process-about1000samples.UseL-Dwithp=2(fortheAR)and10(fortheMA).PlottheARspectraproducedinthetwocaseswithL-D.Listthedirectformandthereflectioncoefficientsinatable.ProfiletheL-D(totalnumberofcomputationsforapthorderL-D).
2.1数据模型
N=1000,根据题目采样数据来验证AR
(2)和MR
(2)模型
uu=wgn(1000,1,1);%生成uu(n)是一个均值为零,功率为一的白噪声
P1=2;%AR计算的阶数,与X的自相关函数的阶数少一
P2=10;%10阶L-D算法对MA产生的信号模型进行谱估计
aa0=1;在计算Dk时初始条件ak,0=1这里当k=0时既有aa0,从而初始Dk
Rxx(n)=mean(X’*X)%自相关函数
an=[100.81]%题目所给的参数,以便后续计算x(n)的功率谱密度用到
ssigma=zeros(1,p);%sigma数组的平方
D=zeros(1,p);%Dk数组初始化
arefa=zeros(1,p);%arefa数组初始化
aa=zeros(p,p);%系数数组初始化
2.2L_D算法模型
Levinson-Durbin算法是一种利用迭代的方法计算Yule-walker方程,从而得到AR模型参数的方法。
求解AR模型的参数只需要求解AR模型的Yule-Walker方程,常用解法,如高斯消元法需要的运算量数量级为p3而Levinson-Durbin算法的运算量数量级为p2。
这是一种按阶次进行递推的算法,即首先以AR(0)和AR
(1)模型参数作为初始条件,计算AR
(2)模型参数,然后根据这些参数计算AR(3)模型参数,等等,一直到计算出AR(p)模型参数为止。
fork=1:
p-1
ssigma(k)=Rxx
(1)+sum(a(k,1:
k).*Rxx(2:
k+1));%sigma的平方的计算
D(k)=aa0*R(k+2)+sum(aa(k,1:
k).*fliplr(Rxx(2:
k+1)));%DK的计算
arefa(k+1)=D(k)/(ssigma(k)+eps);%反射系数r的计算
ssigma(k+1)=(1-arefa(k+1)^2)*ssigma(k);%sigma平方的迭代公式
aa(k+1,1:
k)=aa(k,1:
k)-arefa(k+1)*fliplr(a(k,1:
k));
aa(k+1,k+1)=-arefa(k+1);系数的迭代
End
2.3源代码(见附录)
2.4仿真结果及分析
每次计算所得的反射系数如下表所示
r
p
1
2
3
4
5
6
7
8
9
10
AR(10^-4)
-6
-8152
MA(10^-4)
6748
-2502
-2479
3000
-1053
-1349
2306
-631
-1064
1743
2.4.1对AR模型的估计
图3信号x(n)的波形
图4x(n)的功率谱密度
图5AR2估计后的功率谱密度
由图四中是有噪声信号生成的x(n)的功率谱密度,而图5是2阶AR模型估计后的功率谱密度,由两图比较可以看出经过AR模型后x(n)的功率谱几乎没发生变化,从而说明此时AR(k)模型的阶选的较为合适。
阶选的太低,功率谱收到的平滑太厉害。
2.4.2对MA模型的估计
图6信号x(n)的波形
图7信号x(n)的功率谱密度
图8AR(10)估计后的功率谱密度
由图7所示为有噪声信号生成的x(n)的功率谱密度,而图8是10阶AR模型估计后的功率谱密度,由两图比较可以看出采用10阶的L_D算法去估计MR信号的时候,信号的功率谱会有很大的差距。
说明AR模型的阶的选取问题。
比较图5和图8,两者阶数不同,说明阶数选的太低,功率谱收到的平滑太厉害,阶选的太高,固然会提高谱估计的分辨率,但同时会产生虚假的谱峰或谱的细节。
因此,要估计AR(p)的过程,就应该把AR(k)模型的阶选的等于大于p,即K≧p,但k不能太大。
为验证次想法,选取了p=200来进行仿真,结果如图:
这是一幅p=200的AR模型的频率谱,图中出现了许多虚假峰。
附录
题目2代码
%AR
(2)产生数据,用2阶的LD算法计算,输出估计的功率谱
clearall
closeall
%由已知信号模型产生信号
N=1000;
uu=wgn(1000,1,1);%生成uu(n)是一个均值为零,功率为一的白噪声
x
(1)=uu
(1);
x
(2)=uu
(2);
forkk=3:
N
x(kk)=-0.81*x(kk-2)+uu(kk);
end%产生指定的AR
(2)模型的x(n)
%L-D算法
p=2;%%AR计算的阶数
Rxx=zeros(1,p+1);%自相关矩阵的计算
sample_num=length(x);
forn=1:
p+1
Rxx(n)=mean(x(1:
sample_num-n+1).*x(n:
sample_num));
end
%L-D算法数据初始化,MATLAB运行会更快
ssigma=zeros(1,p);%sigma数组的平方
D=zeros(1,p);%Dk数组初始化
arefa=zeros(1,p);%arefa数组初始化
aa=zeros(p,p);%系数数组初始化
%第一组数据初始化
aa0=1;
D0=aa0*Rxx(1+1);
ssigma0=Rxx
(1);
arefa
(1)=D0/(ssigma0+eps);
aa(1,1)=-arefa
(1);
fork=1:
p-1
%注意,这里的Rxx(k+1)就是课本上的Rxx(k)
ssigma(k)=Rxx
(1)+sum(a(k,1:
k).*Rxx(2:
k+1));
D(k)=aa0*R(k+2)+sum(aa(k,1:
k).*fliplr(Rxx(2:
k+1)));
arefa(k+1)=D(k)/(ssigma(k)+eps);
ssigma(k+1)=(1-arefa(k+1)^2)*ssigma(k);
aa(k+1,1:
k)=aa(k,1:
k)-arefa(k+1)*fliplr(a(k,1:
k));
aa(k+1,k+1)=-arefa(k+1);
end
ap=-a(p,:
);
%计算功率谱密度
w=[-pi:
0.1:
pi];
s=zeros(1,length(w));
id=[1:
1:
p];
forn=1:
length(w)
s(n)=ssigma(p)/(norm(1+ap.*exp(-j*w(n)*id))^2);
end
%计算输入信号x(n)的功率谱
an=[100.81];
forn=1:
N
t(n)=an(2:
3)*exp(-j*w(n)*[1:
2])';
SxxW(n)=1/abs(1+t(n))^2;
end
%画图
figure
(1)
plot(x);
title('信号x(n)');
xlabel('n');
ylabel('x(n)');
plot(w,SxxW);
title('x(n)的功率谱密度');
xlabel('PI');
ylabel('Sxx(w)');
figure
(2)
plot(w,s);
title('AR
(2)估计后的功率谱估计');
xlabel('PI');
ylabel('S(w)');
%对MA模型的估计
clearall
closeall
N=1000;
x=zeros(1,N);
u=randn(1,N);
forkk=3:
N
x(kk)=u(kk)+u(kk-1)+u(kk-2);
end
%L-D算法
p=10;
%L-D算法数据初始化,MATLAB运行会更快
ssigma=zeros(1,p);%sigma数组的平方
D=zeros(1,p);%Dk数组初始化
arefa=zeros(1,p);%arefa数组初始化
aa=zeros(p,p);%系数数组初始化;
%第一组数据初始化
aa1=1;
D1=aa1*Rxx(1,1);
ssigma=Rxx(1,1);
arefa
(1)=D1/(ssigma+eps);
aa(1,1)=-arefa
(1);
[X,Rxx]=corrmtx(x,p);
aa2=-Rxx(1,2)/Rxx(1,1);
ssigma
(2)=(Rxx(1,1)^2-Rxx(1,2)^2)/Rxx(1,1);
fork=2:
p
[X,Rxx]=corrmtx(x,k);
D(k)=Rxx(k+1,:
)*[ak(1:
k)0]';
arefa(k+1)=D(k)/ssigma(k);
ssigma(k+1)=(1-arefa(k+1)^2)*ssigma(k);
aa1=[aa(1:
k)0]-arefa(k+1)*[0aa(k:
-1:
1)];
aa=aa1;
end
w=linspace(-pi,pi,1000);
%计算功率谱密度
sum=0;
forn=2:
p+1
sum=sum+ak(n)*exp(-(n-1)*w*1j);
end
Sw=sigma(p+1)./abs((1+sum).^2);
bn=[111];
forn=1:
N
SxxW(n)=abs(bn*exp(-j*w(n)*[0:
2])')^2;
end
figure
(1)%画图
plot(x);
title('信号x(n)');
xlabel('n');
ylabel('x(n)');
figure
(2)
plot(w/pi,SxxW,'b');
title('x(n)的功率谱密度');
xlabel('PI');
ylabel('Sxx(w)');
figure(3)
plot(w/pi,Sw);
title('AR(10)估计后的功率谱密度,对MA
(2)近似');
xlabel('PI');
ylabel('S(w)');
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- adsp 实验 报告 材料 Lms 算法