報告2AR過程的線性建模與功率譜估計

1、實驗報告2AR過程的線性建模與功率譜估計一、實驗目的1理解AR過程的產(chǎn)生機理，復習實驗1估計自相關(guān)序列的方法。2利用估計出的自相關(guān)序列來求解信號的功率譜，即用周期圖法來估計功率譜。分別采用自相關(guān)法(Yule-Walker法)，協(xié)方差法，Burg法，修正協(xié)方差法來估計功率譜，并與周期圖法進行比較，分析性能孰優(yōu)孰劣。學習matlab在數(shù)字信號處理中的應用。二、實驗過程和分析AR過程的線性建模與功率譜估計?？紤]AR過程：x(n)=a(1)x(n-1)+a(2)x(n-2)+a(3)x(n-3)+a(4)x(n-4)+b(O)v(n)v(n)是單位方差白噪聲。(a)取b(0)=1,a(1)=2.760

2、7,a(2)=-3.8106,a(3)=2.6535,a(4)=-0.9238，產(chǎn)生x(n)的N=256個樣點。用randn(1,N)產(chǎn)生單位方差高斯白噪聲v(n)，用v(n)激勵濾波器產(chǎn)生AR(4)過程，即用x=filter(b,a,v)產(chǎn)生x(n)，b是濾波器分子系數(shù)，這里為b(0)=1，a是濾波器分母系數(shù)，a=1-2.76073.8106-2.65350.9238。計算其自相關(guān)序列的估計r(k)，并與真實的自相關(guān)序列值相比較。x15001000111111111111111:111111111111來估#自相關(guān)真實自相關(guān)11自相關(guān)估計5000-500-1000-1500501001502

3、00250300結(jié)論：真實自相關(guān)序列與估計出的序列對比如上圖所示，兩者在100點之前的形狀相似,幅度有一定差異，而且估計出的自相關(guān)序列有較大的波動，這是因為估計的點數(shù)較少，使得估計精度不夠，另外，估計自相關(guān)序列的下標越大，用來估計的點數(shù)就越少，因而后面的估計值是很不精確的。因為假設(shè)n=500處的數(shù)據(jù)為0,所有估計自相關(guān)不精確,誤差較大。(c)將下(k)的DTFT作為x(n)的功率譜估計，即:xP(ej)=r(k)e-jk=xx丄丨X(ej)|2NOO周期圖法獲取的功率譜與真實功率譜的比較0.10.20.30.40.50.60.70.80.9數(shù)字頻率怎ooO543OOO-k=-N+1結(jié)論：周期圖

4、法求得的功率譜誤差較大，波動也較大，兩個譜峰的分辨也并不清晰，并不精確，這是由于周期圖法對自相關(guān)序列加了一個矩形窗。(d)利用所估計的自相關(guān)值和Yule-Walker法(自相關(guān)法)，估計a(1),a(2),a(3),a和b(0)的值，并討論估計的精度。a(l)=2.5018a(2)=-3.2213a(3)=2.0794a(4)=-0.7025b(0)=2.4198對比實際結(jié)果，可知Yule-Walker法(自相關(guān)法)誤差較大，估計并不理想。(e)用(d)中所估計的a(k)和b(0)來估計功率譜為:P(e.)=1b(0)|2x1+才a(k)e-jkk=iCQP、S3祥o利用Yule-Walker

5、te估計功率譜0.10.20.30.40.50.60.70.80.91數(shù)字頻率ooo543o20oJI-估計出的功率譜上圖所示，難以看出有二個譜峰，這是由于自相關(guān)法在短數(shù)據(jù)量的情況下估計是不精確的。將(c)和(e)的兩種功率譜估計與實際的功率譜進行比較，畫出它們的重疊波形。OO功率譜比較0.10.20.30.40.50.60.70.80.9數(shù)字頻率ooO432OOOJ-結(jié)論：用Yule-Walker法估計得到的譜曲線光滑，但無法識別出有二個譜峰，且b(0)的估計明顯是錯的，使得幅度有較大的偏差；相關(guān)圖的曲線不光滑，不好觀察譜峰，但大致走向與真實譜是一致的。可以發(fā)現(xiàn)，Yule-Walker法(參

6、數(shù)法)比周期圖法(非參數(shù)法)精度更高。但這里Yule-Walker法得到的功率譜分辨率較低，這是由于對自相關(guān)序列加了一個矩形窗引起的。(g)重復上面的(d)(f),只是估計AR參數(shù)分別米用如下方法：(1)協(xié)方差法；(2)Burg方法；(3)修正協(xié)方差法。試比較它們的功率譜估計精度。協(xié)方差法a(1)=2.6843a(2)=-3.6442a(3)=2.5026a(4)=-0.8668b(0)=1.5506mpEn芯徂aAWQO功率譜比較0.10.20.30.40.50.60.70.80.9數(shù)字頻率與Yule-Walker法估計出參數(shù)的誤差相比，有非常大的改進，這是因為協(xié)方差法對自相關(guān)序列估計時不需

7、要對數(shù)據(jù)加窗，因此，對短數(shù)據(jù)記錄，協(xié)方差法一般可獲得比自相關(guān)法更高的分辨率。由上圖可知，譜估計的幅度相差較大，但從形狀和極點定位的角度而言,估計結(jié)果很好的，曲線光滑，且極點的定位非常準確，短觀測數(shù)據(jù)也能對極點精確定位是協(xié)方差法的優(yōu)點。Burg方法a(l)=2.7872a(2)=-3.8728a(3)=2.7115a(4)=-0.9467b(0)=1.3006O功率譜比較0.10.20.30.40.50.60.70.80.9數(shù)字頻率Burg算法最小化前向加反向預測誤差的平方和來求解全極點參數(shù)，為了保證穩(wěn)定性，這一最小化是相對于反射系數(shù)序貫地進行的，其估計精度低于協(xié)方差法。此估計方法與Yule-W

8、alker法估計出參數(shù)的誤差相比，有非常大的改進，但與協(xié)方差法比較，精度沒有其理想。用估計出的參數(shù)來進行譜估計，如上圖所示，藍色虛線就是采用Burg法估計出的功率譜，紅線是真實功率譜，由圖可知，估計結(jié)果較好，曲線光滑，且極點的定位較為準確，Burg法是自相關(guān)法與協(xié)方差法的一個折中。修正協(xié)方差法a(l)=2.7863a(2)=-3.8716a(3)=2.7106a(4)=-0.9467b(0)=1.3006功率譜比較0.10.2030.40.50.60.70.80.9數(shù)字頻率九修正協(xié)方差法是最小化前向加反射預測誤差的平方和，峰值偏移對相位的敏感度較低，比Yule-Walker法精度高，與Burg

9、法的精度相當，但比協(xié)方差法的精度要低一些。用估計出的參數(shù)來進行譜估計，如上圖所示，藍色虛線就是采用修正協(xié)方差法估計出的功率譜，紅線是真實功率譜，由圖可知，估計結(jié)果較好，曲線光滑，且極點的定位較為準確。結(jié)論：協(xié)方差法和修正協(xié)方差法可以分辨出本題中兩個頻率，說明它們相比于自相關(guān)法有更高的分辨率。這是因為它在形成的自相關(guān)估計是不需要對數(shù)據(jù)加窗。Burg算法也沒有對數(shù)據(jù)加窗，因此其AR參數(shù)的估計也比自相關(guān)估計所獲得的估計要精確。附代碼%x(n)=a(1)x(n-1)+a(2)x(n-2)+a(3)x(n-3)+a(4)x(n-4)+b(O)v(n)%v(n)是單位方差白噪聲%Pb(0)=1,a=2.7

10、607,-3.8106,2.6535,-0.9238clc;clear;closeall;N=256;M=1024;B=1；A=1-2.76073.8106-2.65350.9238;a=2.7607,-3.8106,2.6535,-0.9238;b0=1;nf=0:1/1000:1;w=pi*nf;v=randn(1,N);x=filter(B,A,v);figure;stem(x);title(AR過程生成序列x(n);%估計自相關(guān)Nr=N;r1=zeros(1,Nr);fork=0:(Nr-1)forn=1:(N-k)r1(k+1)=r1(k+1)+x(n)*x(n+k);endr1(k

11、+1)=r1(k+1)/N;endfigure;stem(r1,rx);holdon;%計算真是自相關(guān)r=ator(A,b0);fork=length(a)+2:Nr(k)=a(1)*r(k-1)+a(2)*r(k-2)+a(3)*r(k-3)+a(4)*r(k-4);endstem(r,b*);title(自相關(guān)估計);legend(估計自相關(guān),真實自相關(guān));gridon;%利用估計的自相關(guān)求功率譜P=zeros(1,M);P2=zeros(1,M);P_t=zeros(1,M);f=0:1/(M/2-1):1;Y=fft(x,M);P=(Y.*conj(Y)/N;Py=10*log10(P

12、);figure;plot(f,Py(1:length(f),-);holdon;%真實功率譜forh=1:length(nf)temp=0;fork=1:5temp=temp+A(k).*exp(-j*(k-1)*w(h);endP_a(h)=1/(temp.*conj(temp);endP_a=10*log10(abs(P_a);plot(nf,P_a,r);title(周期法獲取的功率譜和真實功率譜的比較);legend(周期圖功率譜，真實功率譜);xlabel(數(shù)字頻率八pi);ylabel(功率譜/dB);gridon;%利用所估計的自相關(guān)，估計AR參數(shù)a(1),a(2),a(3),

13、a(4),b0est_a,epsilon=rtoa(r1(1:5)est_bO=sqrt(epsilon)forh=1:length(nf)temp=0;fork=1:5temp=temp+est_a(k).*exp(-j*(k-1)*w(h);endP_acor(h)=epsilon/(temp.*conj(temp);endP_acor=10*log10(abs(P_acor);figure;plot(nf,P_acor,-);title(利用Yule-Walker法估計功率譜);xlabel(數(shù)字頻率/pi);ylabel(功率譜/dB);gridon;figure;plot(nf,P_

14、a,r);holdon;plot(nf,P_acor,-);title(功率譜比較);legend(真實功率譜,Yule_Walker法);xlabel(數(shù)字頻率(pi);ylabel(功率譜(dB)gridon;%利用協(xié)方差求功率譜p=4;a_cov,err_cov=covm(x,p)forh=1:length(nf)temp=0;fork=1:5temp=temp+a_cov(k)*exp(-j*(k-1)*w(h);endP_cov(h)=x(1)八2/(temp.*conj(temp);endP_cov=10*log10(abs(P_cov);figure(6);plot(nf,P_a

15、,r);gridon;holdon;plot(nf,P_cov,-);title(功率譜比較);legend(真實功率譜,協(xié)方差法);xlabel(數(shù)字頻率/pi);ylabel(PowerSpectrum/dB);%修正協(xié)方差法a_mcov,err_mcov=mcov(x,4)forh=1:length(nf)temp=0;fork=1:5temp=temp+a_mcov(k)*exp(-j*(k-1)*w(h);endP_mcov(h)=x(1)A2/(temp.*conj(temp);endP_mcov=10*log10(abs(P_mcov);figure;plot(nf,P_a,r);gridon;holdon;plot(nf,P_mcov,-);title(功率譜比較);legend(真實功率譜,修正協(xié)方差);xlabel(數(shù)字頻率/pi);ylabel(功率譜/dB);%Burg法求功率譜g_Burg=burg(x,4);a_Burg=gtoa(g_Burg)forh