AR过程的线性建模与功率谱估计一、实验内容：AR 过程的线性建模与功率谱估计。考虑 AR 过程：x(n) = a(1)x(n -1) + a (2) x( n - 2) + a (3) x(n - 3) + a(4) x(n - 4) + b(0) v( n)v(n) 是单位方差白噪声。(a)取 b(0)=1, a(l)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生 x(n)的 N=256 个 样点。 计算其自相关序列的估计产(k)，并与真实的自相关序列值相比较。x(c)将r (k)的DTFT作为x(n)的功率谱估计，即:xP (ej)= 艺 r (k)e-jk = I X (ej) bxxNk =- N+1(d)利用所估计的自相关值和Yule-Walker法(自相关法)，估计a(1), a(2), a(3), a和b(0)的| b(0) |21 + 才 a(k)e-jk值，并讨论估计的精度。(e)用(d)中所估计的a(k)和b(0)来估计功率谱为：P (ej)=xk=1(f)将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。(g)重复上面的(d)，只是估计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(n)的 N=256 个样 点。计算其自相关序列的估计r (k)，并与真实的自相关序列值相比较。x思路：计算真实的自相关值时，采用逆Levinson-Durbin递归方法，由a、b参数得到r (0)，xr (1),，r (p)，其中P为滤波器的阶数，再采用公式r (k)=Ya (l)r (k-1)外推得到xxxp xl=1k p的自相关值；结果分析：仿真参数设置：采样点数为 256自相关序列的估计产(k)与真实自相关序列值的比较见图1,由图可知估计值与真实值存x在一定的误差，但整体变化趋势相差不大。-600-800050100200250300150下标序列08 6 4002值列序关相自00图 1 自相关估计和实际比较问题二：将r (k)的DTFT作为x(n)的功率谱估计，即: xP (ej) = r (k)e- jkw =I X (ej) b。xxNk=- N+1I b (02) I 2利用所估计的自相关值和Yule-Walker法(自相关法)，用所估计的a(k)和b(0)来估计功率P (ejra )=x谱为：1 + 才 a(k)e-jkk=1将周期图法和自相关法的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波 形。迟 r (k)e- jkI k=o x 丿-r(0),x思路：实 际 功 率 谱 P(ej)= r(k)e-jk =r(k)e-jk+xxxk=-N+1k=0N-1r(k)e-jk 可调用 Matlab 中的 FFT 算法得到；自相关序列的估计值采用公式xk=0r (k)= 1 迓x(n) x* (n - k )得到xNn=1结果分析：通过图4比较可知，周期图能辨认出两个峰值，而自相关法不能，说明周期图的分辨 率大于自相关法。同时比较图2和图3可知，周期图的方差大于自相关法。频率兀图2周期图法估计AR(4)过程的功率谱频率兀图3自相关法估计AR(4)过程的功率谱频率/兀图4 周期图法和自相关法得到的功率谱问题三：重复上面的(d)(f),只是估计AR参数分别采用如下方法：(1)协方差法；(2) Burg 方法； (3) 修正协方差法。试比较它们的功率谱估计精度。思路：方法和问题二的思路相似,只是使用的方式不同而已结果分析：图5图7的左图部分分别为采用协方差法、Burg方法、修正协方差法进行功率谱50 次估计的交叠图,右图部分给出了其整体平均及真实的功率谱。由这些图可以看出, 对于这一 AR (4)过程，这三种方法都可以分辨出两个峰值来，而且方差大小差不多。但是和图 2 和图 3 进行综合比较得,除自相关法外,所有估计都能分辨出两个峰值,且 峰值的位置大致相似。同时，周期图法的方差大于其它估计方法。图5协方差法估计AR(4)过程的功率谱图6 Burg法估计AR(4)过程的功率谱2.52修正协方差法50次交叠图0.50 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 频率兀值幅的谱率功0值幅的谱率功图 7 修正协方差法估计 AR(4) 过程的功率谱问题四：分别米用自相关法、协方差法、Burg方法、修正协方差法，估计a(l), a(2), a(3), a(4) 和b(0)的值，并讨论估计的精度。思路：采用平方误差和的计算公式为=亡(a(i)-a(i)，来衡量四中方法的估计精确度。i=1结果分析：从表1可知，自相关法估计的参数与真实值相比相差较大，Burg方法和修正协方差 法参数估计的性能相当，其中协方差法最精确。表 1 各种方法估计的 a 参数和 b 参数aaa(3)ab(0)平方误差真实值-2.76063.8106-2.65350.92381自相关法-2.07032.2862-1.18890.36763.10779.6969协方差法-2.75273.7901-2.63180.913810.0010Burg 法-2.74653.7737-2.61510.906310.0033修正协方差法-2.74723.7747-2.61570.906310.0031三、实验深入当观测数据存在观测噪声时，即y(n) = x(n)+ w(n)，其中w(n)是单位方差的白噪声, 与v (n )不相关，考虑观测噪声对各种谱估计方法的影响。图 8 是在存在观测噪声情况下，周期图法、协方差法、 Burg 方法和修正协方差法对 AR(4)过程的功率谱估计，由图可知，观测噪声对各种谱估计方法的估计性能具有一 定的影响，周期图法辨别能力相对比较强，其次是协方差法还可以辨别出两个峰值。但是 Burg 方法和修正协方差法不能。周期图法它的性能主要依赖于数据序列长度 N，N 越大分辨率越好。当增加频域采样 点数时，Burg方法和修正协方差法也能辨认出两个峰值，但频域采样点数的增加意味着 计算量的增加。仿真参数设置：噪声误差为 1，采样点数为 64，频域采样点数为 128。图 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).* 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=;rt=;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次平均与真实功率谱比较)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);pfl=abs(fft(hl,L).人2;pf=pf;pf1;plot(w,pfl);hold on;endxlabel(频率 pi);ylabel (功率谱的幅值)；title(自相关法50次交叠图)； grid on;hold offag=sum(ag)/t; bg=sum(bg)/t; pf=sum(pf)/t;errfa=sum(ag-a).A2)+(bg-1).A2; 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=abs(fft(hx1,L).A2; px=px;px1;plot(w,px1);hold on;endxlabel(频率 pi);ylabel (功率谱的幅