功率譜估計Levinson遞推法和Burg法

1、數(shù)字信號處理實驗報告姓名：學(xué)號：日期：2015.12.141 .實驗任務(wù)信號為兩個正弦信號加高斯白噪聲，各正弦信號的信噪比均為10dB,長度為N，信號頻率分別為先和f2,初始相位%=52=。，取fjfs=0.2,fjfs取不同的數(shù)值：0.3,0.25。fs為采樣率。(1)分別用Levinson遞推法和Burg法進行功率譜估計，并分析改變數(shù)據(jù)長度、模型階數(shù)對譜估計結(jié)果的影響。(2)當(dāng)正弦信號相位、頻率、信噪比改變后，上述譜估計的結(jié)果有何變化？并作分析說明。2 .原理分析2.1 現(xiàn)代譜估計中的參數(shù)建模根據(jù)參數(shù)模型來描述隨機信號的方法，我們可以知道，如果能確定信號x(n)的信號模型，根據(jù)信號觀測數(shù)據(jù)

2、求出模型參數(shù)，系統(tǒng)函數(shù)用H(z底示，模型輸入白噪聲，其方差為2仃w，信號的功率譜用下式求出：Pxxfejw)=Ow2He按照這種求功率譜的思路，功率譜估計可分為三個步驟：(1)選擇合適的信號模型；(2)根據(jù)x(n府限的觀測數(shù)據(jù)，或者它的有限個自相關(guān)函數(shù)的估計值，估計模型的參數(shù)；(3)計算墨香的輸出功率譜。其中以(1)、(2)兩步最為關(guān)鍵。按照模型的不同，譜估計的方法有許多種，它們共同的特點是對信號觀測區(qū)以外的數(shù)據(jù)不假設(shè)為0,而先根據(jù)信號觀測數(shù)據(jù)估計模型參數(shù)，按照求模型輸出功率的方法估計信號功率譜，回避了數(shù)據(jù)觀測區(qū)以外的數(shù)據(jù)假設(shè)問題。下面分析AR譜估計的兩種方法：自相關(guān)法列文森(Levenson

3、)遞推法和伯格(Burg)遞推法。這兩種方法均為已知信號觀測數(shù)據(jù)，估計功率譜，兩者共同特點是由信號觀測數(shù)據(jù)求模型系數(shù)時采用信號預(yù)測誤差最小的原則。對于長記錄數(shù)據(jù)，這些方法的估計質(zhì)量是相似的，但對于短記錄數(shù)據(jù)，不同方法之間存在差別。2.2 自相關(guān)法列文森(Levenson遞推法自相關(guān)法的出發(fā)點是選擇AR模型參數(shù)使預(yù)測誤差功率最小，預(yù)測誤差功率為一2一18,1gee(n|=一£NNnRp2xnJapiXn-i=1假設(shè)彳t號x(n)的數(shù)據(jù)區(qū)在0MnMN1范圍，有P個預(yù)測系數(shù)，N個數(shù)據(jù)經(jīng)過沖激響應(yīng)NZn=0pXn-二apiXn-ii1為api(i=0,1，p)的濾波器，輸出預(yù)測誤差e(n)的

4、長度為N+p,因此應(yīng)用下式計算:NF'%ennz0e(n)的長度長于數(shù)據(jù)的長度，上式中數(shù)據(jù)x(n)的兩端需補充零點，相當(dāng)于對無窮長的信號加窗處理，得到長度為N的數(shù)據(jù)。上式對系數(shù)api的實部和虛部求微分使預(yù)測誤差功率最小,p得到一朦(0)?x(-1)I朦(0)mm,0(p-1)l?x(p2)i?x(-p+1frap1&p+2)ap2iW)1app_一？(1)1"(2)rXx(p)(1)式中自相關(guān)函數(shù)采用有偏自相關(guān)估計，即i1N-4-ma-Zx*(nX(n+m)m=0,1,2,，pixx(m)=Nn田a*ixx(-m)m=_p+1,_p+2J',_1對比上式，可知

5、式(1)即為已推導(dǎo)出的Yule-Walker方程，因此自相關(guān)法也是基于解Yule-Walker方程的一種方法。但是直接解該方程，需要計算逆矩陣，不方便，因此，基于Yule-Walker方程中自相關(guān)矩陣的性質(zhì)，導(dǎo)出Levinson-Durbin遞推法，這是一種高效的解方程的方法。Levinson-Durbin算法首先由一階AR模型開始：一階AR模型(p=1)的Yule-Walker方程為rxx01rxx1rxx由該方程解出rxx1a)i=-rxx02/2-1-1_a1,1Irxx0.1然后令p=2,3,4,以此類推，可以得到一般遞推公式如下：r“二一<%*=總卜口41左=L2.，夕7cr-

6、=(l-)cr:F、F*p-1口;二:(。)=修了5)222kp稱為反射系數(shù)，kp<1o。1皂。2之之bp,隨著階數(shù)增加，預(yù)測誤差功率將減少或不變。由k=1開始遞推，遞推到k=p,依次得到各階模型參數(shù)，',a11,01J,1a21,a22,二2J,1ap1,ap2,app,二pJAR模型的各個系數(shù)及模型輸入白噪聲方差求出后，信號功率譜用下式計算2(pfpxx(ejw)=。：H(ejwJ=。"/1+£ape，1yJ這種方法計算簡單，但需要預(yù)先估計出信號自相關(guān)函數(shù)，實際中只能按照信號的有限個觀測數(shù)據(jù)估計自相關(guān)函數(shù)。當(dāng)觀測數(shù)據(jù)長度較短時，估計誤差較大，會出現(xiàn)譜峰頻率

7、偏移和譜線分裂(在信號譜峰附近產(chǎn)生虛假譜線)；如數(shù)據(jù)很長，估計自相關(guān)函數(shù)較準(zhǔn)確，但計算量大，應(yīng)適當(dāng)選擇數(shù)據(jù)長度。2.3伯格(Burg)遞推法Levinson-Durbin遞推法需要由觀測數(shù)據(jù)估計自相關(guān)函數(shù)，這是它的缺點。而伯格遞推法則由信號觀測數(shù)據(jù)直接計算AR模型參數(shù)。伯格遞推法利用Levinson-Durbin遞推公式，導(dǎo)出前向預(yù)測誤差與后向預(yù)測誤差，并按照使它們最小的原則求出k,從而實現(xiàn)不用估計自相關(guān)函數(shù)，直接用觀測數(shù)據(jù)得出結(jié)果。pBurg遞推法思想：借助格型預(yù)測誤差濾波器，求前向、后向預(yù)測誤差平均功率，選才ikpp使其最小，求出kp。之后，再利用Levinson-Durbin遞推法求模型

8、參數(shù)和輸入噪聲方差。設(shè)信號x(n)的觀測數(shù)據(jù)區(qū)間：0wnwN1,前向、后向預(yù)測誤差功率分別用Ppf和Ppb表示，預(yù)測誤差平均功率用Pp表示，公式分別為np前向、后向觀測誤差公式分別為epnp十ppfpbfJ.epn)=xn工apkxn-kk1pebpn=xn-pVa*pkxn-pkk1ap0=1上式中，信號項的自變量最大的是n,最小的是n-p,為了保證計算范圍不超出給定的數(shù)據(jù)范圍，在Ppf和Ppb計算公式中，選擇求和范圍為：p<n<N-1odP為求預(yù)測誤差平均功率Pp最小時的反射系數(shù)kp,令一p=0,將前、后向預(yù)測誤差的遞二kp推公式代入得kp=N4zn=pN1-2、efnepn-

9、1n-pBurg遞推法求AR莫型參數(shù)的遞推公式總結(jié)如下:crxx1N2小(0戶工x(n),P(0尸rxx(0)n=0(2)fn)=x(n)n=0,1,2,，N-16b(n)=x(n)n=0,1,2,，N-1N_1-2%e；ne；*n-1Pp=1_kp2限(5) ap,i=a；,ipapAp-Li=1,2,p-1(6) kp=ap,pepn=epnkpe"n-1n=p1,p2,N-1epn)=e"n一1kjepnn=p1,p2,N一13.編程思想(1)編寫程序產(chǎn)生題目要求的信號和噪聲(2)然后分別用兩種方法的遞推流程進行譜估計(3)改變題目中要求的變量參數(shù)，分析結(jié)果的變化4.

10、代碼Levensonclc;clearall;fs=100;%采樣頻率Ts=1/fs;N=2A7;%數(shù)據(jù)長度p1=20;%數(shù)f1=0.2*fs;f2=0.25*fs;%設(shè)置信號頻率pha1=0;pha2=0;%J始相位SNR=2;%設(shè)置信噪比%產(chǎn)生信號w=randn(1,N);Am=sqrt(2*10A(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);xx=x1+x2;x=w+x1+x2;figure(1)subplot(2,1,1);plot(xx);ylim(-20,20);title('

11、;兩個正弦信號波形);gridon;subplot(2,1,2);plot(x);ylim(-20,20);title('加噪信號波形);gridon;%計算自相關(guān)函數(shù)R=;form=1:Ns=0;forn=1:N-(m-1)s=s+x(n)*x(n+m-1);endR(m)=s/N;end發(fā)U文森遞推法a(1,1)=-R(2)/R(1);cov(1)=R(1)+a(1,1)*R(2);forp=2:p1sum=0;fork=1:p-1sum=sum+a(p-1,k)*R(p-k+1);enda(p,p)=-(R(p+1)+sum)/cov(p-1);%#算反射系數(shù)fork=1:p-1

12、a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);endcov(p)=(1-(abs(a(k,k)A2)*cov(p-1);end%計算功率譜，功率譜以2*pi為周期，又信號為實信號，只需輸出0到P1即可；W=0.01:0.01:pi;Z=0;fork=1:p1Z=Z+(a(p1,k).*exp(-j*k*W);endPSD=1./(abs(1+Z)A2);%?出功率譜F=W*fs/(pi*2);%角頻率坐標(biāo)換算成頻率坐標(biāo)figure(2)plot(F,abs(PSD);xlabel('頻率(Hz)');ylabel('功率譜);title('

13、自相關(guān)法-列文森Levenson遞推法的功率譜估計);grid;figure(3)p=1:p1;plot(p,cov(p),'b');xlabel('模型階數(shù)');ylabel('預(yù)測誤差功率大小);title('預(yù)測誤差功率變化趨勢');grid;Burgclc;clearall;fs=100;%采樣頻率Ts=1/fs;N=2A8;煨據(jù)長度p1=20;%數(shù)f1=0.2*fs;f2=0.25*fs;%設(shè)置信號頻率pha1=0;pha2=0;SNR=2;婕置信噪比%產(chǎn)生信號w=randn(1,N);%噪聲為高斯白噪聲，功率為1.Am=sqr

14、t(2*10A(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);x=w+x1+x2;%遞推ef=zeros(p1,N);eb=zeros(p1,N);ef(1,：)=x;eb(1,:)=x;cov(1)=x*x'/N;k(1)=0;a=zeros(p1+1,p1+1);forp=2:p1+1numerator=0;denominator=0;forn=p:Nnumerator=numerator+(-2)*ef(p-1,n)*eb(p-1,n-1);denominator=denominat

15、or+(ef(p-1,n)A2+(eb(p-1,n-1)A2;endk(p)=numerator./(denominator+0.0001);cov(p)=(1-abs(k(p)A2).*cov(p-1);a(p,p)=k(p);forn=p:Nef(p,n)=ef(p-1,n)+k(p).*eb(p-1,n-1);eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n);endendk=k(1,2:p1+1);a=a(2:p1+1,2:p1+1);forp=2:p1fori=1:p-1a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i);endend%功率譜Z=0;W

16、=1:0.01:pi;fori=1:pZ=Z+(a(p,i)*exp(-j*W*i);end;pxx=cov(p1+1)./(abs(1+Z).A2);F=W*fs/(pi*2);%角頻率坐標(biāo)換算成頻率坐標(biāo)figure(1)plot(F,pxx);xlabel('Frequency(Hz),);ylabel('PowerSpectralDensity');title('伯格(Burg)遞推法的功率譜估計);grid;%figure(2)%p=1:p1;%plot(p,cov(p),'b');%xlabel('模型階次');%yla

17、bel('預(yù)測誤差功率大小');%title('預(yù)測誤差功率的變化趨勢');%grid;5.實驗結(jié)果及分析5.1產(chǎn)生信號兩個正弦信號波形加噪信號波形5.2Levensond的推法1.數(shù)據(jù)長度的影響（階數(shù)設(shè)置為20階）N=32N=64160140120港100¥功8060402000111-j1,5101520253035404550頻率（Hz）N=128600自相關(guān)法-列文森Levenson遞推法的功率譜估計500400-rife組率300功2001000011J115101520253035404550頻率（Hz）N=1024300025002000

18、率1500功100050000自相關(guān)法-列文森Levenson遞推法的功率譜估計,jL5101520253035404550頻率（Hz）N=8192600050004000埴率3000功2000100025頻率（Hz）N=163846000自相關(guān)法-列文森Levenson遞推法的功率譜估計5000400030002000100010051015202530頻率（Hz）354045502,模型階數(shù)的影響（數(shù)據(jù)長度設(shè)置為N=1024）預(yù)測誤差功率的變化（從1階到50階）預(yù)測誤差功率變化趨勢20太15I105005101520253035404550模型階數(shù)不同階數(shù)（pl表示）功率譜的圖像如下:P1

19、=570自相關(guān)法-列文森Levenson遞推法的功率譜估計6050JI、W4011t30201000J1510152025頻率30(Hz)35404550P1=10300250自相關(guān)法-列文森Levenson遞推法的功率譜估計200缶150功100I00510152025頻率3035404550(Hz)P1=201000用舉功500005101520253035404550頻率（Hz）P1=30自相關(guān)法-列文森Levenson遞推法的功率譜估計頻率（Hz）P1=50頻率（Hz）3.相位的影響（設(shè)置N=128p1=20）0自相關(guān)法-列文森Levenson遞推法的功率譜估計600500400遭率3

20、00面2001000:J1105101520253035404550頻率(Hz)Pi/6自相關(guān)法-列文森Levenson遞推法的功率譜估計450400350300加250至功20015010050005101520253035404550頻率(Hz)Pi/31/450400350300250200150100500自相關(guān)法-列文森Levenson遞推法的功率譜估計05101520253035404550頻率(Hz)Pi/23503002502001501005010015101520253035404550400自相關(guān)法-列文森Levenson遞推法的功率譜估計頻率(Hz)Pi頻率（Hz）4.

21、頻率的影響(N=128p1=20初始相位為0)Fs=100,f1=0.2*fsf2=0.25*fs頻率(Hz)Fs=100,f1=0.2*fsf2=0.3*fs400350300250200150100501105101520253035404550自相關(guān)法-列文森Levenson遞推法的功率譜估計頻率(Hz)Fs=1000,f1=0.2*fsf2=0.25*fs4002001J50100150200250300350400450500頻率（Hz）50010000Fs=1000,f1=0.2*fsf2=0.3*fs自相關(guān)法-列文森Levenson遞推法的功率譜估計頻率（Hz）5.信噪比的影響S

22、NR=20dB12001000逆800¥功八60040020000111JI5101520253035404550頻率(Hz)SNR=10dB600500自相關(guān)法-列文森Levenson遞推法的功率譜估計400埴率300面200100J1/11051015202530頻率(Hz)35404550SNR=6dB300自相關(guān)法-列文森Levenson遞推法的功率譜估計250200霆150場100501111I/10510152025頻率(Hz)3035404550SNR=2dB80604020111115101520253035404550頻率(Hz)5.3Burg遞推法20010001

23、5I/20253035404550Frequency(Hz)N=64500伯格(Burg)遞推法的功率譜估計450400y350D300De250W2001501100500151J20253035404550Frequency(Hz)N=2567000151111600500400300200100202530354045伯格（Burg）遞推法的功率譜估計50Frequency(Hz)N=102412000100008000600040002000015伯格（Burg）遞推法的功率譜估計2025403035Frequency(Hz)4550N=819245004000350030002500

24、200015001000500015202540伯格（Burg）遞推法的功率譜估計453035Frequency(Hz)50N=1638460005000400030002000100001530Frequency(Hz)2 .模型階數(shù)白影響（N=128）模型階次不同階數(shù)的功率譜015/i1ii1/f20253035404550Frequency(Hz)10QO521P1=10計估譜率功的法推遞UB格伯101010'0'0'0'0ooQQQoO876543201512025404550P1=20伯格(Burg)遞推法的功率譜估計400350300250200150100503035Frequency(Hz)15010050Frequency(Hz)3 .相位的影響(N=256p1=20)0Pi/6Frequency(Hz)5O10000008642yrncouaJpHJBO-Frequency(Hz)5o1oooooooooooooo7654321jjy/mDnuapJgasL/3piFrequency(Hz)5o