AR模型功率譜估計實驗

1、精選優(yōu)質(zhì)文檔-傾情為你奉上實驗二一、實驗內(nèi)容AR過程的線性建模與功率譜估計。考慮AR過程：是單位方差白噪聲。(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個樣點。(b) 計算其自相關(guān)序列的估計，并與真實的自相關(guān)序列值相比較。(c) 將的DTFT作為x(n)的功率譜估計，即：。(d) 利用所估計的自相關(guān)值和Yule-Walker法(自相關(guān)法)，估計和的值，并討論估計的精度。(e) 用(d)中所估計的和來估計功率譜為：。(f) 將(c)和(e)的兩種功率譜估計與實際的功率譜進行比較，畫出它們

2、的重疊波形。(g) 重復(fù)上面的(d)(f)，只是估計AR參數(shù)分別采用如下方法：(1) 協(xié)方差法；(2) Burg方法；(3) 修正協(xié)方差法。試比較它們的功率譜估計精度。二、實驗分析1、 計算真實的自相關(guān)值時，采用逆Levinson-Durbin遞歸方法，由a、b參數(shù)得到，其中為濾波器的階數(shù)，再采用公式外推得到的自相關(guān)值；2、 實際功率譜，可調(diào)用Matlab中的FFT算法得到；3、 自相關(guān)序列的估計值采用公式得到；4、 采用各種功率譜估計方法對功率譜進行估計。三、實驗結(jié)果及分析仿真參數(shù)設(shè)置：采樣點數(shù)為64，頻域采樣點數(shù)為128自相關(guān)序列的估計與真實自相關(guān)序列值的比較見圖1，由圖可知估計值與真實值

3、存在一定的誤差，但整體變化趨勢相差不大。圖1 自相關(guān)序列題目（c）中功率譜的估計方法實際為周期圖法，周期圖法估計的功率譜與自相關(guān)法估計的功率譜的比較見圖2，由圖可知，周期圖能辨認出兩個峰值，而自相關(guān)法不能，說明周期圖的分辨率大于自相關(guān)法。圖2 周期圖法和自相關(guān)法得到的功率譜圖3圖7的（a）部分分別為采用周期圖法、自相關(guān)法、協(xié)方差法、Burg方法、修正協(xié)方差法進行功率譜50次估計的交疊圖，（b）部分給出了其整體平均及真實的功率譜。由這些圖可以看出，對于這一AR（4）過程，除自相關(guān)法外，所有估計都能分辨出兩個峰值，且峰值的位置大致相似。此外，周期圖法的方差大于其它估計方法。（a）(b)圖3 周期圖

4、法估計AR(4)過程的功率譜（a）(b)圖4 自相關(guān)法估計AR(4)過程的功率譜（a）(b)圖5 協(xié)方差法估計AR(4)過程的功率譜（a）(b)圖6 Burg法估計AR(4)過程的功率譜（a）(b)圖7 修正協(xié)方差法估計AR(4)過程的功率譜表1為采用自相關(guān)法、協(xié)方差法、Burg方法、修正協(xié)方差法得到的a參數(shù)和b參數(shù)，表中平方誤差和的計算公式為，從表中可以看出，自相關(guān)法估計的參數(shù)與真實值相比相差較大，協(xié)方差法、Burg方法、修正協(xié)方差法參數(shù)估計的性能相當。a(1)a(2)a(3)a(4)b(0)平方誤差和真實值2.7607-3.81062.6535-0.923810自相關(guān)法1.5371-1.3

5、2120.36348-0.134374.728913.5619協(xié)方差法2.7324-3.73732.5794-0.88852設(shè)為10.0129Burg方法2.7083-3.66992.5098-0.8569設(shè)為10.0477修正協(xié)方差法2.7071-3.66752.5069-0.85566設(shè)為10.0495表1 各種方法估計的a參數(shù)和b參數(shù)四、實驗思考當觀測數(shù)據(jù)存在觀測噪聲時，即，其中是單位方差的白噪聲，與不相關(guān)，考慮觀測噪聲對各種譜估計方法的影響。仿真參數(shù)設(shè)置：采樣點數(shù)為64，頻域采樣點數(shù)為128。圖8圖11為存在觀測噪聲時，各種方法對AR（4）過程的功率譜估計，由圖可知，觀測噪聲對各種譜估

6、計方法的估計性能具有一定的影響，在上述仿真參數(shù)設(shè)置下，周期圖法和協(xié)方差法可辨認出兩個峰值，而Burg方法和修正協(xié)方差法不能。不過，當增加頻域采樣點數(shù)時，Burg方法和修正協(xié)方差法也能辨認出兩個峰值，但頻域采樣點數(shù)的增加意味著計算量的增加。圖8 觀測噪聲下周期圖的功率譜估計圖9 觀測噪聲下協(xié)方差法的功率譜估計圖10 觀測噪聲下Burg方法的功率譜估計圖11 觀測噪聲下修正協(xié)方差法的功率譜估計五、實驗源代碼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=6

7、4;%樣點數(shù)L=128;%頻率采樣點數(shù)times=50;%運行次數(shù)index=0:N-1;r_real=zeros(1,N);%計算自相關(guān)序列的真實值，只有前5個值r_real(1:(p+1)=ator(A,b0);%進行自相關(guān)的外推for k=(p+2):N %教科書p99 公式3.6 r_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.matx=;

8、for i=1:times xt=filter(b0,A,v(i,:); x=x;xt;end%周期圖法Tr_e=;TPx_e=;figure(1)for i=1:times r_e=E_r(x(i,:),N);%教科書p78 公式2.204 Tr_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)/ti

9、mes;TPx_e=sum(TPx_e)/times;figure(2)plot(w,abs(Px_real),r-,w,abs(TPx_e),b-)xlabel(頻率/pi)ylabel(功率譜的幅值)title(周期圖法50次平均與真實功率譜比較)legend(真實值,周期圖法)grid on;%自相關(guān)法figure(3)TA_YW=;Tb0_YW=;TPx_YW=;for i=1:times r_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.;

10、 Tb0_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

11、)plot(w,abs(Px_real),r-,w,abs(TPx_YW),b-)xlabel(頻率/pi)ylabel(功率譜的幅值)title(自相關(guān)法50次平均與真實功率譜比較)legend(真實值,自相關(guān)法)grid on;%協(xié)方差法figure(5)TA_cov=;TPx_cov=;for i=1:times A_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

12、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次平均與真實功率譜比較)legend(真實值,協(xié)方差法)grid on;%Burg方法 figure(7)TA_Bur

13、g=;TPx_Burg=;for i=1:times gamma_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次交疊圖)gr

14、id 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次平均與真實功率譜比較)legend(真實值,Burg方法)grid on;%修正協(xié)方差法figure(9)TA_mcov=;TPx_mcov=;for i=1:times A_mcov,err_mcov=mcov(x(i,:),p); TA_mcov=TA_mcov;A_m

15、cov.; b0_mcov=1; hx_mcov=dimpulse(b0_mcov,A_mcov.,L).; Px_mcov=abs(fft(hx_mcov,L).2; TPx_mcov=TPx_mcov;Px_mcov; plot(w,abs(Px_mcov); hold onendxlabel(頻率/pi)ylabel(功率譜的幅值)title(修正協(xié)方差法50次交疊圖)grid on;hold offTA_mcov=sum(TA_mcov)/times;TPx_mcov=sum(TPx_mcov)/times;figure(10)plot(w,abs(Px_real),r-,w,abs(TPx_mcov),b-)xlabel(頻率/pi)ylabel(功率譜的幅值)title(修正協(xié)方差法50次平均與真實功率譜比較)legend(真實值,修正