自適應(yīng)濾波LMS與RLS的matlab實現(xiàn)

1、MATLAB仿真實現(xiàn)LMS和RLS算法的二階AR模型及仿真結(jié)果分析一、題目概述：二階AR模型如圖1a所示，可以如下差分方程表示： （1）圖1a其中，v(n)是均值為0、方差為0.965的高斯白噪聲序列。a1，a2為描述性參數(shù)，設(shè)x(-1)=x(-2)=0，權(quán)值w10=w20=0，=0.04推導(dǎo)最優(yōu)濾波權(quán)值（理論分析一下）。按此參數(shù)設(shè)置，由計算機仿真模擬權(quán)值收斂曲線并畫出，改變步長在此模擬權(quán)值變化規(guī)律。對仿真結(jié)果進行說明。應(yīng)用RLS算法再次模擬最優(yōu)濾波權(quán)值。解答思路：（1）高斯白噪聲用normrnd函數(shù)產(chǎn)生均值為0、方差為0.965的正態(tài)分布隨機1*N矩陣來實現(xiàn)。隨后的產(chǎn)生的信號用題目中的二階A

2、R模型根據(jù)公式（1）產(chǎn)生，激勵源是之前產(chǎn)生的高斯白噪聲。（2）信號長度N取為2000點，用以觀察濾波器權(quán)值變化從而估計濾波器系數(shù)，得到其收斂值。（3）仿真時分別仿真了單次LMS算法和RLS算法下的收斂性能以及100次取平均后的LMS和RLS算法的收斂性能，以便更好的比較觀察二者的特性。（4）在用不同的分別取3個不同的值仿真LMS算法時，值分別取為0.001,0.003,0.006；用3個不同的值仿真RLS算法時值分別取為1,0.98,0.94，從而分析不同步長因子、不同遺忘因子對相應(yīng)算法收斂效果的影響。二、 算法簡介1自適應(yīng)算法的基本原理自適應(yīng)算法的基本信號關(guān)系如下圖所示：圖 1b 自適應(yīng)濾波

3、器框圖輸入信號x(n)通過參數(shù)可調(diào)的數(shù)字濾波器后產(chǎn)生輸出信號y(n)，將其與參考信號d(n)進行比較，形成誤差信號e(n)。e(n)通過某種自適應(yīng)算法對濾波器參數(shù)進行調(diào)整，最終是e(n)的均方值最小。當(dāng)誤差信號e(n)的均方誤差達到最小的時候，可以證明信號y(n)是信號d(n)的最佳估計。2. LMS算法簡介LMS算法采用平方誤差最小的原則代替最小均方誤差最小的原則，信號基本關(guān)系如下： （2）寫成矩陣型式為： （3）式（3）中，W(n) 為n 時刻自適應(yīng)濾波器的權(quán)值，N為自適應(yīng)濾波器的階數(shù)，本設(shè)計中取為2000；X( n) 為n 時刻自適應(yīng)濾波器的參考輸入矢量，由最近N 個信號采樣值構(gòu)成，；d

4、 ( n) 是期望的輸出值；e ( n) 為自適應(yīng)濾波器的輸出誤差調(diào)節(jié)信號(簡稱失調(diào)信號) ；是控制自適應(yīng)速度與穩(wěn)定性的增益常數(shù)，又叫收斂因子或步長因子。3 RLS算法簡介RLS算法是用二乘方的時間平均的最小化準(zhǔn)則取代最小均方準(zhǔn)則，并按時間進行迭代計算。其基本原理如下所示：遺忘因子，它是小于等于1的正數(shù)。參考信號，也可稱為期望信號。第n次迭代的權(quán)值。均方誤差。RLS算法的準(zhǔn)則為： （4）上式越舊的數(shù)據(jù)對的影響越小。通過計算推導(dǎo)得到系數(shù)的迭代方程為： wn=wn-1+kne*(n) （5）式(5)中，增量k(n)和誤差e*(n)計算公式如下： (6) e*n=d*n-xTn*w(n-1) （7）

5、 式（6）中T(n)= R-1n，也就是當(dāng)前時刻自相關(guān)矩陣的逆。按如下方式更新： = (-* )/ （8）由上邊分析可知，RLS算法遞推的步驟如下：1. 在時刻n，和也已經(jīng)存儲在濾波器的相應(yīng)器件中2. 利用公式（5）、（6）、（7）和（8）計算T(n)、w(n)、k(n)、e*n，并得到濾波器的輸出相應(yīng)和誤差即： （9） （10）3. 進入第次迭代這樣做的優(yōu)點是收斂速度快，而且適用于非平穩(wěn)信號的自適應(yīng)處理缺點是每次迭代時都要知道輸入信號和參考信號，計算量比較大三、 仿真過程仿真過程按照如下過程進行（1） 信號產(chǎn)生：首先產(chǎn)生高斯白噪聲序列w(n)，然后將此通過一個參數(shù)為a1=-0.195，a2=

6、0.95簡單的二階自回歸濾波器生成信號。（2）將step（1）生成的信號通過LMS和RLS自適應(yīng)濾波器進行處理（3）通過改變值對w1,w2收斂速度的影響來分析LMS算法的性能以及通過改變值對w1,w2收斂速度的影響來分析RLS算法的性能。（4）繪制相應(yīng)圖形曲線四、 仿真以及結(jié)果分析 信號和高斯白噪聲波形如圖2所示：圖2 信號和高斯白噪聲波形 圖2中，上邊的圖形為信號波形，下邊的為加性高斯白噪聲。圖3（a）LMS算法下單次收斂曲線圖3(b) LMS算法下百次平均收斂曲線分析1：圖3中，a、b兩幅圖分別為單次實現(xiàn)的LMS算法下最優(yōu)權(quán)值變化過程和100次仿真實現(xiàn)后取平均值做的圖，兩個權(quán)值初始值由已知

7、條件設(shè)置為0，之后收斂到兩個定值。 a圖展現(xiàn)了濾波器權(quán)系數(shù)迭代更新的過程，可以發(fā)現(xiàn)其并不是平滑的變化，而是隨機起伏的，跟最陡下降法不一樣，這充分說明了其權(quán)向量是一個隨機過程向量,梯度噪聲對其產(chǎn)生了一定的影響。b圖給出了100次獨立實驗得到的平均權(quán)向量的估計，即() ，其中是第t次獨立實驗中第n次迭代得到的權(quán)向量，T是獨立實驗次數(shù)?？梢园l(fā)現(xiàn)，多次獨立實驗得到的平均權(quán)向量的估計平滑了隨機梯度引入的梯度噪聲，十分接近理論最優(yōu)權(quán)向量，曲線足夠平滑，噪聲影響很小。圖3中a、b兩圖皆有較好的收斂特性，即使1次實現(xiàn)也能很好的逼近最優(yōu)權(quán)值。圖4（a）RLS算法下單次收斂曲線圖4（b）RLS算法下百次平均收斂曲

8、線 分析2：圖4中，a、b兩幅圖分別為單次實現(xiàn)的RLS算法下最優(yōu)權(quán)值變化過程和100次仿真實現(xiàn)后取平均值做的圖，a圖中可看出權(quán)值已有較好收斂特性，兩個權(quán)值初始值由已知條件設(shè)置為0，之后收斂到兩個定值：-a1=0.195和-a2=-0.95，但是曲線不夠平滑、噪聲較大，b圖經(jīng)過取平均后噪聲影響已經(jīng)很小。單次實現(xiàn)和多次實現(xiàn)的聯(lián)系與區(qū)別也與LMS算法基本相同。但是可以清晰地看出RLS算法收斂極。將兩種算法整合于一個圖中可以看出鮮明的對比。圖5 LMS算法和RLS算法收斂曲線對比分析3：從圖5可以看出，在和一定的情況下，RLS比LMS具有更快的收斂速度。在算法的前期收斂段，RLS算法的收斂速度要明顯高

9、于LMS算法。 總的來說，LMS算法的收斂性和步長有關(guān)，受協(xié)方差矩陣的特征根影響，然而RLS算法一定是收斂的，且收斂速度很快，不過其迭代過程由于要求增益因子和逆矩陣等，使得其算法復(fù)雜度高，計算量比較大。接下來修改步長值，觀察曲線的收斂情況。圖6給出了=0.01, =0.03和=0.06三種情況下的最優(yōu)權(quán)值變化曲線，由此分析不同步長對曲線收斂性產(chǎn)生的影響。圖6 對LMS收斂速度的影響分析4： 可以看出LMS算法中的步長參數(shù)決定抽頭權(quán)向量在每步迭代中的更新量，是影響算法收斂速度的關(guān)鍵參數(shù)，其決定了LMS算法學(xué)習(xí)過程的快慢。 圖6展示了當(dāng)收斂步長值變小時，其均方誤差的收斂速度也相應(yīng)減慢，降低對實變系

10、統(tǒng)的跟蹤速度，不能及時調(diào)整至最優(yōu)權(quán)值。 由前面的理論推導(dǎo)我們可以知道，在這個范圍內(nèi)，越大，均方誤差收斂速度越快，收斂速度和步長因子成正比?？墒侨绻笥谶@個范圍會造成不穩(wěn)定，較大步長會造成較大的穩(wěn)態(tài)誤差，帶來算法的發(fā)散。圖7 對RLS收斂速度的影響 分析5：RLS算法中遺忘因子的作用是對離n時刻越近的誤差加比較大的權(quán)重，而對離n時刻越遠的誤差加比較小的比重。遺忘因子的選擇對RLS算法的性能起決定性的作用。是個小于等于1的數(shù)，如果越小，能量信號就越接近最新的誤差平方，對前面的誤差遺忘的越快，跟蹤效果就越好。但是，遞推RLS算法中的誤差是由期望信號決定的，如果很小則誤差信號對期望信號的依賴性就會很大

11、，所以，輸出信號就很接近期望信號。從圖7可以看出，遺忘因子越小，系統(tǒng)的跟蹤能力越強，收斂的越快（即變平穩(wěn)得越快），但是收斂值比較大；越大，系統(tǒng)跟蹤能力減弱，收斂時均方誤差也越小，學(xué)習(xí)曲線收斂值越小。所以，我們在用普通遞推RLS算法時，一定要對有個準(zhǔn)確的取值，才能保證系統(tǒng)性能的最佳狀態(tài)。 總結(jié)：LMS算法其優(yōu)點是結(jié)構(gòu)簡單，算法復(fù)雜度低，易于實現(xiàn)，穩(wěn)定性高，便于硬件實現(xiàn)，但是這種算法收斂速度慢，對快速變化的信號不適合。RLS算法是基于最小二乘準(zhǔn)則的精確方法，它的收斂速度快，穩(wěn)定性強， 因此常被應(yīng)用于實時系統(tǒng)識別和快速啟動的信道均衡。附件（程序代碼）：%(1)信號序列與高斯白噪聲的產(chǎn)生%參數(shù)初始化a

12、1=-0.195; %生成信號所用AR（2）濾波器的參數(shù)a2=0.95;N=2000; %信號長度 %信號及白噪聲信號序列的初始化x=zeros(1,N)' %信號的初始化sigmasqu=0.0965;%噪聲方差sigma=sqrt(sigmasqu);%噪聲標(biāo)準(zhǔn)差w=normrnd(0,sigma,N)' %高斯白噪聲的初始化，均值為0，方差為0.0965 x(1)=w(1); %信號前兩點的初始賦值x(2)=-a1*x(1)+w(2); %信號序列的產(chǎn)生for i=3:N x(i)=-a1*x(i-1)-a2*x(i-2)+w(i);%信號由AR（2）產(chǎn)生end %繪制信

13、號和高斯白噪聲波形figure(1)subplot(2,1,1),plot(1:N,x,'b-')axis(0,2000,-5,5)title('基于AR（2）模型產(chǎn)生的信號x');xlabel('信號長度n');ylabel('x(n)');subplot(2,1,2),plot(1:N,w,'r-');axis(0,2000,-5,5)title('基于AR（2）模型產(chǎn)生的高斯白噪聲w(n)');xlabel('信號長度n');ylabel('w(n)'); %-

14、LMS100次疊加% N=2000;M=100;%計算的次數(shù)w1=zeros(N,M);w2=zeros(N,M);I1=eye(2);E1=zeros(N,M);wa1=zeros(N,1);wa2=zeros(N,1);en=zeros(N,1); for k=1:M e=zeros(1,N)'%定義誤差向量 %根據(jù)RLS算法進行遞推 x=zeros(1,N)' %信號的初始化 sigmasqu=0.0965;%噪聲方差 sigma=sqrt(sigmasqu);%噪聲標(biāo)準(zhǔn)差 w=normrnd(0,sigma,N)' %高斯白噪聲的初始化，均值為0，方差為0.09

15、65 x(1)=w(1); %信號前兩點的初始賦值 x(2)=-a1*x(1)+w(2); %信號序列的產(chǎn)生 for i=3:N x(i)=-a1*x(i-1)-a2*x(i-2)+w(i);%信號由AR（2）產(chǎn)生 end L=2; %濾波器長度 u=0.04; %LMS算法下自適應(yīng)增益常數(shù)初始化 wL=zeros(L,N);%LMS濾波器的權(quán)值初始化 for i=(L+1):N X=x(i-1:-1:(i-L); y(i)=X'*wL(:,i); %i時刻輸出信號 e(i)=x(i)-y(i); %i時刻誤差信號 wL(:,(i+1)=wL(:,i)+2*u*e(i)*X; %i時刻

16、濾波器的權(quán)值 end;% a1R=wR(1,1:n); % a1在RLS算法下值的變化% a2R=wR(2,1:n); % a2在RLS算法下值的變化 w1(:,k)=wL(1,1:N)' w2(:,k)=wL(2,1:N)'%將每次計算得到的權(quán)矢量值儲存 E1(:,k)=e(:,1);%將每次計算得到的誤差儲存 %求權(quán)矢量和誤差的M次的平均值 wa1(:,1)=wa1(:,1)+w1(:,k); wa2(:,1)=wa2(:,1)+w2(:,k); en(:,1)=en(:,1)+E1(:,k); end %-止LMS100次疊加% %繪制LMS算法下a1、a2單次及100次

17、平均收斂曲線figure(2)plot(1:N,w1(1:N,M)','r-',1:N,w2(1:N,M)','b-',1:N,-a2,'k-',1:N,-a1,'k-');legend('LMS-W1變化','LMS-W2變化','W1收斂值','W2收斂值',0); % 圖例title('LMS單次實現(xiàn)的權(quán)值收斂曲線');xlabel('信號長度n');ylabel('權(quán)值W'); figure(3)

18、plot(1:N,wa1(1:N,1)./M,'r-',1:N,wa2(1:N,1)./M,'b-',1:N,-a2,'k-',1:N,-a1,'k-');%作出100次計算權(quán)矢量的平均變化曲線xlabel('N');ylabel('W(N)');title('LMS算法下100次仿真W1(n)和W2(n)的平均變化曲線') % %RLS濾波 %-RLS100次疊加% N=2000;M=100;%計算的次數(shù)w3=zeros(N,M);w4=zeros(N,M);I2=eye(2);E

19、2=zeros(N,M);wa3=zeros(N,1);wa4=zeros(N,1);en2=zeros(N,1); for k=1:M e=zeros(1,N)'%定義誤差向量 %根據(jù)RLS算法進行遞推 x=zeros(1,N)' %信號的初始化 sigmasqu=0.0965;%噪聲方差 sigma=sqrt(sigmasqu);%噪聲標(biāo)準(zhǔn)差 w=normrnd(0,sigma,N)' %高斯白噪聲的初始化，均值為0，方差為0.0965 x(1)=w(1); %信號前兩點的初始賦值 x(2)=-a1*x(1)+w(2); %信號序列的產(chǎn)生 for i=3:N x(i

20、)=-a1*x(i-1)-a2*x(i-2)+w(i);%信號由AR（2）產(chǎn)生 end L=2; %濾波器長度 lam=0.98; %RLS算法下lambda取值 wR=zeros(L,N);%權(quán)系數(shù)，初值為0 T=eye(L,L)*10;% %RLS算法下T參數(shù)的初始化,T初始值為10 for i=(L+1):N X=x(i-1:-1:(i-L); K=(T*X)/(lam+X'*T*X);%i時刻增益值 e1=x(i)-wR(:,i-1)'*X; wR(:,i)=wR(:,i-1)+K*e1; %i時刻權(quán)值 y(i)=wR(:,i)'*X;%輸出信號 e(i)=x(

21、i)-y(i);%預(yù)測誤差 T=(T-K*X'*T)/lam; %i時刻的維納解 end; w3(:,k)=wR(1,:)' w4(:,k)=wR(2,:)'%將每次計算得到的權(quán)矢量值儲存 E2(:,k)=e(:,1);%將每次計算得到的誤差儲存 %求權(quán)矢量和誤差的M次的平均值 wa3(:,1)=wa3(:,1)+w3(:,k); wa4(:,1)=wa4(:,1)+w4(:,k); en2(:,1)=en2(:,1)+E2(:,k); end %-止RLS100次疊加% figure(4)plot(1:N,w3(1:N,M)','r-',1:N

22、,w4(1:N,M)','b-',1:N,-a2,'k-',1:N,-a1,'k-');legend('RLS-W1變化','RLS-W2變化','W1收斂值','W2收斂值',0); % 圖例title('RLS單次實現(xiàn)權(quán)值收斂曲線');xlabel('信號長度n');ylabel('權(quán)值W'); figure(5)plot(1:N,wa3(1:N,1)./M,'r-',1:N,wa4(1:N,1)./M,&#

23、39;b-',1:N,-a2,'k-',1:N,-a1,'k-');%作出100次計算權(quán)矢量的平均變化曲線xlabel('N');ylabel('W(N)');title('RLS算法下100次仿真W1(n)和W2(n)的平均變化曲線線') %繪制LMS與RLS算法下a1、a2收斂曲線figure(6)plot(1:N,w1(1:N,1)','r-.',1:N,w2(1:N,1)','b:',1:N,w3(1:N,1)','g-',1:

24、N,w4(1:N,1)','m-',1:N,-a2,'k-',1:N,-a1,'k-');legend('LMS-W1變化','LMS-W2變化','RLS-W1變化','RLS-W2變化','W1收斂值','W2收斂值',0); % 圖例title('LMS與RLS算法單次實現(xiàn)的對比');xlabel('信號長度n');ylabel('權(quán)值W'); %(3)LMS算法下不同u值的參數(shù)收斂曲線wL

25、=zeros(L,N,3);eL=zeros(N,3); % LMS算法下誤差初始化yL=zeros(N,3); % LMS算法下濾波器輸出初始化u=0.01,0.03,0.06; %不同的u值 for j=1:3 for i=(L+1):N yL(i,j)=x(i-1:-1:i-2)'*wL(1:L,i-1,j); eL(i,j)=x(i)-yL(i,j); wL(1:L,i,j)=wL(1:L,i-1,j)+2*u(j)*eL(i,j)*x(i-1:-1:i-L); endenda1L1=wL(1,1:N,1);a1L2=wL(1,1:N,2);a1L3=wL(1,1:N,3);a

26、2L1=wL(2,1:N,1);a2L2=wL(2,1:N,2);a2L3=wL(2,1:N,3); figure(7)plot(1:N,a1L1,'b-',1:N,a1L2,'r:',1:N,a1L3,'g-.',1:N,-a1,'k',1:N,a2L1,'b-',1:N,a2L2,'r:',1:N,a2L3,'g-.',1:N,-a2,'k') % 畫圖顯示不同u值下LMS算法性能差別legend('u=0.01','u=0.03','u=0.06','W收斂值',0) % 圖例title('LMS算法下不同的u值對W