ADSP仿真報告-華中科技大學

1、 ADSP仿真作業(yè) 專業(yè)：電磁場與微波技術(shù)學生成員：劉 航(M201571827) 完成時間：2022年2月21日 目錄目錄1圖目錄21題目131.1題目31.2算法模型31.2.1自適應(yīng)濾波原理31.2.2LMS算法簡介31.2.3LSL算法簡介41.3本題模型51.4仿真過程及結(jié)果分析61.4.1仿真過程61.4.2結(jié)果分析82題目292.1題目92.2算法模型92.2.1特征分解頻率估計原理92.2.2Music算法簡介102.2.3PHD算法簡介102.3仿真過程及結(jié)果分析112.3.1仿真過程112.3.2結(jié)果分析133附錄143.1題目1源程序143.2題目2源程序15圖目錄圖 1

2、.1自適應(yīng)濾波器框圖3圖 1.2 題目一模型5圖 1.3 LSL和LMS算法估計信號模型參數(shù)的性能比較8圖 1.4 初始預(yù)測誤差對LSL收斂性能的影響9圖 2.1 MUSIC算法得到的譜峰131 題目11.1 題目 Implement LSL algorithm and LMS algorithm based on figure 3.30(P92) and figure 3.31(P93). Model and parameters see page 91.P91模型如下：二階自回歸隨機過程x(n)=a1x(n-1)+a2x(n-2)+w(n)，對應(yīng)的濾波器參數(shù)a1=1.558，a2 =-0.

3、81。分別用LMS算法和LSL算法預(yù)測信號x(n)，從而得到對信號模型的兩個參數(shù)值的估計值。1.2 算法模型1.2.1 自適應(yīng)濾波原理自適應(yīng)濾波器由參數(shù)可調(diào)的數(shù)字濾波器和自適應(yīng)算法兩部分組成，其原理如下圖所示： 圖 1.1自適應(yīng)濾波器框圖輸入信號x(n)通過參數(shù)可調(diào)的數(shù)字濾波器后產(chǎn)生輸出信號y(n)，將其與期望輸出信號d(n)進行比較，形成誤差信號e(n)。e(n)通過某種自適應(yīng)算法對濾波器參數(shù)進行調(diào)整，最終使得e(n)的均方值最小。1.2.2 LMS算法簡介在自適應(yīng)算法中，最理想的情況是沿著性能曲面最陡的方向向下搜索曲面的最低點，每次迭代都需要知道性能曲面上某點的梯度值，而實際應(yīng)用中梯度值只

4、能根據(jù)觀測數(shù)據(jù)進行估計。LMS（Least mean square）最小均方算法，是一種更為簡便的方法，它的核心思想是用平方誤差代替均方誤差，從而使得計算量大大減小，易于實現(xiàn)。LMS算法中，梯度用下式近似：可以得到：實際上，只是單個平方誤差序列的梯度，而則是多個平方誤差序列統(tǒng)計平均的梯度，所以LMS算法就是用前者作為后者的近似。得到LMS算法的基本關(guān)系式此時，全系數(shù)的調(diào)整過程是有“噪聲”的，其估計路徑不可能準確的沿著理想的最陡下降的路徑。LMS算法中，信號基本關(guān)系總結(jié)如下：1.2.3 LSL算法簡介基于最小均方誤差(MMSE)準則的算法，如最陡下降法、LMS算法，其缺點是：收斂速度慢；對非平穩(wěn)

5、信號的適應(yīng)性較差。為克服以上缺點，引入“最小二乘(LS)”準則。最小二乘準則，是以誤差的平方和最小作為最佳估計的一種誤差準則。最小二乘濾波(LS濾波)可描述為：調(diào)整濾波器的權(quán)矢量，使得在每個時刻對所有已輸入的信號而言，濾波器輸出的誤差平方和最小。它與LMS算法的區(qū)別是：LMS濾波是以集合平均為基礎(chǔ)的，屬于統(tǒng)計分析方法；LS濾波是一種瞬時分析方法, 即在每個時刻對所有已輸入信號，都要重新評估其平方和，并通過調(diào)整權(quán)矢量使其達到最小。LSL是最小二乘格型算法（Least square lattice），用矢量空間法描述前向和后向線性預(yù)測誤差，用預(yù)測誤差濾波器的格型結(jié)構(gòu)來遞推更新誤差參數(shù)。算法的具體推

6、導比較復(fù)雜，詳見教材的3.14節(jié)及3.15節(jié)。以下歸納LSL自適應(yīng)算法的計算過程與流程：(1) 初始化：；為一個小的正數(shù)，它是給定的前向和后向誤差能量的初值, 若此值未預(yù)先給出,則可任意設(shè)定。(2) 迭代計算（按時間）；(3) 迭代計算（按階數(shù)）；（這里注意同階的DM+1嵌套著按時間的迭代運算）；。(4) 計算各階計算前向和后向反射系數(shù)：；。1.3 本題模型對于本題，其模型框圖如下：圖 1.2 題目一模型根據(jù)照題中所給的2階自回歸模型，由單位方差高斯白噪聲w(n)激勵一個線性時不變?nèi)珮O點系統(tǒng)產(chǎn)生信號x(n)。預(yù)測誤差為：y(n)=x(n)-d(n)，然后再按照某種準則控制預(yù)測誤差，從而自適應(yīng)的

7、調(diào)節(jié)FIR濾波器的權(quán)系數(shù)，使之最終達到最優(yōu)。本題用LMS算法以及LSL算法對預(yù)測誤差進行處理：LMS算法直接計出w1(n)及w2(n)；LSL算法計算出格型濾波器的1階和2階前向和后向反射系數(shù)，然后由這些系數(shù)計算出參數(shù)估計值，具體公式見教材3.244-3.277式。1.4 仿真過程及結(jié)果分析1.4.1 仿真過程(1) 初始化參數(shù)，由AR(2)模型生成x(n)。如下截圖：其中信號點數(shù)取500，AR(2)模型由filter函數(shù)產(chǎn)生，filter函數(shù)具體用法參見matlab幫助。(2) LMS算法實現(xiàn)，程序如下截圖：程序中for循環(huán)的公式見1.2.2節(jié)最后總結(jié)的信號基本關(guān)系的公式。矩陣w的第一行為參

8、數(shù)a1的估計，第二行為參數(shù)a2的估計。(3)LSL算法實現(xiàn)，程序截圖如下：程序中要注意下標值，matlab數(shù)組下標從1開始，對于2階情況，按公式迭代是從m=0,1，故程序中是從m=1,2迭代。下圖中for循環(huán)實現(xiàn)的公式見1.2.3節(jié)總結(jié)的算法流程。計算得到的a1(n)和a2(n)為參數(shù)估計值。(3) 畫圖程序，比較兩個算法的性能。程序截圖如下：將兩個算法得到的參數(shù)估計值的收斂軌跡畫在同一張圖中，用不同顏色表示出來。(4) 對LSL算法，前向和后向預(yù)測誤差初始值選取不同值，畫出收斂圖。這部分直接在LSL算法程序中修改初始參數(shù)，選取=0.1,1.0,10三種情況，得到a1的估計值，然后在同一張圖中

9、畫出來。1.4.2 結(jié)果分析(1) LMS和LSL算法性能比較，程序圖形如下：圖 1.3 LSL和LMS算法估計信號模型參數(shù)的性能比較如圖，4條顏色不同的線表示參數(shù)的收斂軌跡，黃色線為a1和a2的標準值。從圖中可以看出，LSL算法和LMS算法計算得到的a1(n)都收斂于1.558，a2(n)都收斂于-0.81，但是LSL算法明顯地比LMS算法收斂的更快。(2)如下圖，對LSL算法，畫出了對參數(shù)a1的估計a1(n)的收斂軌跡，分別取0.1,1.0,10三種情況?？梢钥闯?，在=0.1和1兩種情況下（圖中藍色線和黑色線），a1(n)很快收斂于1.558，但是在軌跡的開始階段，有一個下沖；=10時（紅

10、色線），收斂速度稍微慢一些，但開始階段沒有明顯地下沖。圖 1.4 初始預(yù)測誤差對LSL收斂性能的影響2 題目22.1 題目Implement PHD algorithm and MUSIC algorithm based on P164 4.25(4)&(5).4.25題目：計算機模擬實驗，給定：x(n) = exp(j2 * 0.5n) + exp(j(2 * 0.52n + /4) + v(n), n = 0,1,2,24v(n)是方差為v2的復(fù)白噪聲。信噪比定義為SNR = 10log10(1/v2)(dB)。(4)用PHD方法估計兩個頻率。使用M+1=3階自相關(guān)矩陣，其元素是有

11、偏自相關(guān)函數(shù)值。(5)用MUSIC方法估計頻率，選取N=12，M=2。計算式(4.155)時要采用更細的頻率間隔。2.2 算法模型2.2.1 特征分解頻率估計原理特征分解的主要思想是：把自相關(guān)矩陣的信息空間分成兩個子空間，即信號子空間和噪聲子空間；這兩個子空間中的矢量的函數(shù)在正弦波頻率上有尖銳的風，據(jù)此即可以估計正弦波的頻率。總結(jié)如下：已知求正弦波頻率估計。，其中Rs矩陣的秩為M，它有N-M個零特征值。所以Rx可以分解為兩個子空間：噪聲子空間（N-M維），它有N-M正交特征矢量，對應(yīng)的特征值均為（即噪聲功率）；信號子空間，有M個正交特征特征矢量，對應(yīng)著M個大特征值。這里有一個重要性質(zhì)：主特征矢

12、量張成的子空間與由信號矢量張成的子空間是相同的。所以，信號矢量與噪聲子空間的所有矢量是正交的，即書上4.148式：這個性質(zhì)是噪聲子空間頻率估計的基礎(chǔ)。2.2.2 Music算法簡介MUSIC算法就是基于上面講的4.148式，定義下面這個譜函數(shù)（即書上4.155式）：式子中，f表示頻率，e表示信號矢量，vi表示噪聲子空間的特征矢量。即在處出現(xiàn)峰值，這些峰值所對應(yīng)的M個f，就是所要估計的正弦信號頻率。2.2.3 PHD算法簡介在N=M+1的特殊情況下，是(M+1)xM矩陣，對Rx的最小特征值對應(yīng)的特征矢量，有如下關(guān)系成立：所以，只要我們得到，就可以求出上式的解，這些解就是待估計的頻率。2.3 仿真

13、過程及結(jié)果分析2.3.1 仿真過程(1) 初始化參數(shù)，產(chǎn)生復(fù)白噪聲和待估計信號如圖，根據(jù)題目要求，取樣點數(shù)N_sample取25，信號個數(shù)M=2，由題目定義的信噪比產(chǎn)生復(fù)白噪聲信號。(2) PHD算法實現(xiàn)A、 根據(jù)取樣的觀測信號x(n)，求出其N=M+1=3階自相關(guān)矩陣。如下截圖：這里使用xcorr函數(shù)和toeplitz函數(shù)。實際用的自相關(guān)矩陣的元素，用的是自相關(guān)函數(shù)的無偏估計（參見xcorr函數(shù)用法，第二個參數(shù)選unbiased）。用toeplitz函數(shù)，產(chǎn)生3階自相關(guān)矩陣，保證得到的自相關(guān)矩陣具有理論自相關(guān)矩陣的性質(zhì)（詳見toeplitz函數(shù)的具體用法）。B、 求出Rx的特征值和特征向量，

14、并找出最小特征值對應(yīng)的特征向量C、 根據(jù)書上式4.152或者式4.154，求解方程，得到帶估計頻率值這里要注意，angle函數(shù)求出的角度范圍是0, 2*pi),當求出的角頻率w <0時，需要加上一個2*pi以調(diào)整到0, 2*pi)之間。其中roots函數(shù)具體用法參見matlab幫助。(3) MUSIC算法實現(xiàn)A、 第一步同PHD算法，求解自相關(guān)矩陣，并得到其特征矢量。其中，自相關(guān)矩陣的階數(shù)N根據(jù)題目要求取N=12。B、 第二步，根據(jù)書上式4.155，用作圖法，求出待估計頻率代碼中用到ginput函數(shù)，來認為捕捉圖像的兩個尖峰，然后輸出了估計得到的兩個頻率。2.3.2 結(jié)果分析(1) 用P

15、HD算法，運行得到頻率的估計值：運行4次，得到的估計頻率如下：，可以看出，PHD算法能基本分辨出0.5和0.52兩個頻率，但誤差較大。更改待估計信號的頻率值，將兩個信號的頻率差距拉大，可以發(fā)現(xiàn)，得到的頻率更準。如下（改成0.5和0.55兩個頻率）：,(2) 用MUSIC算法，捕捉圖像的兩個尖峰，得到頻率的估計值：圖 2.1 MUSIC算法得到的譜峰多運行幾次，得到的頻率值如下：，可以看出，MUSIC算法基本能夠分辨出0.5和0.52兩個頻率。相比于PHD算法，其誤差小一些，估計的更準一些。3 附錄3.1 題目1源程序% 2階自回歸模型參數(shù)估計LSM&LSL算法by劉航 20151108

16、clear;clf;clc;% 初始化參數(shù)%N是信號點數(shù)；N=500;na=1:N;%二階自回歸信號模型a1=1.558;a2=-0.81;x=zeros(N,1);%Wa=randn(N,1);Wa=randn(1,N)'%白噪聲x=filter(1,1 -a1 -a2,Wa);% LMS算法過程L=2; %濾波器長度w=zeros(L,N); %LMS濾波器的系數(shù)mu=0.002; %mu是迭代步長；d=zeros(1,N);for i=L+1:N X=x(i-1) x(i-2); y(i)=X*w(:,i); e(i)=x(i)-y(i); w(:,i+1)=w(:,i)+2*m

17、u*e(i)*X'end% LSL算法過程%初始化相關(guān)的參量M=3;Deta(1:M,1)=0;gama(1:M,1)=1.0;sef(1:M,1)=1;seb(1:M,1)=1;%前向和后向預(yù)測誤差初始值for m=1:M-1 %按階數(shù)迭代（階數(shù)為2,matlab下標從1開始，故M設(shè)為3） for n = 2 : N %包含按時間迭代 eb(1,n)=x(n); ef(1,n)=x(n); seb(1,n)=sef(1,n-1)+x(n)*x(n); sef(1,n)=sef(1,n-1)+x(n)*x(n); gama(1,n)=1; Deta(m+1,n)=Deta(m+1,n-

18、1)+eb(m,n-1)*ef(m,n)/gama(m,n-1); ef(m+1,n)=ef(m,n)-Deta(m+1,n)*eb(m,n-1)/seb(m,n-1); eb(m+1,n)=eb(m,n-1)-Deta(m+1,n)*ef(m,n)/sef(m,n); sef(m+1,n)=sef(m,n)-Deta(m+1,n)*Deta(m+1,n)/seb(m,n-1); seb(m+1,n)=seb(m,n-1)-Deta(m+1,n)*Deta(m+1,n)/seb(m,n); gama(m+1,n-1)=gama(m,n-1)-eb(m,n-1)*eb(m,n-1)/seb(m,

19、n-1); kf(m+1,n)=Deta(m+1,n)/sef(m,n); kb(m+1,n)=Deta(m+1,n)/seb(m,n-1); endendkb(:,1)=0;kf(:,1)=0;for n=1:N a1n(n)=kb(2,n)-kf(2,n)*kb(3,n); a2n(n)=kb(3,n);end% 畫圖程序figure(1),plot(na,w(1,na),'k-','linewidth',2),hold on; plot(na,w(2,na),'b-','linewidth',2),hold on;plot(

20、na,a1n,'g-','linewidth',2);hold on;plot(na,a2n,'r-','linewidth',2);legend('表示LMS中a1','表示LMS中a2','表示LSL中a1','表示LSL中a2',0);hold on;plot(1.558*ones(n),'y-');plot(-0.81*ones(n),'y-');text(0,1.558,'1.558','color&#

21、39;,'r'),text(0,-0.81,'-0.81','color','r');title('LSL 與 LMS 性能參數(shù)比較') ;grid on;xlabel('n');ylabel('a1(n),a2(n)');% LSL不同初始值收斂性能% %更改sef(1:M,1)=1;seb(1:M,1)=1的值為0.1和10(得到a1_1n,a1_2n)figure(2),plot(na,a1n,'k-','linewidth',2),hold o

22、n;plot(na,a1_1n,'r-','linewidth',2),hold on;plot(na,a1_2n,'b-','linewidth',2),hold on;legend('1.0','10','0.1',0);grid on;%3.2 題目2源程序% PHD算法by劉航20151108% 初始化參數(shù)，產(chǎn)生復(fù)白噪聲噪聲和待估計信號clc;clear all;N_sample = 25;% 取樣自相關(guān)點數(shù)M = 2; % 正弦信號的個數(shù)sigma = 0.1; v_rea

23、l = sigma * randn(N_sample, 1) / sqrt(2); % 白噪聲的實部v_imag = sigma * randn(N_sample, 1) / sqrt(2); % 白噪聲的虛部v = v_real + 1i*v_imag; % 復(fù)白噪聲n = 1 : N_sample;n = n's = exp(1i * 2 * pi * 0.5 * n) + exp(1i * (2 * pi * 0.55 * n + pi/4); x = s + v; % 觀察信號% PHD算法實現(xiàn)%根據(jù)信號取樣求出自相關(guān)序列，進而求出3階相關(guān)矩陣N = 3; % 題(4)中自相關(guān)

24、矩陣Rx的階數(shù)x_corr,a = xcorr(x,'unbiased'); %求相關(guān)函數(shù)無偏估計Rx = (toeplitz(x_corr(N_sample : N_sample + N -1).' %求自相關(guān)矩陣% 求Rx的特征值和特征向量，并找出最小特征值對應(yīng)的特征向量V,D = eig(Rx);% 求特征值和特征向量vmin = V(:,1); % 最小特征值對應(yīng)的特征向量N=M+1% 解正交方程4.154，直接解出頻率值Z = roots(vmin.'); % 方程(4.152)求解，得到w = angle(Z); % 零點的角度for i = 1 : M if w(i) < 0 w(i) = w(i) + 2 * pi% angle函數(shù)求出的角度范圍是-pi, pi),當求出的角頻 %率w <時，需要加上一個2*pi以調(diào)整到0, 2*pi)之間 endendf_PHD = w / (2*pi)% MUSCI算法by劉航20151109% 初始化clc;clear all;N = 12; % 題(5)中自相關(guān)矩陣Rx的階數(shù)N_sample = 25;% 取樣自相關(guān)點數(shù)M = 2; % 正弦信號的個數(shù)% 產(chǎn)生復(fù)白