AR模型功率譜估計(jì)實(shí)驗(yàn)

1、實(shí)驗(yàn)二一、實(shí)驗(yàn)內(nèi)容AR 過(guò)程的線性建模與功率譜估計(jì)?？紤]AR 過(guò)程：x(n)a(1)x(n1)a(2) x(n2)a(3) x(n3)a(4) x(n4)b(0)v( n)v(n) 是單位方差白噪聲。(a) 取 b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238 ，產(chǎn)生 x(n)的 N=64 個(gè)樣點(diǎn)。(b) 計(jì)算其自相關(guān)序列的估計(jì)r?x( k ) ，并與真實(shí)的自相關(guān)序列值相比較。(c)將 ?的DTFT作為x(n)的功率譜估計(jì)，即：r (k)x?jN 1jk1j2?。Px( e )(k)e| X (e ) |rxNkN 1利用

2、所估計(jì)的自相關(guān)值和 Yule-Walker 法( 自相關(guān)法 )，估計(jì) a(1), a(2), a(3), a(4) 和 b(0) 的值，并討論估計(jì)的精度。(e) 用 (d)中所估計(jì)的 a(k ) 和 b(0)?( ej| b(0) |2來(lái)估計(jì)功率譜為： Px)2 。41a(k )e jkk 1將 (c)和 (e)的兩種功率譜估計(jì)與實(shí)際的功率譜進(jìn)行比較，畫出它們的重疊波形。(g)重復(fù)上面的 (d)(f) ，只是估計(jì)AR 參數(shù)分別采用如下方法：(1)協(xié)方差法； (2) Burg 方法； (3)修正協(xié)方差法。試比較它們的功率譜估計(jì)精度。二、實(shí)驗(yàn)分析1、 計(jì)算真實(shí)的自相關(guān)值時(shí)，采用逆Levinson-

3、Durbin遞歸方法，由a、 b 參數(shù)得到 rx (0) ， rx (1)，L ， rx ( p) ，其中 p 為濾波器的階數(shù)，再采用公式prxkap lrx kl 外推得到 k fp 的自l 1相關(guān)值；N 1N 1N12、 實(shí)際功率譜 Px (ej )rx (k )e jkrx (k )e jkrx (k )e jkrx 0 ，k N 1k 0k0N 1rx (k )e jk 可調(diào)用 Matlab 中的 FFT 算法得到；k 03rxkx n x nk得到；、 自相關(guān)序列的估計(jì)值采用公式?1 NN n14、 采用各種功率譜估計(jì)方法對(duì)功率譜進(jìn)行估計(jì)。三、實(shí)驗(yàn)結(jié)果及分析仿真參數(shù)設(shè)置：采樣點(diǎn)數(shù)為6

4、4，頻域采樣點(diǎn)數(shù)為128自相關(guān)序列的估計(jì)r?x( k) 與真實(shí)自相關(guān)序列值的比較見圖1，由圖可知估計(jì)值與真實(shí)值存在一定的誤差，但整體變化趨勢(shì)相差不大。800下 下下600下 下下400200下下下0下下下 -200 -400 -600 -800010203040506070下下下下圖 1自相關(guān)序列題目（ c）中功率譜的估計(jì)方法實(shí)際為周期圖法，周期圖法估計(jì)的功率譜與自相關(guān)法估計(jì)的功率譜的比較見圖 2，由圖可知，周期圖能辨認(rèn)出兩個(gè)峰值，而自相關(guān)法不能，說(shuō)明周期圖的分辨率大于自相關(guān)法。18000下 下 下16000下下下下下下下下1400012000下下 10000下下8000下下600040002

5、00000.20.40.60.811.21.41.61.820下 下 /圖 2周期圖法和自相關(guān)法得到的功率譜圖 3圖 7 的（ a）部分分別為采用周期圖法、自相關(guān)法、協(xié)方差法、Burg 方法、修正協(xié)方差法進(jìn)行功率譜 50 次估計(jì)的交疊圖， （b）部分給出了其整體平均及真實(shí)的功率譜。由這些圖可以看出，對(duì)于這一 AR （ 4）過(guò)程，除自相關(guān)法外，所有估計(jì)都能分辨出兩個(gè)峰值，且峰值的位置大致相似。此外，周期圖法的方差大于其它估計(jì)方法。x 104下下下下50下下下下3.532.5下2下下下下 1.5下10.5000.20.40.60.811.21.41.61.82下 下 /a）圖 3周期圖法估計(jì)x 1

6、04下下下下50下下下下4.543.53下下 2.5下下2下下1.510.5000.20.40.60.811.21.41.61.82下 下 /a）圖 4自相關(guān)法估計(jì)x 106下下下下50下下下下32.52下下下 1.5下下下10.5000.20.40.60.811.21.41.61.82下 下 /a）圖 5協(xié)方差法估計(jì)18000下下下下50下下下下下下下下下下下下 下 下16000下下下下1400012000下下 10000下下 8000下下600040002000000.20.40.60.811.21.41.61.82下 下 /(b)AR(4) 過(guò)程的功率譜下下下下50下下下下下下下下下下下

7、18000下 下 下16000下下下下1400012000下下 10000下下 8000下下600040002000000.20.40.60.811.21.41.61.82下 下 /(b)AR(4) 過(guò)程的功率譜x 104下下下下50下下下下下下下下下下下7下 下 下下下下下65下下4下下下 3下21000.20.40.60.811.21.41.61.82下 下 /(b)AR(4) 過(guò)程的功率譜x 104Burg 下 下 50 下 下下 下Burg 下 下 50 下 下 下 下 下 下 下 下 下 下 下180006下 下 下16000Burg 下下514000412000下下下下 10000

8、下 3下下下8000下下下下2600014000200000.20.40.60.811.21.41.61.82000.20.40.60.811.21.41.61.82下 下 /0下 下 /（a）(b)圖 6 Burg 法估計(jì) AR(4) 過(guò)程的功率譜x 104下下下下下下50下下下下下下下下下下50下下下下下下下下下下下618000下下 下16000下下下下下下514000412000下下下下 10000下 3下下下8000下下下下26000140002000000.20.40.60.811.21.41.61.8200.20.40.60.811.21.41.61.82下 下 /0下 下 /（a

9、）(b)圖 7修正協(xié)方差法估計(jì)AR(4) 過(guò)程的功率譜表 1 為采用自相關(guān)法、協(xié)方差法、Burg 方法、修正協(xié)方差法得到的a 參數(shù)和 b 參數(shù)，表中平方4?2誤差和的計(jì)算公式為，從表中可以看出，自相關(guān)法估計(jì)的參數(shù)與真實(shí)值相比相a iaii1差較大，協(xié)方差法、Burg 方法、修正協(xié)方差法參數(shù)估計(jì)的性能相當(dāng)。a(1)a(2)a(3)a(4)b(0)平方誤差和真實(shí)值2.7607-3.81062.6535-0.923810自相關(guān)法1.5371-1.32120.36348-0.134374.728913.5619協(xié)方差法2.7324-3.73732.5794-0.88852設(shè)為 10.0129Burg

10、方法2.7083-3.66992.5098-0.8569設(shè)為 10.0477修正協(xié)方差法2.7071-3.66752.5069-0.85566設(shè)為 10.0495表 1 各種方法估計(jì)的a 參數(shù)和 b 參數(shù)四、實(shí)驗(yàn)思考當(dāng)觀測(cè)數(shù)據(jù)存在觀測(cè)噪聲時(shí)，即y nx nw n ，其中 w n 是單位方差的白噪聲，與v n 不相關(guān)，考慮觀測(cè)噪聲對(duì)各種譜估計(jì)方法的影響。仿真參數(shù)設(shè)置：采樣點(diǎn)數(shù)為64，頻域采樣點(diǎn)數(shù)為128。圖 8圖 11 為存在觀測(cè)噪聲時(shí)，各種方法對(duì) AR （ 4）過(guò)程的功率譜估計(jì)，由圖可知，觀測(cè)噪聲對(duì)各種譜估計(jì)方法的估計(jì)性能具有一定的影響，在上述仿真參數(shù)設(shè)置下，周期圖法和協(xié)方差法可辨認(rèn)出兩個(gè)峰值

11、，而Burg 方法和修正協(xié)方差法不能。不過(guò)，當(dāng)增加頻域采樣點(diǎn)數(shù)時(shí)，Burg 方法和修正協(xié)方差法也能辨認(rèn)出兩個(gè)峰值，但頻域采樣點(diǎn)數(shù)的增加意味著計(jì)算量的增加。下下下下50下下下下下下下50下下下60001800016000500014000400012000下下下下 10000下 3000下下下8000下下下下2000600040001000200000.20.40.60.811.21.41.61.8200.20.40.60.811.21.41.61.8200下 下 /下 下 /圖 8觀測(cè)噪聲下周期圖的功率譜估計(jì)圖 9觀測(cè)噪聲下協(xié)方差法的功率譜估計(jì)800Burg 下 下 50下 下 下900下下下

12、下下下50下下下700800600700500600下下下下 500下 400下下下400下下下 300下30020020010010000.20.40.60.811.21.41.61.8200.20.40.60.811.21.41.61.8200下 下 /下 下 /圖 10觀測(cè)噪聲下 Burg 方法的功率譜估計(jì)圖 11觀測(cè)噪聲下修正協(xié)方差法的功率譜估計(jì)五、實(shí)驗(yàn)源代碼clc;clear;clf;%a 參數(shù)a1=2.7607;a2=-3.8106;a3=2.6535;a4=-0.9238;p=4;%A=1 -a1 -a2 -a3 -a4;b0=1;%b 參數(shù)N=64;% 樣點(diǎn)數(shù)L=128;% 頻

13、率采樣點(diǎn)數(shù)times=50;% 運(yùn)行次數(shù)index=0:N-1;r_real=zeros(1,N);%計(jì)算自相關(guān)序列的真實(shí)值，只有前5 個(gè)值r_real(1:(p+1)=ator(A,b0);%進(jìn)行自相關(guān)的外推for k=(p+2):N%教科書 p99 公式 3.6r_real(k)=-sum(A(2:p+1).* seqreverse(r_real(k-p:k-1);endPx_real=fft(r_real,L)+conj(fft(r_real,L)-r_real(1);w=(1:L)-1)/L*2;v=randn(times,N);load v1.mat x=;for i=1:times

14、 xt=filter(b0,A,v(i,:); x=x;xt;end%周期圖法Tr_e=;TPx_e=;figure(1)for i=1:timesr_e=E_r(x(i,:),N);% 教科書 p78 公式 2.204Tr_e=Tr_e;r_e;Px_e=fft(r_e,L)+conj(fft(r_e,L)-r_e(1);TPx_e=TPx_e;Px_e;plot(w,abs(Px_e);hold onendxlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title( 周期圖法 50 次交疊圖 )grid on;hold offTr_e=sum(Tr_e)/times;TPx_

15、e=sum(TPx_e)/times;figure(2)plot(w,abs(Px_real),r-,w,abs(TPx_e),b-)xlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title( 周期圖法 50 次平均與真實(shí)功率譜比較)legend(真實(shí)值 ,周期圖法 )grid on;%自相關(guān)法figure(3)TA_YW=;Tb0_YW=;TPx_YW=;for i=1:timesr_e=E_r(x(i,:),N);r_ep=r_e(1:(p+1);A_YW,err_YW=rtoa(r_ep);b0_YW=sqrt(err_YW);TA_YW=TA_YW;A_YW.;Tb0_

16、YW=Tb0_YW;b0_YW;%先求出系統(tǒng)的沖擊響應(yīng)，然后得到頻域幅度，從而求出功率譜hx_YW=dimpulse(b0_YW,A_YW.,L).;Px_YW=abs(fft(hx_YW,L).2;TPx_YW=TPx_YW;Px_YW;plot(w,abs(Px_YW);hold onendxlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title( 自相關(guān)法 50 次交疊圖 )grid on;hold offTA_YW=sum(TA_YW)/times;Tb0_YW=sum(Tb0_YW)/times;TPx_YW=sum(TPx_YW)/times;figure(4)pl

17、ot(w,abs(Px_real),r-,w,abs(TPx_YW),b-)xlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title( 自相關(guān)法 50 次平均與真實(shí)功率譜比較)legend(真實(shí)值 ,自相關(guān)法 )grid on;%協(xié)方差法figure(5)TA_cov=;TPx_cov=;for i=1:timesA_cov,err_cov=covm(x(i,:),p);TA_cov=TA_cov;A_cov.;b0_cov=1;hx_cov=dimpulse(b0_cov,A_cov.,L).;Px_cov=abs(fft(hx_cov,L).2;TPx_cov=TPx_co

18、v;Px_cov;plot(w,abs(Px_cov);hold onendxlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title( 協(xié)方差法 50 次交疊圖 )grid on;hold offTA_cov=sum(TA_cov)/times;TPx_cov=sum(TPx_cov)/times;figure(6)plot(w,abs(Px_real),r-,w,abs(TPx_cov),b-)xlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title( 協(xié)方差法 50 次平均與真實(shí)功率譜比較)legend(真實(shí)值 ,協(xié)方差法 )grid on;%Burg 方法

19、figure(7)TA_Burg=;TPx_Burg=;for i=1:timesgamma_Burg,err_Burg = burg(x(i,:),p);A_Burg=gtoa(gamma_Burg);TA_Burg=TA_Burg;A_Burg.;b0_Burg=1;hx_Burg=dimpulse(b0_Burg,A_Burg.,L).;Px_Burg=abs(fft(hx_Burg,L).2;TPx_Burg=TPx_Burg;Px_Burg;plot(w,abs(Px_Burg);hold onendxlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title(Burg 方法 50 次交疊圖 )grid on;hold offTA_Burg=sum(TA_Burg)/times;TPx_Burg=sum(TPx_Burg)/times;figure(8)plot(w,abs(Px_real),r-,w,abs(TPx_Burg),b-)xlabel( 頻率 /pi)ylabel( 功率譜的幅值 )title(Burg 方法 50 次平均與真實(shí)功率譜比較)legend(真實(shí)值 ,Burg 方法 )grid on;%修正協(xié)方差法figure(9)TA_mcov=;TPx_mcov=;for i=1:timesA_mcov,err_mcov