AR模型功率谱估计实验

实验二一、实验内容AR过程的线性建模与功率谱估计。考虑AR过程：()(1)(1)(2)(2)(3)(3)(4)(4)(0)()xnaxnaxnaxnaxnbvn()vn是单位方差白噪声。(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)计算其自相关序列的估计ˆ()xrk，并与真实的自相关序列值相比较。(c)将ˆ()xrk的DTFT作为x(n)的功率谱估计，即：1211ˆˆ()()|()|NjjkjxxkNPerkeXeN。(d)利用所估计的自相关值和Yule-Walker法(自相关法)，估计(1),(2),(3),(4)aaaa和(0)b的值，并讨论估计的精度。(e)用(d)中所估计的()ak和(0)b来估计功率谱为：2241|(0)|ˆ()1()jxjkkbPeake。(f)将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。(g)重复上面的(d)~(f)，只是估计AR参数分别采用如下方法：(1)协方差法；(2)Burg方法；(3)修正协方差法。试比较它们的功率谱估计精度。二、实验分析1、计算真实的自相关值时，采用逆Levinson-Durbin递归方法，由a、b参数得到(0)xr，(1)xr，，()xrp，其中p为滤波器的阶数，再采用公式1pxpxlrkalrkl外推得到kp的自相关值；2、实际功率谱111100()()()()0NNNjjkjkjkxxxxxkNkkPerkerkerker，10()Njkxkrke可调用Matlab中的FFT算法得到；3、自相关序列的估计值采用公式11ˆNxnrkxnxnkN得到；4、采用各种功率谱估计方法对功率谱进行估计。三、实验结果及分析仿真参数设置：采样点数为64，频域采样点数为128自相关序列的估计ˆ()xrk与真实自相关序列值的比较见图1，由图可知估计值与真实值存在一定的误差，但整体变化趋势相差不大。010203040506070-800-600-400-2000200400600800下标序列自相关序列值真实值估计值图1自相关序列题目（c）中功率谱的估计方法实际为周期图法，周期图法估计的功率谱与自相关法估计的功率谱的比较见图2，由图可知，周期图能辨认出两个峰值，而自相关法不能，说明周期图的分辨率大于自相关法。00.20.40.60.811.21.41.61.82020004000600080001000012000140001600018000频率/功率谱的幅值真实值周期图法自相关法图2周期图法和自相关法得到的功率谱图3～图7的（a）部分分别为采用周期图法、自相关法、协方差法、Burg方法、修正协方差法进行功率谱50次估计的交叠图，（b）部分给出了其整体平均及真实的功率谱。由这些图可以看出，对于这一AR（4）过程，除自相关法外，所有估计都能分辨出两个峰值，且峰值的位置大致相似。此外，周期图法的方差大于其它估计方法。00.20.40.60.811.21.41.61.8200.511.522.533.5x104频率/功率谱的幅值周期图法50次交叠图（a）00.20.40.60.811.21.41.61.82020004000600080001000012000140001600018000频率/功率谱的幅值周期图法50次平均与真实功率谱比较真实值周期图法(b)图3周期图法估计AR(4)过程的功率谱00.20.40.60.811.21.41.61.8200.511.522.533.544.5x104频率/功率谱的幅值自相关法50次交叠图（a）00.20.40.60.811.21.41.61.82020004000600080001000012000140001600018000频率/功率谱的幅值自相关法50次平均与真实功率谱比较真实值自相关法(b)图4自相关法估计AR(4)过程的功率谱00.20.40.60.811.21.41.61.8200.511.522.53x106频率/功率谱的幅值协方差法50次交叠图（a）00.20.40.60.811.21.41.61.8201234567x104频率/功率谱的幅值协方差法50次平均与真实功率谱比较真实值协方差法(b)图5协方差法估计AR(4)过程的功率谱00.20.40.60.811.21.41.61.820123456x104频率/功率谱的幅值Burg方法50次交叠图（a）00.20.40.60.811.21.41.61.82020004000600080001000012000140001600018000频率/功率谱的幅值Burg方法50次平均与真实功率谱比较真实值Burg方法(b)图6Burg法估计AR(4)过程的功率谱00.20.40.60.811.21.41.61.820123456x104频率/功率谱的幅值修正协方差法50次交叠图（a）00.20.40.60.811.21.41.61.82020004000600080001000012000140001600018000频率/功率谱的幅值修正协方差法50次平均与真实功率谱比较真实值修正协方差法(b)图7修正协方差法估计AR(4)过程的功率谱表1为采用自相关法、协方差法、Burg方法、修正协方差法得到的a参数和b参数，表中平方误差和的计算公式为241ˆiaiai，从表中可以看出，自相关法估计的参数与真实值相比相差较大，协方差法、Burg方法、修正协方差法参数估计的性能相当。a(1)a(2)a(3)a(4)b(0)平方误差和真实值2.7607-3.81062.6535-0.923810自相关法1.5371-1.32120.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参数四、实验思考当观测数据存在观测噪声时，即ynxnwn，其中wn是单位方差的白噪声，与vn不相关，考虑观测噪声对各种谱估计方法的影响。仿真参数设置：采样点数为64，频域采样点数为128。图8～图11为存在观测噪声时，各种方法对AR（4）过程的功率谱估计，由图可知，观测噪声对各种谱估计方法的估计性能具有一定的影响，在上述仿真参数设置下，周期图法和协方差法可辨认出两个峰值，而Burg方法和修正协方差法不能。不过，当增加频域采样点数时，Burg方法和修正协方差法也能辨认出两个峰值，但频域采样点数的增加意味着计算量的增加。00.20.40.60.811.21.41.61.820100020003000400050006000频率/功率谱的幅值周期图法50次平均图8观测噪声下周期图的功率谱估计00.20.40.60.811.21.41.61.82020004000600080001000012000140001600018000频率/功率谱的幅值协方差法50次平均图9观测噪声下协方差法的功率谱估计00.20.40.60.811.21.41.61.820100200300400500600700800频率/功率谱的幅值Burg方法50次平均图10观测噪声下Burg方法的功率谱估计00.20.40.60.811.21.41.61.820100200300400500600700800900频率/功率谱的幅值修正协方差法50次平均图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=64;%样点数L=128;%频率采样点数times=50;%运行次数index=0:N-1;r_real=zeros(1,N);%计算自相关序列的真实值，只有前5个值r_real(1:(p+1))=ator(A,b0);%进行自相关的外推fork=(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);%loadv1.matx=[];fori=1:timesxt=filter(b0,A,v(i,:));x=[x;xt];end%周期图法Tr_e=[];TPx_e=[];figure(1)fori=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));holdonendxlabel('频率/\pi')ylabel('功率谱的幅值')title('周期图法50次交叠图')gridon;holdoffTr_e=sum(Tr_e)/times;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('真实值','周期图法')gridon;%自相关法figure(3)TA_YW=[];Tb0_YW=[];TPx_YW=[];fori=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_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));holdonendxlabel('频率/\pi')ylabel('功率谱的幅值')title('自相关法50次交叠图')gridon;holdoffTA_YW=sum(TA_YW)/times;Tb0_YW=sum(Tb0_YW)/times;TPx_YW=sum(TPx_YW)/times;figure(4)plot(w,abs(Px_real),'r-',w,abs(TPx_YW),'b--')xlabel('频率/\pi')ylabel('功率谱的幅值')title('自相关法50次平均与真实功率谱比较')legend('真实值','自相关法')gridon;%协方差法figure(5)TA_cov=[];TPx_cov=[];fori=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_cov;Px_cov];plot(w,abs(Px_