華中科技大學(xué)ADSP仿真作業(yè)(LMSLSLPHDMUSIC)

1、ADSP 實(shí)驗(yàn)報告實(shí)驗(yàn)報告院系：電子與通信學(xué)院作者：*學(xué)號：老師：日期：2015/11/20Huazhong University of Science and Technology 目目 錄錄ADSP 實(shí)驗(yàn)報告實(shí)驗(yàn)報告實(shí)驗(yàn)一實(shí)驗(yàn)一.11.1題目.11.2數(shù)據(jù)模型及參數(shù)選取.11.3算法模型.11.3.1自適應(yīng)濾波器.11.3.2LMS算法簡介.21.3.3LSL算法簡介.21.4實(shí)驗(yàn)仿真.31.4.1仿真結(jié)果.31.4.2結(jié)果分析.51.5實(shí)驗(yàn)小結(jié).5實(shí)驗(yàn)二實(shí)驗(yàn)二.62.1題目.62.2數(shù)據(jù)模型及參數(shù)選取.62.3算法模型.72.3.1PHD算法簡介.72.3.2MUSIC算法簡介.82.4

2、實(shí)驗(yàn)仿真.92.4.1仿真結(jié)果.92.4.2結(jié)果分析.102.5實(shí)驗(yàn)小結(jié).11附錄：源代碼附錄：源代碼.12附錄 A：實(shí)驗(yàn)一 .12附錄 B：實(shí)驗(yàn)二.141實(shí)驗(yàn)一實(shí)驗(yàn)一1.1題目Implement LSL algorithm and LMS algorithm based on figure 3.30(P92) and figure 3.31(P93) .Model and parameters see page 91.課本 91 頁模型及參數(shù)為：序列由 2 階自回歸模型產(chǎn)生：( )x n，是均值為 0、方差為 1 的高斯白)()2()1()(21nwnxanxanx( )w n噪聲序列。用

3、LMS 算法和 LSL 算法來估計(jì)模型參數(shù)、121.558,0.81aa 1a，對兩種算法的性能進(jìn)行比較，并分析 LSL 算法中初始預(yù)測誤差剩余對收2a斂性能的影響。1.2數(shù)據(jù)模型及參數(shù)選取1）高斯白噪聲可以用中的函數(shù)產(chǎn)生均值為 0、方差為 1 的標(biāo)matlabrandn準(zhǔn)正態(tài)分布隨機(jī)矩陣來實(shí)現(xiàn)，用模型產(chǎn)生信號，激勵源是產(chǎn)生的高(2)ARrandn斯白噪聲；2）信號點(diǎn)數(shù) N=500，用 500 個信號來估計(jì)濾波器系數(shù)；3）依據(jù)題意，LMS 算法中選取，并取=0.1,1, 10來觀察初始預(yù)0.005測誤差剩余對 LSL 收斂性能的影響。1.3算法模型1.3.1自適應(yīng)濾波器自適應(yīng)濾波器自適應(yīng)濾波器

4、由參數(shù)可調(diào)的數(shù)字濾波器和自適應(yīng)算法兩部分組成，如Error! Reference source not found.所示。輸入信號通過參數(shù)可調(diào)的數(shù)字濾( )x n波器后產(chǎn)生輸出信號，將其與參考信號進(jìn)行比較，形成誤差信號。( )y n( )d n( )e n通過某種自適應(yīng)算法對濾波器參數(shù)進(jìn)行調(diào)整，最終使的均方值最小。( )e n( )e n此時得到的是對的最佳估計(jì)。( )y n( )d n2自適應(yīng)算法參數(shù)可調(diào)數(shù)字濾波器x(n)d(n)y(n)e(n)-+1.3-1 自適應(yīng)濾波器原理1.3.2LMS 算法簡介算法簡介LMS 算法采用平方誤差最小的原則代替最小均方誤差最小的原則，使用LMS 算法來得

5、到權(quán)系數(shù)的公式如下：(1.3-1) nXnenwnw21其中：(1.3-2) nyndnenXnwnyT*式中，為時刻自適應(yīng)濾波器的權(quán)矢量， w nn，為自適應(yīng)濾波器的階數(shù)；為時刻自 TNnwnwnwnw120,.,N X nn適應(yīng)濾波器的參考輸入矢量，由最近個信號采樣值構(gòu)成，N；是期望的輸出值；為自適應(yīng) TNnxnxnxnX1,.1,( )d n( )e n濾波器的參考響應(yīng)與輸出相應(yīng)之差，即誤差信號；是控制自適應(yīng)速度與穩(wěn)定性的增益常數(shù)，又叫收斂因子或步長因子。1.3.3LSL 算法簡介算法簡介LSL 算法步驟如下：1）初始化0( )0,(0)0,(0)1,(0)(0)bffmmmme n這里

6、是前向和后向預(yù)測誤差剩余的初始值。如果這個值沒有給出，可以任意選擇。2）迭代計(jì)算（按時間） ；1,2n 0020000( )( )( )( )( )(1)( )( )1bfbffe nenx ne nennxnn33）迭代計(jì)算（按階）0,1,1mM11(1)( )( )(1)(1)bfmmmmmenennnn 11( )(1)( )( )(1)bffmmmmbmn enenenn11( )( )( )(1)( )fbbmmmmfmn enenenn211( )( )( )(1)ffmmmbmnnnn211( )( )(1)( )bbmmmfmnnnn21(1)(1)(1)(1)bmmmbmen

7、nnn同階的嵌套著按時間進(jìn)行迭代計(jì)算。1m各階前向和后向反射系數(shù)可分別由以下二式計(jì)算：11( )( )( )fmmfmnknn11( )( )(1)bmmbmnknnM 是給定的濾波器的階。采用 LSL 自適應(yīng)算法，利用給出的數(shù)據(jù)可計(jì)算出格型濾波器的 1 階和 2 階前向和后向反射系數(shù)、和，然后由這些反射系數(shù)計(jì)算1( )fkn1( )bkn2( )fkn2( )bkn出估計(jì)、。經(jīng)推導(dǎo)得1 a2 a 1112 ( )( )( )( )bfba nknkn kn （1.3-3） （1.3-4）22 ( )( )ba nkn1.4實(shí)驗(yàn)仿真1.4.1仿真結(jié)果仿真結(jié)果4圖 1.4-1 AR(2)所產(chǎn)生信

8、號圖 1.4-2 自適應(yīng)權(quán)系數(shù)1 a、2 a的收斂軌跡 5圖 1.4-3 初始預(yù)測誤差剩余對 LSL 收斂性能的影響1.4.2結(jié)果分析結(jié)果分析1）由圖 1.4-2 可看出，一是 LSL 算法和 LMS 算法計(jì)算得到的都收斂1 ( )a n于 1.558，都收斂于-0.81，但是 LSL 算法很明顯地比 LMS 算法收斂得更2 ( )a n快，二是 LMS 算法得到的、的收斂軌跡有所波動，而 LSL 算法所得到收1 a2 a斂軌跡更加平穩(wěn)；2）由圖 1.4-3，在=0.1 和 1 兩種情況下，很快就收斂到，而當(dāng)1 ( )a n1a選為 10 時，收斂稍許慢了一些，但三種情況的收斂速度都要比 LM

9、S 算法快的多。1.5實(shí)驗(yàn)小結(jié)通過本次實(shí)驗(yàn)，對 LMS 算法和 LSL 算法的原理有了更深刻的理解，對影響 LMS 和 LSL 算法性能的一些參數(shù)也有了一定的認(rèn)識，懂得了如何通過調(diào)節(jié)參數(shù)來獲得想要的性能。6實(shí)驗(yàn)二題目給定，是方差為的復(fù)白噪聲。信(2 *0.52)2 *0.54( )( )jnjnx neev n0,1,24n ( )v n2v噪比定義為211010log()vSNRdB1）用 PHD 方法估計(jì)兩個頻率，使用 M+1=3 階自相關(guān)矩陣；2）MUSIC 方法估計(jì)頻率，選取 N=12，M=2。2.2數(shù)據(jù)模型及參數(shù)選取1）依題意，可通過提供的函數(shù)可對復(fù)正弦信號( )x nmatlaba

10、wgn加高斯白噪聲得到。(2 *0.52)2 *0.54( )jnjns nee2）PHD算法取樣自相關(guān)點(diǎn)數(shù)_25Nsample 正弦信號個數(shù)2M 自相關(guān)矩陣階數(shù)13NM 復(fù)白噪聲功率 = 0.012v信噪比211010log20()vSNRdB待估計(jì)頻率 0.5, 0.52f 3）MUSIC算法取樣自相關(guān)點(diǎn)數(shù)_25Nsample 正弦信號個數(shù)2M 自相關(guān)矩陣階數(shù)112NM 復(fù)白噪聲功率 = 0.042v信噪比SNR =13.9794dB待估計(jì)頻率 0.5, 0.52f 72.3算法模型2.3.1PHD 算法簡介算法簡介PHD（Pisarenko Harmonic Decomposition）

11、即諧波分解方法。在的特殊情況下，1NM是矩陣，是 M 階對角陣，信號的自相關(guān)矩陣1 2MEeee(1)MM12mPdiag PPP可表示為SR HsREPE （2.3-1） 設(shè)是中的任一矢量，表示由映射x1MR1MR （2.3-2)sR xy得到的一組矢量，令y （2.3-3）HPE x于是 （2.3-4）1Mi iIye由此看出，的秩最多是 M，若可任意選擇，則秩恰為 M。由于是矩陣，這SRi1M 就是要求他在中。利用式（2.3-3） ，只需證明是滿秩的或=是線性獨(dú)立的。MRHPE()HHPEEP但 （2.3-5）1Mii iiEPPe它不能是零，除非所有是零（假設(shè)所有是非零且所有的頻率不同

12、） 。這就證明了的iPsR秩是且有唯一的零特征值。由于M （2.3-6）2xsvRRI與有相同的特征矢量，而的特征值是，所以的特征值是sRsR12 ,0M xR?，F(xiàn)考慮對于最小特征值的特征矢量，有222212,vvMvvxR2v1Mv （2.3-7）211xMvMR vv由式（2.3-6） ，得到8 （2.3-8）1100sMHMR vEPE v左乘，得到1HMv （2.3-9）1111()()0HHHHHMMMMvEPE vE vP E v由于 P 是正定的，于是得出 （2.3-10）10HME v, （2.3-11）10HiMe v1,2,iM上式可寫成展開形式, （2.3-12）10(1

13、)exp()0MMinvnjwn1,2,iM這就是說，如果給定，就能夠求出多項(xiàng)式1Mv, （2.3-13）10(1),exp()MnMinvnzzjw1,2,iM位于單位圓上的零點(diǎn)，這些零點(diǎn)的角度就是正弦波的頻率。2.3.2MUSIC 算法簡介算法簡介MUSIC（Multiple Signal Classification）即多信號分類算法，這種算法是以下列譜估計(jì) （2.3-14）211( )|MUSICNHii MPfe v的峰作為正弦波頻率的估計(jì)的。式中，表示頻率， 表示信號矢量。理論上，當(dāng)時，feiff，。由于存在估計(jì)誤差，故在正弦波頻率上或者附近將有一個峰，與式iee( )MUSICP

14、f （2.3-14）的 M 個最大峰對應(yīng)的頻率便是正弦波的頻率估計(jì)。式（2.3-14）的最大值也可以利用信號子空間的特征矢量來進(jìn)行計(jì)算。為此，注意到該式分母最小時最大，但分母可表示成 ( )MUSICPf （2.3-15）211 |NNHHHiiii Mi Me vevv e9該式利用了。式（2.3-14）的最大化等效于下式1 NHiiiIvv （2.3-16）21|NHii Me v的最大化?？梢?，無論是原先提出的噪聲子空間特征矢量，還是式（2.3-15）的信號子空間矢量，都可以用來確定正弦波的頻率估計(jì)。2.4實(shí)驗(yàn)仿真2.4.1仿真結(jié)果仿真結(jié)果PHD 算法中解方程結(jié)果（運(yùn)行五次程序）：1f0

15、.45170.45910.45030.45200.45492f0.52550.53690.52470.52570.5298表 2.4-1 解方程求得的頻率估計(jì)值圖 2.4-1 作圖法估計(jì)頻率值101f0.43660.45330.44590.45050.46202f0.51500.52650.52600.52650.5403表 2.4-2 PHD 作圖法捕捉到的頻率值圖 2.4-2 MUSIC 求估計(jì)頻率值1f0.46940.47080.47350.47250.47122f0.51180.51230.51040.51500.5134表 2.4-3 MUSIC 方法捕捉到的頻率值2.4.2結(jié)果分析

16、結(jié)果分析1）PHD 算法中，解方程的方法，求解方程（2.4-12）位于單位圓上的零點(diǎn)，并計(jì)算Z的角度，然后把變換到，再將角頻率除以，就可以得到估計(jì)的頻率值。Zww0,2 )w2f由表 2.1 可以看出，估計(jì)出來的兩個頻率值與實(shí)際值 0.5 和 0.52 存在一定的誤差。PHD 的另一種估計(jì)頻率的方法是，計(jì)算式（2.4-10）的值，在兩個頻率附近，值應(yīng)該為0，但實(shí)際上由于估計(jì)誤差，不會精確地等于 0，如果把計(jì)算的結(jié)果倒數(shù)，那么在估計(jì)的頻率附近，將會出現(xiàn)尖峰，捕捉尖峰所對應(yīng)的橫軸值，就可以得到估計(jì)的頻率值。由圖 2.4-1、表2.4-4 與表 2.4-2 可以看出，兩個峰值頻點(diǎn)與實(shí)際都存在一定的誤

17、差，但都在可以接受的范圍11之內(nèi)。2）圖 2.4-2 為 MUSIC 方法估計(jì)頻率，其方法與 PHD 的第二種方法類似，由表 2.4-5 可看出，該方法估計(jì)的頻率值與標(biāo)準(zhǔn)值 0.5 和 0.52 之間也存在誤差。2.5實(shí)驗(yàn)小結(jié)通過本次實(shí)驗(yàn)，對頻率估計(jì)方法有了一定的了解，掌握了 PHD 和 MUSIC 算法估計(jì)頻率方法。12附錄：源代碼附錄：源代碼附錄 A：實(shí)驗(yàn)一clear;close all;% 由AR(2)模型生成輸入信號% 初始化參數(shù)a1 = 1.558;a2 = -0.81; N = 500; %信號點(diǎn)數(shù) %信號及白噪聲信號序列的初始化x = zeros(1,N); %信號的初始化Noi

18、se = randn(1,N); %白噪聲的初始化，均值為0，方差為1x(1) = Noise(1);x(2) = a1*x(1)+Noise(2);for n = 3:N x(n) = a1*x(n-1)+a2*x(n-2)+Noise(n);end% LMS算法估計(jì)a1、a2的收斂曲線 u = 0.005; %LMS算法下自適應(yīng)增益常數(shù)初始化L = 2;LMS_w = zeros(L,n);for i = (L+1):n X = x(i-1:-1:(i-2); y(i) = X*LMS_w(:,i); %i時刻輸出信號 e(i) = x(i)-y(i); LMS_w(:,(i+1) = L

19、MS_w(:,i)+2*u*e(i)*X; %i時刻濾波器的權(quán)值end LMS_a1 = LMS_w(1,1:n); %LMS算法中權(quán)系數(shù)a1的提取 LMS_a2 = LMS_w(2,1:n);% LSL算法估計(jì)a1、a2的收斂曲線 %初始化相關(guān)的參量M = 3;Delta(1:M,1) = 0; %前后向預(yù)測誤差相關(guān)系數(shù)Sf(1:M,1) = 1; %初始前向預(yù)測誤差剩余Sb(1:M,1) = 1; %初始后向預(yù)測誤差剩余gama(1:M,1) = 1; %角參量%迭代計(jì)算for m = 1:M-113 for n = 2:N eb(1,n) = x(n); ef(1,n) = x(n);

20、Sb(1,n) = Sf(1,n-1)+x(n)*x(n); Sf(1,n) = Sf(1,n-1)+x(n)*x(n); gama(1,n) = 1; Delta(m+1,n) = Delta(m+1,n-1)+eb(m,n-1)*ef(m,n)/gama(m,n-1); ef(m+1,n) = ef(m,n)-Delta(m+1,n)*eb(m,n-1)/Sb(m,n-1); eb(m+1,n) = eb(m,n-1)-Delta(m+1,n)*ef(m,n)/Sf(m,n); Sf(m+1,n) = Sf(m,n)-Delta(m+1,n)*Delta(m+1,n)/Sb(m,n-1);

21、 Sb(m+1,n) = Sb(m,n-1)-Delta(m+1,n)*Delta(m+1,n)/Sb(m,n); gama(m+1,n-1) = gama(m,n-1)-eb(m,n-1)*eb(m,n-1)/Sb(m,n-1); kf(m+1,n) = Delta(m+1,n)/Sf(m,n); kb(m+1,n) = Delta(m+1,n)/Sb(m,n-1); endendkb(:,1) = 0;%初始kb，kf的值都是0，根據(jù)遞推關(guān)系可以算出kf(:,1) = 0;for n=1:N LSL_a1(n)= kb(2,n)-kf(2,n)*kb(3,n); LSL_a2(n) = k

22、b(3,n);end% 繪圖figure;plot(x);xlabel(n);ylabel(x(n);title(x(n) generated by AR model);figure();hold onplot(1:n,LMS_a1,g,1:n,LMS_a2,g);plot(1:n,LSL_a1,r,1:n,LSL_a2,r);plot(a1*ones(1,n),-);plot(a2*ones(1,n),-); text(0,1.65,a1=1.558);text(0,-0.85,a2=-0.81);text(100,1.2,LMS);text(400,1.7,LSL,delta=1);xla

23、bel(n);title(LSL與LMS算法估計(jì)信號模型參數(shù)的性能比較);axis(0 500 -1.5 2);hold off% 初始預(yù)測誤差剩余對LSL收斂性能的影響D = 0.1,1,10; %初始預(yù)測誤差剩余a1 = zeros(3,N); %不同初始值對應(yīng)的a1的矩陣14%迭代計(jì)算for i = 1:3 Sf(1:M,1) = D(i); Sb(1:M,1) = D(i);for m = 1:M-1 for n = 2:N eb(1,n) = x(n); ef(1,n) = x(n); Sb(1,n) = Sf(1,n-1)+x(n)*x(n); Sf(1,n) = Sf(1,n-1

24、)+x(n)*x(n); gama(1,n) = 1; Delta(m+1,n) = Delta(m+1,n-1)+eb(m,n-1)*ef(m,n)/gama(m,n-1); ef(m+1,n) = ef(m,n)-Delta(m+1,n)*eb(m,n-1)/Sb(m,n-1); eb(m+1,n) = eb(m,n-1)-Delta(m+1,n)*ef(m,n)/Sf(m,n); Sf(m+1,n) = Sf(m,n)-Delta(m+1,n)*Delta(m+1,n)/Sb(m,n-1); Sb(m+1,n) = Sb(m,n-1)-Delta(m+1,n)*Delta(m+1,n)/

25、Sb(m,n); gama(m+1,n-1) = gama(m,n-1)-eb(m,n-1)*eb(m,n-1)/Sb(m,n-1); kf(m+1,n) = Delta(m+1,n)/Sf(m,n); kb(m+1,n) = Delta(m+1,n)/Sb(m,n-1); endendkb(:,1) = 0;%初始kb，kf的值都是0，根據(jù)遞推關(guān)系可以算出kf(:,1) = 0;for n=1:N LSL_a1(n) = kb(2,n)-kf(2,n)*kb(3,n); a2(n) = kb(3,n);end a1(i,:) = LSL_a1;endfigure;hold onplot(a1

26、(1,:),g);plot(a1(2,:),b);plot(a1(3,:),r);axis(0 500 -0.5 2);xlabel(n);title(初始預(yù)測誤差剩余delta對LSL收斂性能的影響);legend(delta=0.1,delta=1,delta=10);%附錄 B：實(shí)驗(yàn)二1、PHD 算法clc; clear; % 初始值設(shè)定15N = 3; %自相關(guān)矩陣 Rx 的階數(shù) N_sample = 25; % 取樣點(diǎn)數(shù) M = 2; % 正弦信號的個數(shù) sigma = 0.1; %白噪聲的方差 SNR = 10 * log10(1/(sigma . 2); %求信噪比 % 產(chǎn)生輸入

27、信號n = (0:N_sample-1); %數(shù)據(jù)取樣點(diǎn);s = exp(1i*2*pi*0.5*n)+exp(1i*(2*pi*0.52*n+0.25*pi); %正弦信號x = awgn(s,SNR); %向信號中加入白噪聲% 實(shí)現(xiàn) PHD 算法Rxx = xcorr(x); %求自相關(guān)函數(shù) Rx = toeplitz(Rxx(N_sample:N_sample+N-1); %求自相關(guān)矩陣 V,D = eig(Rx); %求特征值和特征向量 V_min = V(:,1); %最小特征值對應(yīng)的特征向量 %方法 1：解正交方程（2.3 12） ，直接求出頻率值 Z = roots(V_min.); %方程（2.3 12）位于單位圓上的零點(diǎn)w = angle(Z); %零點(diǎn)的角度 for i = 1 : M if w(i) 0 w(i) = w(i) + 2 * pi % angle 函數(shù)求出的角度范圍是-pi, pi), %當(dāng) w 時，需要加 2*pi,將范圍調(diào)整到0, 2*pi)之間 endendf_PHD = w / (2*pi) %求得的頻率值%方法 2：信號矢量的頻率估計(jì),通過方程（2.3 10）用作圖法求估計(jì)頻率值K = 1000; % 在 0-1 之間分 1000 步搜索 f = linspace(0, 1, K); e = zeros(N, K); for