AR过程的线性建模与功率谱估计.docx
《AR过程的线性建模与功率谱估计.docx》由会员分享,可在线阅读,更多相关《AR过程的线性建模与功率谱估计.docx(35页珍藏版)》请在冰豆网上搜索。
![AR过程的线性建模与功率谱估计.docx](https://file1.bdocx.com/fileroot1/2023-7/27/6ec87aad-7021-4df8-a044-f01064e52a2f/6ec87aad-7021-4df8-a044-f01064e52a2f1.gif)
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)的两种功率谱估计与实际的功率谱进行比较,画出它们的重叠波形。
(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参数得到
,
,
,
,其中
为滤波器的阶数,再采用公式
外推得到
的自相关值;
结果分析:
仿真参数设置:
采样点数为256
自相关序列的估计
与真实自相关序列值的比较见图1,由图可知估计值与真实值存在一定的误差,但整体变化趋势相差不大。
图1自相关估计和实际比较
问题二:
将
的DTFT作为x(n)的功率谱估计,即:
。
利用所估计的自相关值和Yule-Walker法(自相关法),用所估计的
和
来估计功率谱为:
。
将周期图法和自相关法的两种功率谱估计与实际的功率谱进行比较,画出它们的重叠波形。
思路:
实际功率谱
,
可调用Matlab中的FFT算法得到;自相关序列的估计值采用公式
得到
结果分析:
通过图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)过程的功率谱
图6Burg法估计AR(4)过程的功率谱
图7修正协方差法估计AR(4)过程的功率谱
问题四:
分别采用自相关法、协方差法、Burg方法、修正协方差法,估计
和
的值,并讨论估计的精度。
思路:
采用平方误差和的计算公式为
,来衡量四中方法的估计精确度。
结果分析:
从表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.0010
Burg法
-2.7465
3.7737
-2.6151
0.9063
1
0.0033
修正协方差法
-2.7472
3.7747
-2.6157
0.9063
1
0.0031
三、实验深入
当观测数据存在观测噪声时,即
,其中
是单位方差的白噪声,与
不相关,考虑观测噪声对各种谱估计方法的影响。
图8是在存在观测噪声情况下,周期图法、协方差法、Burg方法和修正协方差法对AR(4)过程的功率谱估计,由图可知,观测噪声对各种谱估计方法的估计性能具有一定的影响,周期图法辨别能力相对比较强,其次是协方差法还可以辨别出两个峰值。
但是Burg方法和修正协方差法不能。
周期图法它的性能主要依赖于数据序列长度N,N越大分辨率越好。
当增加频域采样点数时,Burg方法和修正协方差法也能辨认出两个峰值,但频域采样点数的增加意味着计算量的增加。
仿真参数设置:
噪声误差为1,采样点数为64,频域采样点数为128。
图8加噪声情况下的功率估计方法比较
附录:
MATLAB代码
clearall;
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);
fork=(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=[];
fori=1:
t;
xt=filter(b0,a,v(i,:
));
x=[x;xt];
end;
rg=[];
fori=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('自相关序列值')
gridon;
%周期图法
pt=[];
rt=[];
figure
(2)
fori=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));
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('周期图法50次交叠图');
gridon;
holdoff
pt=sum(pt)/t;
rt=sum(rt)/t;
figure(3)
plot(w,abs(pz),'b--',w,abs(pt),'r--');
xlabel('频率/\pi')
ylabel('功率谱的幅值')
title('周期图法50次平均与真实功率谱比较')
legend('真实值','周期图法')
gridon;
%自相关法
ag=[];
bg=[];
pf=[];
figure(4)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('自相关法50次交叠图');
gridon;
holdoff
ag=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('真实值','自相关法')
gridon;
%协方差法
axg=[];
px=[];
figure(6)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('协方差法50次交叠图');
gridon;
holdoff
axg=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('真实值','协方差法')
gridon;
%burg法
abg=[];
pb=[];
figure(8)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('Burg法50次交叠图');
gridon;
holdoff
abg=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法')
gridon;
%修正协方差法
axxg=[];
pxx=[];
figure(10)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('修正协方差法50次交叠图');
gridon;
holdoff
axxg=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('功率谱的幅值')
title('修正协方差法50次平均与真实功率谱比较')
legend('真实值','修正协方差法')
gridon;
1、计算真实的自相关值时,采用逆Levinson-Durbin递归方法,由a、b参数得到
,
,
,
,其中
为滤波器的阶数,再采用公式
外推得到
的自相关值;
仿真参数设置:
采样点数为256
自相关序列的估计
与真实自相关序列值的比较见图1,由图可知估计值与真实值存在一定的误差,但整体变化趋势相差不大。
clearall;
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);
>>b0=1;
>>rz(1:
(p+1))=ator(a,b0);
>>fork=(p+2):
N;
rz(k)=-sum(a(2:
p+1).*seqreverse(rz(k-p:
k-1)));
end;
>>v=randn(t,N);
>>x=[];
>>fori=1:
t;
xt=filter(b0,a,v(i,:
));
x=[x;xt];
end;
>>rg=[];
>>fori=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('自相关序列值')
自相关法
a1.0773********2.30912828300994-1.220441717925570.387568716048877
b.0819********
Error5.062586351911801
>>clear
clear
clearall;
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=[];
fori=1:
t;
xt=filter(b0,a,v(i,:
));
x=[x;xt];
end;
fori=1:
t;
end;
az=[];
ax=[];
ab=[];
ax=[];
b=[];
fori=1:
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];
end
az=sum(az)/t;
b=sum(b)/t;
>>err1=sum((az-a).^2);
%周期图法
pt=[];
rt=[];
figure
(1)
fori=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));
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('周期图法50次交叠图');
gridon;
holdoff
pt=sum(pt)/t;
rt=sum(rt)/t;
figure
(2)
plot(w,abs(pz),'b--',w,abs(pt),'r--');
xlabel('频率/\pi')
ylabel('功率谱的幅值')
title('周期图法50次平均与真实功率谱比较')
legend('真实值','周期图法')
gridon;
1.0760********2.28961393515513-1.197476207357190.374809926614702
.021*********
9.29004642552237
>>%自相关法
ag=[];
bg=[];
rf=[];
pf=[];
figure
(1)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('自相关法50次交叠图');
gridon;
holdoff
ag=sum(ag)/t;
bg=sum(bg)/t;
pf=sum(pf)/t;
errfa=sum((ag-a).^2)+(bg-1).^2;
figure
(2)
plot(w,abs(pz),'b--',w,abs(pf),'r--');
xlabel('频率/\pi')
ylabel('功率谱的幅值')
title('自相关法50次平均与真实功率谱比较')
legend('真实值','自相关法')
gridon;
1-2.749997994542033.78018793385172-2.621693548862280.909229911939516
1
0.00226336448829608
%协方差法
>>axg=[];
rf=[];
px=[];
figure
(1)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('协方差法50次交叠图');
gridon;
holdoff
axg=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('真实值','协方差法')
gridon;
>>
1-2.739513196063153.75693194746984-2.599889635130360.900789901116956
1
0.00673267639550165
>>%burg法
abg=[];
rf=[];
pb=[];
figure
(1)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('Burg法50次交叠图');
gridon;
holdoff
abg=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法')
gridon;
1-2.741233305265393.75889893476294-2.600291166034340.900129384211403
1
0.00644343041433378
>>%修正协方差法
axxg=[];
rf=[];
pxx=[];
figure
(1)
fori=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);
holdon;
end
xlabel('频率\pi');
ylabel('功率谱的幅值');
title('修正协方差法50次交叠图');
gridon;
holdoff
axxg=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('真实值','修正协方差法')
gridon;