AR过程的线性建模与功率谱估计.docx

1、AR过程的线性建模与功率谱估计AR过程的线性建模与功率谱估计 一、实验内容：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=256个样点。(b) 计算其自相关序列的估计，并与真实的自相关序列值相比较。(c) 将的DTFT作为x(n)的功率谱估计，即：。(d) 利用所估计的自相关值和Yule-Walker法(自相关法)，估计和的值，并讨论估计的精度。(e) 用(d)中所估计的和来估计功率谱为：。(f) 将(c)和(e)的两种功率谱估

2、计与实际的功率谱进行比较，画出它们的重叠波形。(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(n)的N=256个样点。计算其自相关序列的估计，并与真实的自相关序列值相比较。思路：计算真实的自相关值时，采用逆Levinson-Durbin递归方法，由a、b参数得到，其中为滤波器的阶数，再采用公式外推得到的自相关值；结果分析：仿真参数设

3、置：采样点数为256自相关序列的估计与真实自相关序列值的比较见图1，由图可知估计值与真实值存在一定的误差，但整体变化趋势相差不大。图1 自相关估计和实际比较问题二：将的DTFT作为x(n)的功率谱估计，即：。利用所估计的自相关值和Yule-Walker法(自相关法)，用所估计的和来估计功率谱为： 。将周期图法和自相关法的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。思路：实际功率谱， 可调用Matlab中的FFT算法得到；自相关序列的估计值采用公式得到结果分析：通过图4比较可知，周期图能辨认出两个峰值，而自相关法不能，说明周期图的分辨率大于自相关法。同时比较图2和图3可知，周期图的

4、方差大于自相关法。 图2 周期图法估计AR(4)过程的功率谱图3 自相关法估计AR(4)过程的功率谱 图4 周期图法和自相关法得到的功率谱问题三：重复上面的(d)(f)，只是估计AR参数分别采用如下方法：(1) 协方差法；(2) Burg方法；(3) 修正协方差法。试比较它们的功率谱估计精度。思路： 方法和问题二的思路相似，只是使用的方式不同而已结果分析：图5图7的左图部分分别为采用协方差法、Burg方法、修正协方差法进行功率谱50次估计的交叠图，右图部分给出了其整体平均及真实的功率谱。由这些图可以看出，对于这一AR（4）过程，这三种方法都可以分辨出两个峰值来，而且方差大小差不多。但是和图2和

5、图3进行综合比较得，除自相关法外，所有估计都能分辨出两个峰值，且峰值的位置大致相似。同时，周期图法的方差大于其它估计方法。图5 协方差法估计AR(4)过程的功率谱图6 Burg法估计AR(4)过程的功率谱图7 修正协方差法估计AR(4)过程的功率谱问题四：分别采用自相关法、协方差法、Burg方法、修正协方差法，估计和的值，并讨论估计的精度。思路：采用平方误差和的计算公式为，来衡量四中方法的估计精确度。结果分析：从表1可知，自相关法估计的参数与真实值相比相差较大， Burg方法和修正协方差法参数估计的性能相当，其中协方差法最精确。表1 各种方法估计的a参数和b参数a(1)a(2)a(3)a(4)

6、b(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三、实验深入当观测数据存在观测噪声时，即，其中是单位方差的白噪声，与不相关，考虑观测噪声对各种谱估计方法的影响。图8是在存在观测噪声情况下，周期图法、协方差法、Burg方法和修正协方差法对AR（4）过程的功率谱估

7、计，由图可知，观测噪声对各种谱估计方法的估计性能具有一定的影响，周期图法辨别能力相对比较强，其次是协方差法还可以辨别出两个峰值。但是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,

8、-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(

9、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:t r1=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=s

10、um(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); pf1=

11、abs(fft(h1,L).2; pf=pf;pf1; plot(w,pf1); 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).2)+(bg-1).2;figure(5)plot(w,abs(pz),b-,w,abs(pf),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(自相关法50次平均与真实功率谱比较)legend(真实值,自相关法)grid on;%协方

12、差法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).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

13、(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(频率

14、pi);ylabel(功率谱的幅值);title(Burg法50次交叠图);grid on;hold offabg=sum(abg)/t;pb=sum(pb)/t;errba=sum(abg-a).2);figure(9)plot(w,abs(pz),b-,w,abs(pb),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(Burg法法50次平均与真实功率谱比较)legend(真实值,Burg法)grid on;%修正协方差法axxg=;pxx=;figure(10)for i=1:t; axxg1,errxx=mcov(x(i,:),p); bxx=1; axxg=

15、axxg;axxg1; % 冲击响应计算功率谱 hxx1=dimpulse(bxx,axxg1,L); pxx1=abs(fft(hxx1,L).2; pxx=pxx;pxx1; plot(w,pxx1); hold on;endxlabel(频率pi);ylabel(功率谱的幅值);title(修正协方差法50次交叠图);grid on;hold offaxxg=sum(axxg)/t;pxx=sum(pxx)/t;errxxa=sum(axxg-a).2);figure(11)plot(w,abs(pz),b-,w,abs(pxx),r-);xlabel(频率/pi)ylabel(功率谱的

16、幅值)title(修正协方差法50次平均与真实功率谱比较)legend(真实值,修正协方差法)grid on;1、 计算真实的自相关值时，采用逆Levinson-Durbin递归方法，由a、b参数得到，其中为滤波器的阶数，再采用公式外推得到的自相关值；仿真参数设置：采样点数为256自相关序列的估计与真实自相关序列值的比较见图1，由图可知估计值与真实值存在一定的误差，但整体变化趋势相差不大。clear all;a1=2.7607;a2=-3.8106;a3=2.6535;a4=-0.9238;t=50;p=4;%a=1,-a1,-a2,-a3,-a4;N=256;rz=zeros(1,N); b

17、0=1; 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; v=randn(t,N); x=; for i=1:t; xt=filter(b0,a,v(i,:); x=x;xt;end; rg=; for i=1:t; rt=E_r(x(i,:),N); rg=rg;rt;end; rg=sum(rg)/t; k1=1:1:N; plot(k1,rz,b-,k1,rg,r-);legend(真实值,估计值); xlabel(下标序列); ylabel(自相关序列值)自相关

18、法 a 1.0773* 2.30912828300994 -1.22044171792557 0.387568716048877b .0819*Error 5.062586351911801 clearclearclear all;a1=2.7607;a2=-3.8106;a3=2.6535;a4=-0.9238;t=50;p=4;%a=1,-a1,-a2,-a3,-a4;N=256;b0=1;v=randn(t,N);x=;for i=1:t; xt=filter(b0,a,v(i,:); x=x;xt;end;for i=1:t;end;az=;ax=;ab=;ax=;b=;for i=1

19、:t; re=E_r(x(i,:),N); rp=re(1:(p+1);a11,err1=rtoa(rp); b1=sqrt(err1); a11=a11; az=az;a11; b=b;b1;endaz=sum(az)/t;b=sum(b)/t; err1=sum(az-a).2);%周期图法pt=;rt=;figure(1)for i=1:t r1=E_r(x(i,:),N); rt=rt;r1; p=fft(r1,L)+conj(fft(r1,L)-r1(1); pt=pt;p; plot(w,abs(p); hold on;endxlabel(频率pi);ylabel(功率谱的幅值);

20、title(周期图法50次交叠图);grid on;hold offpt=sum(pt)/t;rt=sum(rt)/t;figure(2)plot(w,abs(pz),b-,w,abs(pt),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(周期图法50次平均与真实功率谱比较)legend(真实值,周期图法)grid on;1.0760* 2.28961393515513 -1.19747620735719 0.374809926614702.021*9.29004642552237 %自相关法ag=;bg=;rf=;pf=;figure(1)for i=1:t; r

21、e=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=sum(ag)/t;bg=sum(bg)/t;pf=sum(pf)/t;errfa=sum(ag-a).2)+(bg-

22、1).2;figure(2)plot(w,abs(pz),b-,w,abs(pf),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(自相关法50次平均与真实功率谱比较)legend(真实值,自相关法)grid on;1 -2.74999799454203 3.78018793385172 -2.62169354886228 0.90922991193951610.00226336448829608%协方差法 axg=;rf=;px=;figure(1)for i=1:t; axg1,errx=covm(x(i,:),p); bx=1; axg=axg;axg1; %

23、冲击响应计算功率谱 hx1=dimpulse(bx,axg1,L); px1=abs(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(2)plot(w,abs(pz),b-,w,abs(px),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(协方差法50次平均与真实功率谱比较)legend

24、(真实值,协方差法)grid on;1 -2.73951319606315 3.75693194746984 -2.59988963513036 0.90078990111695610.00673267639550165 %burg法abg=;rf=;pb=;figure(1)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); ho

25、ld on;endxlabel(频率pi);ylabel(功率谱的幅值);title(Burg法50次交叠图);grid on;hold offabg=sum(abg)/t;pb=sum(pb)/t;errba=sum(abg-a).2);figure(2)plot(w,abs(pz),b-,w,abs(pb),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(Burg法法50次平均与真实功率谱比较)legend(真实值,Burg法)grid on;1 -2.74123330526539 3.75889893476294 -2.60029116603434 0.9001

26、2938421140310.00644343041433378 %修正协方差法axxg=;rf=;pxx=;figure(1)for i=1:t; axxg1,errxx=mcov(x(i,:),p); bxx=1; axxg=axxg;axxg1; % 冲击响应计算功率谱 hxx1=dimpulse(bxx,axxg1,L); pxx1=abs(fft(hxx1,L).2; pxx=pxx;pxx1; plot(w,pxx1); hold on;endxlabel(频率pi);ylabel(功率谱的幅值);title(修正协方差法50次交叠图);grid on;hold offaxxg=sum(axxg)/t;pxx=sum(pxx)/t;errxxa=sum(axxg-a).2);figure(2)plot(w,abs(pz),b-,w,abs(pxx),r-);xlabel(频率/pi)ylabel(功率谱的幅值)title(修正协方差法50次平均与真实功率谱比较)legend(真实值,修正协方差法)grid on;