AR过程的线性建模与功率谱估计.docx
《AR过程的线性建模与功率谱估计.docx》由会员分享，可在线阅读，更多相关《AR过程的线性建模与功率谱估计.docx（35页珍藏版）》请在冰豆网上搜索。
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;
升级会员 