《AR模型功率谱估计实验(共9页).doc》由会员分享，可在线阅读，更多相关《AR模型功率谱估计实验(共9页).doc（9页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上实验二一、实验内容AR过程的线性建模与功率谱估计。考虑AR过程：是单位方差白噪声。(a) 取b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生x(n)的N=64个样点。(b) 计算其自相关序列的估计，并与真实的自相关序列值相比较。(c) 将的DTFT作为x(n)的功率谱估计，即：。(d) 利用所估计的自相关值和Yule-Walker法(自相关法)，估计和的值，并讨论估计的精度。(e) 用(d)中所估计的和来估计功率谱为：。(f) 将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们
2、的重叠波形。(g) 重复上面的(d)(f)，只是估计AR参数分别采用如下方法：(1) 协方差法；(2) Burg方法；(3) 修正协方差法。试比较它们的功率谱估计精度。二、实验分析1、 计算真实的自相关值时，采用逆Levinson-Durbin递归方法，由a、b参数得到，其中为滤波器的阶数，再采用公式外推得到的自相关值；2、 实际功率谱，可调用Matlab中的FFT算法得到；3、 自相关序列的估计值采用公式得到；4、 采用各种功率谱估计方法对功率谱进行估计。三、实验结果及分析仿真参数设置：采样点数为64，频域采样点数为128自相关序列的估计与真实自相关序列值的比较见图1，由图可知估计值与真实值
3、存在一定的误差，但整体变化趋势相差不大。图1 自相关序列题目（c）中功率谱的估计方法实际为周期图法，周期图法估计的功率谱与自相关法估计的功率谱的比较见图2，由图可知，周期图能辨认出两个峰值，而自相关法不能，说明周期图的分辨率大于自相关法。图2 周期图法和自相关法得到的功率谱图3图7的（a）部分分别为采用周期图法、自相关法、协方差法、Burg方法、修正协方差法进行功率谱50次估计的交叠图，（b）部分给出了其整体平均及真实的功率谱。由这些图可以看出，对于这一AR（4）过程，除自相关法外，所有估计都能分辨出两个峰值，且峰值的位置大致相似。此外，周期图法的方差大于其它估计方法。（a）(b)图3 周期图
4、法估计AR(4)过程的功率谱（a）(b)图4 自相关法估计AR(4)过程的功率谱（a）(b)图5 协方差法估计AR(4)过程的功率谱（a）(b)图6 Burg法估计AR(4)过程的功率谱（a）(b)图7 修正协方差法估计AR(4)过程的功率谱表1为采用自相关法、协方差法、Burg方法、修正协方差法得到的a参数和b参数，表中平方误差和的计算公式为，从表中可以看出，自相关法估计的参数与真实值相比相差较大，协方差法、Burg方法、修正协方差法参数估计的性能相当。a(1)a(2)a(3)a(4)b(0)平方误差和真实值2.7607-3.81062.6535-0.923810自相关法1.5371-1.3
5、2120.36348-0.134374.728913.5619协方差法2.7324-3.73732.5794-0.88852设为10.0129Burg方法2.7083-3.66992.5098-0.8569设为10.0477修正协方差法2.7071-3.66752.5069-0.85566设为10.0495表1 各种方法估计的a参数和b参数四、实验思考当观测数据存在观测噪声时，即，其中是单位方差的白噪声，与不相关，考虑观测噪声对各种谱估计方法的影响。仿真参数设置：采样点数为64，频域采样点数为128。图8图11为存在观测噪声时，各种方法对AR（4）过程的功率谱估计，由图可知，观测噪声对各种谱估
6、计方法的估计性能具有一定的影响，在上述仿真参数设置下，周期图法和协方差法可辨认出两个峰值，而Burg方法和修正协方差法不能。不过，当增加频域采样点数时，Burg方法和修正协方差法也能辨认出两个峰值，但频域采样点数的增加意味着计算量的增加。图8 观测噪声下周期图的功率谱估计图9 观测噪声下协方差法的功率谱估计图10 观测噪声下Burg方法的功率谱估计图11 观测噪声下修正协方差法的功率谱估计五、实验源代码clc;clear;clf;%a参数a1=2.7607;a2=-3.8106;a3=2.6535;a4=-0.9238;p=4;%A=1 -a1 -a2 -a3 -a4;b0=1;%b参数N=6
7、4;%样点数L=128;%频率采样点数times=50;%运行次数index=0:N-1;r_real=zeros(1,N);%计算自相关序列的真实值，只有前5个值r_real(1:(p+1)=ator(A,b0);%进行自相关的外推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;%自相关法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; %先求出系统的冲击响应，然后得到频域幅度，从而求出功率谱 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(自相关法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(自相关法50次平均与真实功率谱比较)legend(真实值,自相关法)grid on;%协方差法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(协方差法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(协方差法50次平均与真实功率谱比较)legend(真实值,协方差法)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;%修正协方差法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(修正协方差法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(
16、TPx_mcov),b-)xlabel(频率/pi)ylabel(功率谱的幅值)title(修正协方差法50次平均与真实功率谱比较)legend(真实值,修正协方差法)grid on;figure(11);plot(index,r_real,r-,index,Tr_e,b-)xlabel(下标序列)ylabel(自相关序列值)legend(真实值,估计值)grid on; figure(12)plot(w,abs(Px_real),r-,w,abs(TPx_e),b-,w,abs(TPx_YW),g:)xlabel(频率/pi)ylabel(功率谱的幅值)legend(真实值,周期图法,自相关法)grid on; 专心-专注-专业