有代碼功率譜估計Levinson遞推法Burg遞推法隨機(jī)信號處理

1、隨機(jī)信號處理實驗報告功率譜估計隨機(jī)信號處理 學(xué)號： 姓名： 實驗三 功率譜估計1實驗內(nèi)容信號為兩個正弦信號加高斯白噪聲，各正弦信號的信噪比均為10dB，長度為N，信號頻率分別為和，初始相位，取，取不同數(shù)值：0.3，0.25。為采樣頻率。分別用Levinson遞推法和Burg法進(jìn)行功率譜估計，并分析改變數(shù)據(jù)長度、模型階數(shù)對譜估計結(jié)果的影響。2實驗原理2.1 Levinson遞推法：自相關(guān)法列文森（Lenvison）遞推法是已知信號觀測數(shù)據(jù)，估計功率譜。它的出發(fā)點是選擇AR模型參數(shù)使預(yù)測誤差功率最小。假設(shè)信號的數(shù)據(jù)區(qū)在范圍，有P個預(yù)測系數(shù)，N個數(shù)據(jù)經(jīng)過沖激響應(yīng)為的濾波器，輸出預(yù)測誤差的長度為，因此

2、有預(yù)測誤差功率為的長度長于數(shù)據(jù)的長度，上式中數(shù)據(jù)在以外補充零點，相當(dāng)于對無窮長的信號加窗處理，會引入誤差。上式對系數(shù)的實部和虛部求微分使預(yù)測誤差功率最小，得（Yule-Walker方程）式中自相關(guān)函數(shù)采用有偏自相關(guān)估計，即Levinson-Durbin算法：使一種按階次遞推的算法。它以和模型參數(shù)作為初始條件，計算模型參數(shù)；再用模型參數(shù)計算模型參數(shù)，k階模型參數(shù)由k-1階模型參數(shù)計算得到。一直計算出模型參數(shù)為止。一階AR模型的Yule-Walker方程為由該方程解出然后令，以此類推，可以得到一般遞推公式如下：稱為反射系數(shù)，。，隨著階數(shù)增加，預(yù)測誤差功率將減少或不變。由k=1開始遞推，遞推到k=p

3、，依次得到各階模型參數(shù)，AR模型的各個系數(shù)及模型輸入白噪聲方差求出后，信號功率譜用下式計算這種方法遞推效率高，當(dāng)階數(shù)變化時，無需從頭計算。但需要預(yù)先估計出信號自相關(guān)函數(shù)，當(dāng)觀測數(shù)據(jù)長度較短時，估計誤差較大，會出現(xiàn)譜峰頻率偏移和譜線分裂；如數(shù)據(jù)很長，估計自相關(guān)函數(shù)較準(zhǔn)確。2.2 Burg遞推法：Levinson-Durbin遞推法需要由觀測數(shù)據(jù)估計自相關(guān)函數(shù)，這是它的缺點。而伯格遞推法則由信號觀測數(shù)據(jù)直接計算AR模型參數(shù)。伯格遞推法利用Levinson-Durbin遞推公式，導(dǎo)出前向預(yù)測誤差與后向預(yù)測誤差，并按照使它們最小的原則求出，從而實現(xiàn)不用估計自相關(guān)函數(shù)，直接用觀測數(shù)據(jù)得出結(jié)果。Burg遞

4、推法思想：借助格型預(yù)測誤差濾波器，求前向、后向預(yù)測誤差平均功率，選擇使其最小，求出。之后，再利用Levinson-Durbin遞推法求模型參數(shù)和輸入噪聲方差。設(shè)信號的觀測數(shù)據(jù)區(qū)間：,前向、后向預(yù)測誤差功率分別用和表示，預(yù)測誤差平均功率用表示，公式分別為前向、后向觀測誤差公式分別為 上式中，信號項的自變量最大的是n,最小的是n-p，為了保證計算范圍不超出給定的數(shù)據(jù)范圍，在和計算公式中，選擇求和范圍為： 。為求預(yù)測誤差平均功率最小時的反射系數(shù)，令,將前、后向預(yù)測誤差的遞推公式代入得Burg遞推法求AR模型參數(shù)的遞推公式總結(jié)：(1) (2) (3) (4) (5) (6) (7) 3實驗結(jié)果及分析3

5、.1原始信號，觀測信號這里取，。3.2 Levenson遞推法3.2.1 取，或，階數(shù)不變，實驗不同數(shù)據(jù)長度對功率譜估計的影響1) 信號長度信號長度N=35，階數(shù)M=20的功率譜估計2) 信號長度 信號長度N=145，階數(shù)M=20的功率譜估計3) 信號長度 信號長度N=2000，階數(shù)M=20的功率譜估計分析：由以上三個實驗對比，可以看出當(dāng)觀測數(shù)據(jù)長度較短時，估計誤差較大，會出現(xiàn)譜峰頻率偏移與譜線分裂；當(dāng)數(shù)據(jù)很長時，估計自相關(guān)函數(shù)較準(zhǔn)確，但計算量較大。3.2.2取，信號長度不變，實驗不同模型階數(shù)對功率譜估計的影響1）階數(shù)M=22）階數(shù)M=43）階數(shù)M=8 4)階數(shù)M=16分析：由以上幾個實驗對比

6、，可以看出當(dāng)階次較低，會使譜估計產(chǎn)生偏移，降低分辨率；當(dāng)階次越高，分辨率越高；當(dāng)階次太高，會使估計誤差加大，譜峰分裂。3.3 Burg遞推法3.3.1取，或，階數(shù)不變，實驗不同數(shù)據(jù)長度對功率譜估計的影響1) 信號長度 信號長度N=35，階數(shù)M=20的功率譜估計2) 信號長度 信號長度N=145，階數(shù)M=20的功率譜估計3) 信號長度 信號長度N=2000，階數(shù)M=20的功率譜估計分析：由以上三個實驗對比，可以看出當(dāng)觀測數(shù)據(jù)長度較短時，估計誤差較大，會出現(xiàn)譜峰頻率偏移與譜線分裂；當(dāng)數(shù)據(jù)很長時，估計自相關(guān)函數(shù)較準(zhǔn)確，但計算量較大。頻率越靠近的譜估計，需要的階數(shù)越高。3.3.2 取，信號長度不變，實

7、驗不同模型階數(shù)對功率譜估計的影響1) 階數(shù)M=42) 階數(shù)M=83) 階數(shù)M=164) 階數(shù)M=28分析：由以上幾個實驗對比，可以看出當(dāng)階次較低，會使譜估計產(chǎn)生偏移，降低分辨率；當(dāng)階次越高，分辨率越高；當(dāng)階次太高，會使估計誤差加大，譜峰分裂。3.4實驗總結(jié)本次試驗采用分別用Levinson遞推法和Burg遞推法進(jìn)行功率譜估計，并分析改變數(shù)據(jù)長度、模型階數(shù)對譜估計結(jié)果的影響。通過實驗，學(xué)習(xí)了Levinson遞推法和Burg遞推法的基本原理和一般流程，和如何選擇AR模型的階次，并使用Matlab語言，編寫源代碼，完成實驗過程。在實驗過程中，分別設(shè)計了不同信號長度、不同的AR模型的階次和不同頻率組合

8、而成的4組實驗，并在實驗后，分析對比了實驗結(jié)果。5源代碼5.1 Levenson遞推法% clear all;Clear；tic;% 產(chǎn)生信號fs=1;%設(shè)采樣頻率為1N=100;%數(shù)據(jù)長度 改變數(shù)據(jù)長度會導(dǎo)致分辨率的變化f1=0.2*fs;%第一個sin信號的頻率,f1/fs=0.2f2=0.3*fs;%第二個sin信號的頻率，f1/fs=0.2或者0.3M = 60;%濾波器階數(shù)的最大取值,超過則認(rèn)為代價太大而放棄L = 2*N;%有限長序列進(jìn)行離散傅里葉變換前，序列補零的長度n=1:N; s = sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s為原始信號x =

9、 awgn(s,10);%x為觀測信號，即對原始信號加入白噪聲，信噪比10dB%畫出原始信號和觀測信號figure(1);subplot(2,1,1);plot(s,b),xlabel(時間),ylabel(幅度),title(原始信號s);grid;subplot(2,1,2);plot(x,r),xlabel(時間),ylabel(幅度),title(觀測信號x);grid;%計算自相關(guān)函數(shù)rxx = xcorr(x,x,M,biased);%計算有偏估計自相關(guān)函數(shù)，長度為-M到M，共%2M+1r0 = rxx(M+1); %r0為零點上的自相關(guān)函數(shù)，相對于-M，第M+1個點為零點R =

10、rxx(M+2:2*M+1);% R為從1到第M個點的自相關(guān)函數(shù)矩陣%Levinson遞推算法%確定矩陣大小a = zeros(M,M);FPE = zeros(1,M);%FPE：最終預(yù)測誤差，用來估計模型的階次var = zeros(1,M);%求初值a(1,1) = -R(1)/r0;%一階模型參數(shù)var(1) = (1-(abs(a(1,1)2)*r0;%一階方差FPE(1) = var(1)*(M+2)/(M);%遞推for p=2:M sum=0; for k=1:p-1%求a(p,p) sum=sum+a(p-1,k)*R(p-k); end a(p,p)=-(R(p)+sum)

11、/var(p-1); for k=1:p-1 %求a(p,k) a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k); end var(p)=(1-a(p,p)2)*var(p-1); %求方差 FPE(p)=var(p)*(M+1+p)/(M+1-p);%求最終預(yù)測誤差end %確定AR模型的最佳階數(shù)min=FPE(1); %求出FPE最小時對應(yīng)的階數(shù)p = 1;for k=2:M if FPE(k)2 for i=1:p-2 a(p-1,i)=a(p-2,i)+k(p-1)*a(p-2,p-1-i); end end a(p-1,p-1)=k(p-1);% 求解前向預(yù)測誤差 for n=p+1:N ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1); end% 求解后向預(yù)測誤差 for n=p:N-1 eb(p,n)=eb(p-1,n-1)+k(p-1)*ef(p-1,n); endend % 計算功率譜for j=1:N sum3=0; sum4=0; for i=1:p-1 sum3=sum3+a(p-1,i)*cos(2*pi*i*j/N); end sum3=1+sum3; for i=1:p-1 sum4=sum4+a(p-1,i)*sin(2*pi*i*j/N);