1、AR 过程的线性建模与功率谱估计一、实验内容：AR 过程的线性建模与功率谱估计。考虑 AR 过程：()1()(2)(3)(4)(0)xnaaxnaxnaxnbvn是单位方差白噪声。v(a) 取 b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生 x(n)的 N=256 个样点。(b) 计算其自相关序列的估计 ，并与真实的自相关序列值相比较。()xrk(c) 将 的 DTFT 作为 x(n)的功率谱估计，即：()xrk。1 21()|()|Nj jkjxkPereXe(d) 利用所估计的自相关值和 Yule-Walker
2、 法(自相关法)，估计 和(1), 2(3), 4aa的值，并讨论估计的精度。(0)b(e) 用(d)中所估计的 和 来估计功率谱为： 。()ak0b 241|0|()()jx jkkbPeae(f) 将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。(g) 重复上面的(d)(f)，只是估计 AR 参数分别采用如下方法：(1) 协方差法；(2) Burg方法；(3) 修正协方差法。试比较它们的功率谱估计精度。二、实验结果及分析问题一：取 b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生 x(
3、n)的 N=256 个样点。计算其自相关序列的估计 ，并与真实的自相关序列值相比较。()xrk思路：计算真实的自相关值时，采用逆 Levinson-Durbin 递归方法，由 a、b 参数得到 ，(0)xr， ， ，其中 为滤波器的阶数，再采用公式 外推得(1)xr ()xrp1pxxlrkrkl到 的自相关值；k结果分析：仿真参数设置：采样点数为 256自相关序列的估计 与真实自相关序列值的比较见图 1，由图可知估计值与真实值()xrk存在一定的误差，但整体变化趋势相差不大。0 50 100 150 200 250 300-800-600-400-2000200400600800下下下下下下
4、下下下下下下下下下下图 1 自相关估计和实际比较问题二：将 的 DTFT 作为 x(n)的功率谱估计，即：()xrk。1 21()|()|Nj jkjxkPereXe利用所估计的自相关值和 Yule-Walker 法(自相关法) ，用所估计的 和 来估计功()ak0b率谱为： 。241|(0)|jx jkkbeae将周期图法和自相关法的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。思路：实际功率谱 ，11100()()()()0NNNj jkjkjkxxxxxkPerererer可调用 Matlab 中的 FFT 算法得到；自相关序列的估计值采用公式10()Njkxkr得到1Nxn
5、xk结果分析：通过图 4 比较可知，周期图能辨认出两个峰值，而自相关法不能，说明周期图的分辨率大于自相关法。同时比较图 2 和图 3 可知，周期图的方差大于自相关法。0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2012345678x 104下下下下下下下下下下下下50下下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 104下下/下下下下下下下下下下50下下下下下下下下下下下下下下下下下下图 2 周期图法估计 AR(4)过程的功率谱0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 201
6、23456x 104下下下下下下下下下下下下50下下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 104下下/下下下下下下下下下下50下下下下下下下下下下下下下下下下下下图 3 自相关法估计 AR(4)过程的功率谱0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2-0.500.511.522.5x 104下下/下下下下下下下下下下下下下下下下下图 4 周期图法和自相关法得到的功率谱问题三：重复上面的(d)(f)，只是估计 AR 参数分别采用如下方法：(1) 协方差法；(2) Burg方法；(3) 修正协方差法。
7、试比较它们的功率谱估计精度。思路： 方法和问题二的思路相似，只是使用的方式不同而已结果分析：图 5图 7 的左图部分分别为采用协方差法、Burg 方法、修正协方差法进行功率谱50 次估计的交叠图，右图部分给出了其整体平均及真实的功率谱。由这些图可以看出，对于这一 AR（4）过程，这三种方法都可以分辨出两个峰值来，而且方差大小差不多。但是和图 2 和图 3 进行综合比较得，除自相关法外，所有估计都能分辨出两个峰值，且峰值的位置大致相似。同时，周期图法的方差大于其它估计方法。0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.533.54x 105下下下
8、下下下下下下下下下50下下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.533.5x 104下下/下下下下下下下下下下50下下下下下下下下下下下下下下下下下下图 5 协方差法估计 AR(4)过程的功率谱0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 105下下下下下下下下下下下下下下50下下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 104下下/下下下下下下Burg下下50下下下下下下下下下下下下下下Burg下图 6 Burg
9、 法估计 AR(4)过程的功率谱0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 105下下下下下下下下下下下下下下50下下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 104下下/下下下下下下下下下下下下50下下下下下下下下下下下下下下下下下下下下图 7 修正协方差法估计 AR(4)过程的功率谱问题四：分别采用自相关法、协方差法、Burg 方法、修正协方差法，估计和 的值，并讨论估计的精度。(1), 2(3), 4aa(0)b思路：采用平方误差和的计算公式为 ，来衡量四中方法的估计
10、精确度。241iai结果分析：从表 1 可知，自相关法估计的参数与真实值相比相差较大， Burg 方法和修正协方差法参数估计的性能相当，其中协方差法最精确。表 1 各种方法估计的 a 参数和 b 参数a(1) a(2) a(3) a(4) b(0) 平方误差真实值 -2.7606 3.8106 -2.6535 0.9238 1自相关法 -2.0703 2.2862 -1.1889 0.3676 3.1077 9.6969协方差法 -2.7527 3.7901 -2.6318 0.9138 1 0.0010Burg 法 -2.7465 3.7737 -2.6151 0.9063 1 0.0033
11、修正协方差法 -2.7472 3.7747 -2.6157 0.9063 1 0.0031三、实验深入当观测数据存在观测噪声时，即 ，其中 是单位方差的白噪ynxwnn声，与 不相关，考虑观测噪声对各种谱估计方法的影响。vn图 8 是在存在观测噪声情况下，周期图法、协方差法、Burg 方法和修正协方差法对AR（4 ）过程的功率谱估计，由图可知，观测噪声对各种谱估计方法的估计性能具有一定的影响，周期图法辨别能力相对比较强，其次是协方差法还可以辨别出两个峰值。但是 Burg 方法和修正协方差法不能。周期图法它的性能主要依赖于数据序列长度 N，N 越大分辨率越好。当增加频域采样点数时，Burg 方法
12、和修正协方差法也能辨认出两个峰值，但频域采样点数的增加意味着计算量的增加。仿真参数设置：噪声误差为 1，采样点数为 64，频域采样点数为 128。0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20100020003000400050006000下下/下下下下下下下下下下50下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 200.511.522.5x 104下下/下下下下下下下下下下50下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20100200300400500600700800下下/下下下下下下Burg
13、下下50下下下0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20100200300400500600700800900下下/下下下下下下下下下下下下50下下下图 8 加噪声情况下的功率估计方法比较附录：MATLAB 代码clear all;a1=2.7607;a2=-3.8106;a3=2.6535;a4=-0.9238;b0=1;t=50;p=4;%a=1,-a1,-a2,-a3,-a4;N=256;L=512;%频率采样点rz=zeros(1,N);rz(1:(p+1)=ator(a,b0);for k=(p+2):N;rz(k)=-sum(a(2:p+1).*
14、seqreverse(rz(k-p:k-1);end;pz=fft(rz,L)+conj(fft(rz,L)-rz(1);w=(1:L)-1)/L*2;v=randn(t,N);x=;for i=1:t;xt=filter(b0,a,v(i,:);x=x;xt;end;rg=;for i=1:t;rt1=E_r(x(i,:),N);rg=rg;rt1;end;rg=sum(rg)/t;k1=1:1:N;figure(1)plot(k1,rz,b-,k1,rg,r-);legend(真实值 ,估计值);xlabel(下标序列 );ylabel(自相关序列值 )grid on;%周期图法pt=;r
15、t=;figure(2)for i=1:tr1=E_r(x(i,:),N);rt=rt;r1;pt1=fft(r1,L)+conj(fft(r1,L)-r1(1);pt=pt;pt1;plot(w,abs(pt1);hold on;endxlabel(频率pi);ylabel(功率谱的幅值 );title(周期图法 50 次交叠图);grid on;hold offpt=sum(pt)/t;rt=sum(rt)/t;figure(3)plot(w,abs(pz),b-,w,abs(pt),r-);xlabel(频率/pi)ylabel(功率谱的幅值 )title(周期图法 50 次平均与真实功
16、率谱比较)legend(真实值 ,周期图法)grid on;%自相关法ag=;bg=;pf=;figure(4)for i=1:t;re=E_r(x(i,:),N);rp=re(1:(p+1);ag1,err=rtoa(rp);b1=sqrt(err);ag=ag;ag1;bg=bg;b1;% 冲击响应计算功率谱h1=dimpulse(b1,ag1,L);pf1=abs(fft(h1,L).2;pf=pf;pf1;plot(w,pf1);hold on;endxlabel(频率pi);ylabel(功率谱的幅值 );title(自相关法 50 次交叠图);grid on;hold offag=
17、sum(ag)/t;bg=sum(bg)/t;pf=sum(pf)/t;errfa=sum(ag-a).2)+(bg-1).2;figure(5)plot(w,abs(pz),b-,w,abs(pf),r-);xlabel(频率/pi)ylabel(功率谱的幅值 )title(自相关法 50 次平均与真实功率谱比较)legend(真实值 ,自相关法)grid on;%协方差法axg=;px=;figure(6)for i=1:t;axg1,errx=covm(x(i,:),p);bx=1;axg=axg;axg1;% 冲击响应计算功率谱hx1=dimpulse(bx,axg1,L);px1=a
18、bs(fft(hx1,L).2;px=px;px1;plot(w,px1);hold on;endxlabel(频率pi);ylabel(功率谱的幅值 );title(协方差法 50 次交叠图);grid on;hold offaxg=sum(axg)/t;px=sum(px)/t;errxa=sum(axg-a).2);figure(7)plot(w,abs(pz),b-,w,abs(px),r-);xlabel(频率/pi)ylabel(功率谱的幅值 )title(协方差法 50 次平均与真实功率谱比较)legend(真实值 ,协方差法)grid on;%burg 法abg=;pb=;figure(8)for i=1:t;gamma,errb=burg(x(i,:),p);abg1=gtoa(gamma);bb=1;abg=abg;abg1;% 冲击响应计算功率谱hb1=dimpulse(bb,abg1,L);pb1=abs(fft(hb1,L).2;pb=pb;pb1;plot(w,pb1);hold on;endxlabel(频率pi);ylabel(功率谱的幅值 );title(Burg 法 50 次交叠图);grid on;hold offabg=sum(abg)/t;
