LMS算法波束形成的基礎(chǔ)仿真分析(總6頁

1、LMS算法應(yīng)用于波束形成的仿真分析1 實驗原理LMS最小均方誤差的方法是由最速下降法推導而出。最速下降法需要求出其梯度的精確值，要求輸入信號和期望信號平穩(wěn)，且 (Rak=抽頭輸入向量u(n)與期望響應(yīng)d(n)的互相關(guān)向量；Rxx=抽頭輸入向量u(n)的相關(guān)矩陣；W=抽頭權(quán)向量)要首先估計 和 ，這給具體實現(xiàn)帶來很大困難，因此該算法還不是真正意義的自適應(yīng)濾波算法，但討論最陡下降法是有意義，由最陡下降法可以很直觀地導出一類自適應(yīng)濾波算法 - LMS算法。LMS算法的基本思想：調(diào)整濾波器自身參數(shù)，使濾波器能夠自適應(yīng)地跟蹤這種輸入信號的變化，實現(xiàn)最優(yōu)濾波。當橫向濾波器運行在實數(shù)據(jù)的情況下，該算法大體上

2、可描述為：抽頭權(quán)向量更新值=老的抽頭權(quán)向量值+學習速率參數(shù)*抽頭輸入向量*誤差信號其中誤差信號定義為期望向量與抽頭輸入向量所產(chǎn)生的橫向濾波器的實際向量之差設(shè)輸入信號為u（n）,LMS算法理論推導過程如下：濾波器輸出y(n)為： （1）由誤差定義得： （2）使用最小均方誤差法，得代價函數(shù)為均方誤差為： （3）式(3)中J 是濾波器的系數(shù)k w (k = 0,1,2,)的函數(shù)。通過選擇最優(yōu)的系數(shù)，使J 達到最小值。定義梯度向量為J ， （4） 另外，最陡下降迭代方程為： （5）LMS 是直接利用單次采樣數(shù)據(jù)獲得的e2(n)來代替均方誤差J(n)，從而進行梯度估計，每次迭代時計梯度估計為： （6）式

3、（6）代入式（5），得到系數(shù)向量自適應(yīng)迭代法： （7）式（7）稱最小均方自適應(yīng)算法LMS。2 實驗內(nèi)容假定陣元間距為的陣元天線陣，接收信號到達角為，干擾信號的到達角為。利用MATLAB編寫LMS程序求解期望的權(quán)值。假定期望接收的信號向量為，其中，。假定干擾信號向量為，其中。兩個信號在時間間隔幾乎是正交的。令期望信號為。設(shè)初始的天線陣權(quán)值全為零。采用100次迭代。使用MATLAB：(1) 令步長；(2) 計算100次迭代的8個天線陣權(quán)值；(3) 畫出權(quán)值幅度與迭代次數(shù)的關(guān)系；(4) 畫出期望信號和天線陣輸出信號；(5) 畫出平方誤差；(6) 利用計算出的最終權(quán)值畫出陣因子方向圖。實驗得出最有天線

4、陣的權(quán)值為：權(quán)的幅度與迭代次數(shù)的關(guān)系如圖1所示，可以看出隨著迭代次數(shù)的增加，權(quán)值的幅度逐漸上升，在迭代次數(shù)60次左右后，權(quán)值的幅度逐漸趨于穩(wěn)定。圖2所示是經(jīng)過約60次迭代后，天線陣輸出信號如何獲得并跟蹤期望信號的，在迭代過程中，開始天線陣輸出信號并沒有跟蹤上期望信號，隨著迭代次數(shù)的增加，天線陣輸出信號逐漸跟蹤上期望信號，獲得很好的輸出效果。圖3表示，在經(jīng)過60次迭代之后，均方誤差收斂到接近于零。開始的時候均方誤差值達到最大，隨著迭代次數(shù)的增加，均方誤差值逐漸下降。最終在迭代大概60次后，均方誤差值趨近于0。圖4表明加權(quán)后的天線陣在30的期望方向上達到了峰值，在-60的干擾方向上為零值。圖1 天

5、線陣權(quán)值的幅度圖2 期望信號的獲得與跟蹤圖3 平方誤差圖4 陣因子方向圖3 項目源碼%- Givens -%N = 8; d = .5; thetaS = 30*pi/180; thetaI = -60*pi/180;%- Desired Signal & Interferer -%T = 1E-3; t = (1:100)*T/100; it = 1:100;S = cos(2*pi*t/T);I = randn(1,100); %- Create Array Factors for each users signal for linear array -%vS = ; vI = ;i =

6、1:N;vS = exp(1j*(i-1)*2*pi*d*sin(thetaS).;vI = exp(1j*(i-1)*2*pi*d*sin(thetaI).;%- Solve for Weights using LMS -%w = zeros(N,1);X = vS + vI;Rxx = X*X;mu = 1/(4*abs(trace(Rxx);wi = zeros(N,max(it);for n = 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 +

7、 mu*conj(e)*x; wi(:,n) = w;endw = w/w(1); % normalize results to first weight%- Display Weights -%disp( )disp(%-%)disp( )disp( The weights for the N = ,num2str(N), ULA are:)disp( )for m = 1:length(w) disp( w,num2str(m), = ,num2str(w(m),3)enddisp( )%- Plot Results -% 1.) Plot Weights vs. Iteration No

8、.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(迭代次數(shù)), ylabel(|權(quán)值|)% 2.) Plot Signal acquisition and trackingfigure(2)plot(it,S,k,it,real(y),k-)xlabel(迭代次數(shù)), ylabel(信號)legend(期望信號,陣列輸出)% 3.) Plot MSEfigure(3), plot(it,esave,k)xlabel(迭代次數(shù)), ylabel(|e|2)% 4.) Plot Array Factortheta = -pi/2:.01:pi/2;AF = 0;for i = 1:N AF = AF + w(i).*exp(j*(i-1)*2*pi*d*sin(theta);endfigure(4),