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

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

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

4、式計算:NF'%ennz0e(n)的長度長于數據的長度，上式中數據x(n)的兩端需補充零點，相當于對無窮長的信號加窗處理，得到長度為N的數據。上式對系數api的實部和虛部求微分使預測誤差功率最小,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)式中自相關函數采用有偏自相關估計，即i1N-4-ma-Zx*(nX(n+m)m=0,1,2,，pixx(m)=Nn田a*ixx(-m)m=_p+1,_p+2J',_1對比上式，可知式(1)即為已推導出的Y

5、ule-Walker方程，因此自相關法也是基于解Yule-Walker方程的一種方法。但是直接解該方程，需要計算逆矩陣，不方便，因此，基于Yule-Walker方程中自相關矩陣的性質，導出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-=(l-)cr:F、F*

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

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

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

9、AR莫型參數的遞推公式總結如下: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限ap,i=a；,ipapAp-Li=1,2,p-1kp=ap,pepn=epnkpe"n-1n=p1,p2,N-1epn)=e"n一1kjepnn=p1,p2,N一13.編程思想(1)編寫程序產生題目要求的信號和噪聲(2)然后分別用兩種方法的遞推流程進行譜估計(3)改變題目中要求的變量參數，分析結果的變化4.代碼Levensonclc;cleara

10、ll;fs=100;%采樣頻率Ts=1/fs;N=2A7;%數據長度p1=20;%數f1=0.2*fs;f2=0.25*fs;%設置信號頻率pha1=0;pha2=0;%J始相位SNR=2;%設置信噪比%產生信號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('兩個正弦信號波形);gridon;su

11、bplot(2,1,2);plot(x);ylim(-20,20);title('加噪信號波形);gridon;%計算自相關函數R=;form=1:Ns=0;forn=1:N-(m-1)s=s+x(n)*x(n+m-1);endR(m)=s/N;end發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);%#算反射系數fork=1:p-1a(p,k)=a(p-1,k)+a(p,

12、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);%角頻率坐標換算成頻率坐標figure(2)plot(F,abs(PSD);xlabel('頻率(Hz)');ylabel('功率譜);title('自相關法-列文森Levenson遞推法的

13、功率譜估計);grid;figure(3)p=1:p1;plot(p,cov(p),'b');xlabel('模型階數');ylabel('預測誤差功率大小);title('預測誤差功率變化趨勢');grid;Burgclc;clearall;fs=100;%采樣頻率Ts=1/fs;N=2A8;煨據長度p1=20;%數f1=0.2*fs;f2=0.25*fs;%設置信號頻率pha1=0;pha2=0;SNR=2;婕置信噪比%產生信號w=randn(1,N);%噪聲為高斯白噪聲，功率為1.Am=sqrt(2*10A(SNR/10);k=1:

14、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=denominator+(ef(p-1,n)A2+(eb(

15、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=1:0.01:pi;fori=1:pZ

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

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

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

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

20、20253035404550頻率(Hz)Pi/6自相關法-列文森Levenson遞推法的功率譜估計450400350300加250至功20015010050005101520253035404550頻率(Hz)Pi/31/450400350300250200150100500自相關法-列文森Levenson遞推法的功率譜估計05101520253035404550頻率(Hz)Pi/23503002502001501005010015101520253035404550400自相關法-列文森Levenson遞推法的功率譜估計頻率(Hz)Pi頻率（Hz）4.頻率的影響(N=128p1=20初始相位

21、為0)Fs=100,f1=0.2*fsf2=0.25*fs頻率(Hz)Fs=100,f1=0.2*fsf2=0.3*fs400350300250200150100501105101520253035404550自相關法-列文森Levenson遞推法的功率譜估計頻率(Hz)Fs=1000,f1=0.2*fsf2=0.25*fs4002001J50100150200250300350400450500頻率（Hz）50010000Fs=1000,f1=0.2*fsf2=0.3*fs自相關法-列文森Levenson遞推法的功率譜估計頻率（Hz）5.信噪比的影響SNR=20dB12001000逆800&

22、#165;功八60040020000111JI5101520253035404550頻率(Hz)SNR=10dB600500自相關法-列文森Levenson遞推法的功率譜估計400埴率300面200100J1/11051015202530頻率(Hz)35404550SNR=6dB300自相關法-列文森Levenson遞推法的功率譜估計250200霆150場100501111I/10510152025頻率(Hz)3035404550SNR=2dB80604020111115101520253035404550頻率(Hz)5.3Burg遞推法200100015I/20253035404550Fre

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

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