LMS算法波束形成的基础仿真分析

LMS算法应用于波束形成的仿真分析1实验原理LMS最小均方误差的方法是由最速下降法推导而出。最速下降法需要求出其梯度的精确值，要求输入信号和期望信号平稳，且22jxxdxRWR(Rak=抽头输入向量u(n)与期望响应d(n)的互相关向量；Rxx=抽头输入向量u(n)的相关矩阵；W=抽头权向量)要首先估计xxR和dxR，这给具体实现带来很大困难，因此该算法还不是真正意义的自适应滤波算法，但讨论最陡下降法是有意义，由最陡下降法可以很直观地导出一类自适应滤波算法---LMS算法。LMS算法的基本思想：调整滤波器自身参数，使滤波器能够自适应地跟踪这种输入信号的变化，实现最优滤波。当横向滤波器运行在实数据的情况下，该算法大体上可描述为：抽头权向量更新值=老的抽头权向量值+学习速率参数*抽头输入向量*误差信号其中误差信号定义为期望向量与抽头输入向量所产生的横向滤波器的实际向量之差设输入信号为u（n）,LMS算法理论推导过程如下：滤波器输出y(n)为：10()()Nkkynwunk0,1,2...n（1）由误差定义得：()()-()endnyn（2）使用最小均方误差法，得代价函数为均方误差为：2[()]JEen（3）式(3)中J是滤波器的系数kw(k=0,1,2,…)的函数。通过选择最优的系数，使J达到最小值。定义梯度向量为∇J，∇2[()]()2[()]2[()()]kkkJEenenJEenEunken0,1,2...k（4）另外，最陡下降迭代方程为：()1wnwnJn（5）LMS是直接利用单次采样数据获得的e2(n)来代替均方误差J(n)，从而进行梯度估计，每次迭代时计梯度估计为：22()()[()()()()()2()()()]()()TTTenJndnwnununwndnunwnwnwn2()()()2()()Tununwndnun2[()()()]()2()()Tdnunwnunenun（6）式（6）代入式（5），得到系数向量自适应迭代法：1()2wnwnJnwnenun（7）式（7）称最小均方自适应算法LMS。2实验内容假定阵元间距为/2d的8M阵元天线阵，接收信号到达角为030，干扰信号的到达角为160。利用MATLAB编写LMS程序求解期望的权值。假定期望接收的信号向量为0()()skskxa，其中()cos(2**()/)skpitkT，1msT，(1:100)*/100tT。假定干扰信号向量为1()()ikikxa，其中()randn(1,100)ik。两个信号在时间间隔T几乎是正交的。令期望信号为()()dksk。设初始的天线阵权值全为零。采用100次迭代。使用MATLAB：(1)令步长0.02；(2)计算100次迭代的8个天线阵权值；(3)画出权值幅度与迭代次数的关系；(4)画出期望信号()sk和天线阵输出信号()yk；(5)画出平方误差2e；(6)利用计算出的最终权值画出阵因子方向图。实验得出最有天线阵的权值为：权的幅度与迭代次数的关系如图1所示，可以看出随着迭代次数的增加，权值的幅度逐渐上升，在迭代次数60次左右后，权值的幅度逐渐趋于稳定。图2所示是经过约60次迭代后，天线阵输出信号如何获得并跟踪期望信号的，在迭代过程中，开始天线阵输出信号并没有跟踪上期望信号，随着迭代次数的增加，天线阵输出信号逐渐跟踪上期望信号，获得很好的输出效果。图3表示，在经过60次迭代之后，均方误差收敛到接近于零。开始的时候均方误差值达到最大，随着迭代次数的增加，均方误差值逐渐下降。最终在迭代大概60次后，均方误差值趋近于0。图4表明加权后的天线阵在30°的期望方向上达到了峰值，在-60°的干扰方向上为零值。图1天线阵权值的幅度图2期望信号的获得与跟踪图3平方误差图4阵因子方向图3项目源码%%-----Givens-----%%N=8;d=.5;thetaS=30*pi/180;thetaI=-60*pi/180;%%-----DesiredSignal&Interferer-----%%T=1E-3;t=(1:100)*T/100;it=1:100;S=cos(2*pi*t/T);I=randn(1,100);%%-----CreateArrayFactorsforeachuser'ssignalforlineararray-----%%vS=[];vI=[];i=1:N;vS=exp(1j*(i-1)*2*pi*d*sin(thetaS)).';vI=exp(1j*(i-1)*2*pi*d*sin(thetaI)).';%%-----SolveforWeightsusingLMS-----%%w=zeros(N,1);X=vS+vI;Rxx=X*X';mu=1/(4*abs(trace(Rxx)));wi=zeros(N,max(it));forn=1:length(S)x=S(n)*vS+I(n)*vI;y(n)=w'*x;e=conj(S(n))-y(n);esave(n)=abs(e)^2;w=w+mu*conj(e)*x;wi(:,n)=w;endw=w/w(1);%normalizeresultstofirstweight%%-----DisplayWeights-----%%disp('')disp('%------------------------------------------------------------------------%')disp('')disp(['TheweightsfortheN=',num2str(N),'ULAare:'])disp('')form=1:length(w)disp(['w',num2str(m),'=',num2str(w(m),3)])enddisp('')%%-----PlotResults-----%%%1.)PlotWeightsvs.IterationNo.wi=wi.';figure(1),plot(it,abs(wi(:,1)),'kx',it,abs(wi(:,2))...,'ko',it,abs(wi(:,3)),'ks',it,abs(wi(:,4)),'k+',...it,abs(wi(:,5)),'kd','markersize',2)xlabel('迭代次数'),ylabel('|权值|')%2.)PlotSignalacquisitionandtrackingfigure(2)plot(it,S,'k',it,real(y),'k--')xlabel('迭代次数'),ylabel('信号')legend('期望信号','阵列输出')%3.)PlotMSEfigure(3),plot(it,esave,'k')xlabel('迭代次数'),ylabel('|e|^2')%4.)PlotArrayFactortheta=-pi/2:.01:pi/2;AF=0;fori=1:NAF=AF+w(i)'.*exp(j*(i-1)*2*pi*d*sin(theta));endfigure(4),plot(theta*180/pi,abs(AF)/max(abs(AF)),'k')xlabel('AOA(deg)'),ylabel('|AF_n|')axis([-909001]),gridonset(gca,'xtick',[-90-60-300306090])