报告2AR过程的线性建模与功率谱估计.docx
《报告2AR过程的线性建模与功率谱估计.docx》由会员分享，可在线阅读，更多相关《报告2AR过程的线性建模与功率谱估计.docx（15页珍藏版）》请在冰豆网上搜索。
报告2AR过程的线性建模与功率谱估计
实验报告2——AR过程的线性建模与功率谱估计
一、实验目的
1.理解AR过程的产生机理，复习实验1估计自相关序列的方法。
2.利用估计出的自相关序列来求解信号的功率谱，即用周期图法来估计功率谱。
3.分别采用自相关法（Yule-Walker法），协方差法，Burg法，修正协方差法来估计功率谱，并与周期图法进行比较，分析性能孰优孰劣。
4.学习matlab在数字信号处理中的应用。
二、实验过程和分析
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个样点。
用randn（1,N）产生单位方差高斯白噪声v（n），用v（n）激励滤波器产生AR（4）过程，即用x=filter（b,a,v）产生x（n），b是滤波器分子系数，这里为b（0）=1，a是滤波器分母系数，a=[1-2.76073.8106-2.65350.9238]。
（b）计算其自相关序列的估计
，并与真实的自相关序列值相比较。
结论：
真实自相关序列与估计出的序列对比如上图所示，两者在100点之前的形状相似，幅度有一定差异，而且估计出的自相关序列有较大的波动，这是因为估计的点数较少，使得估计精度不够，另外，估计自相关序列的下标越大，用来估计的点数就越少，因而后面的估计值是很不精确的。
因为假设n>=500处的数据为0，所有估计自相关不精确，误差较大。
（c）将
的DTFT作为x（n）的功率谱估计，即：
。
结论：
周期图法求得的功率谱误差较大，波动也较大，两个谱峰的分辨也并不清晰，并不精确，这是由于周期图法对自相关序列加了一个矩形窗。
（d）利用所估计的自相关值和Yule-Walker法（自相关法），估计
和
的值，并讨论估计的精度。
a
（1）=2.5018a
（2）=-3.2213a（3）=2.0794a（4）=-0.7025b（0）=2.4198
对比实际结果，可知Yule-Walker法（自相关法）误差较大，估计并不理想。
（e）用（d）中所估计的
和
来估计功率谱为：
。
估计出的功率谱上图所示，难以看出有二个谱峰，这是由于自相关法在短数据量的情况下估计是不精确的。
（f）将（c）和（e）的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。
结论：
用Yule-Walker法估计得到的谱曲线光滑，但无法识别出有二个谱峰，且b（0）的估计明显是错的，使得幅度有较大的偏差；相关图的曲线不光滑，不好观察谱峰，但大致走向与真实谱是一致的。
可以发现，Yule-Walker法（参数法）比周期图法（非参数法）精度更高。
但这里Yule-Walker法得到的功率谱分辨率较低，这是由于对自相关序列加了一个矩形窗引起的。
（g）重复上面的（d）~（f），只是估计AR参数分别采用如下方法：
（1）协方差法；
（2）Burg方法；（3）修正协方差法。
试比较它们的功率谱估计精度。
协方差法
a
（1）=2.6843a
（2）=-3.6442a（3）=2.5026a（4）=-0.8668b（0）=1.5506
与Yule-Walker法估计出参数的误差相比，有非常大的改进，这是因为协方差法对自相关序列估计时不需要对数据加窗，因此，对短数据记录，协方差法一般可获得比自相关法更高的分辨率。
由上图可知，谱估计的幅度相差较大，但从形状和极点定位的角度而言，估计结果很好的，曲线光滑，且极点的定位非常准确，短观测数据也能对极点精确定位是协方差法的优点。
Burg方法
a
（1）=2.7872a
（2）=-3.8728a（3）=2.7115a（4）=-0.9467b（0）=1.3006
Burg算法最小化前向加反向预测误差的平方和来求解全极点参数，为了保证稳定性，这一最小化是相对于反射系数序贯地进行的，其估计精度低于协方差法。
此估计方法与Yule-Walker法估计出参数的误差相比，有非常大的改进，但与协方差法比较，精度没有其理想。
用估计出的参数来进行谱估计，如上图所示，蓝色虚线就是采用Burg法估计出的功率谱，红线是真实功率谱，由图可知，估计结果较好，曲线光滑，且极点的定位较为准确，Burg法是自相关法与协方差法的一个折中。
修正协方差法
a
（1）=2.7863a
（2）=-3.8716a（3）=2.7106a（4）=-0.9467b（0）=1.3006
修正协方差法是最小化前向加反射预测误差的平方和，峰值偏移对相位的敏感度较低，比Yule-Walker法精度高，与Burg法的精度相当，但比协方差法的精度要低一些。
用估计出的参数来进行谱估计，如上图所示，蓝色虚线就是采用修正协方差法估计出的功率谱，红线是真实功率谱，由图可知，估计结果较好，曲线光滑，且极点的定位较为准确。
结论：
协方差法和修正协方差法可以分辨出本题中两个频率，说明它们相比于自相关法有更高的分辨率。
这是因为它在形成的自相关估计是不需要对数据加窗。
Burg算法也没有对数据加窗，因此其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）是单位方差白噪声
%Pb（0）=1,a=[2.7607,-3.8106,2.6535,-0.9238]
clc;
clear;
closeall;
N=256;
M=1024;
B=[1];
A=[1-2.76073.8106-2.65350.9238];
a=[2.7607,-3.8106,2.6535,-0.9238];
b0=1;
nf=0:
1/1000:
1;
w=pi*nf;
v=randn（1,N）;
x=filter（B,A,v）;
figure;
stem（x）;
title（'AR过程生成序列x（n）'）;
%估计自相关
Nr=N;
r1=zeros（1,Nr）;
fork=0:
（Nr-1）
forn=1:
（N-k）
r1（k+1）=r1（k+1）+x（n）*x（n+k）;
end
r1（k+1）=r1（k+1）/N;
end
figure;
stem（r1,'rx'）;
holdon;
%计算真是自相关
r=ator（A,b0）;
fork=length（a）+2:
N
r（k）=a
（1）*r（k-1）+a
（2）*r（k-2）+a（3）*r（k-3）+a（4）*r（k-4）;
end
stem（r,'b*'）;
title（'自相关估计'）;
legend（'估计自相关','真实自相关'）;
gridon;
%利用估计的自相关求功率谱
P=zeros（1,M）;
P2=zeros（1,M）;
P_t=zeros（1,M）;
f=0:
1/（M/2-1）:
1;
Y=fft（x,M）;
P=（Y.*conj（Y））/N;
Py=10*log10（P）;
figure;
plot（f,Py（1:
length（f））,'--'）;
holdon;
%真实功率谱
forh=1:
length（nf）
temp=0;
fork=1:
5
temp=temp+A（k）.*exp（-j*（k-1）*w（h））;
end
P_a（h）=1/（temp.*conj（temp））;
end
P_a=10*log10（abs（P_a））;
plot（nf,P_a,'r'）;
title（'周期法获取的功率谱和真实功率谱的比较'）;
legend（'周期图功率谱','真实功率谱'）;
xlabel（'数字频率/\pi'）;
ylabel（'功率谱/dB'）;
gridon;
%利用所估计的自相关，估计AR参数a
（1）,a
（2）,a（3）,a（4）,b0
[est_a,epsilon]=rtoa（r1（1:
5））
est_b0=sqrt（epsilon）
forh=1:
length（nf）
temp=0;
fork=1:
5
temp=temp+est_a（k）.*exp（-j*（k-1）*w（h））;
end
P_acor（h）=epsilon/（temp.*conj（temp））;
end
P_acor=10*log10（abs（P_acor））;
figure;
plot（nf,P_acor,'--'）;
title（'利用Yule-Walker法估计功率谱'）;
xlabel（'数字频率/\pi'）;
ylabel（'功率谱/dB'）;
gridon;
figure;
plot（nf,P_a,'r'）;
holdon;
plot（nf,P_acor,'--'）;
title（'功率谱比较'）;
legend（'真实功率谱','Yule_Walker法'）;
xlabel（'数字频率（\pi）'）;
ylabel（'功率谱（dB）'）
gridon;
%利用协方差求功率谱
p=4;
[a_cov,err_cov]=covm（x,p）
forh=1:
length（nf）
temp=0;
fork=1:
5
temp=temp+a_cov（k）*exp（-j*（k-1）*w（h））;
end
P_cov（h）=x
（1）^2/（temp.*conj（temp））;
end
P_cov=10*log10（abs（P_cov））;
figure（6）;
plot（nf,P_a,'r'）;
gridon;
holdon;
plot（nf,P_cov,'--'）;
title（'功率谱比较'）;
legend（'真实功率谱','协方差法'）;
xlabel（'数字频率/\pi'）;
ylabel（'PowerSpectrum/dB'）;
%修正协方差法
[a_mcov,err_mcov]=mcov（x,4）
forh=1:
length（nf）
temp=0;
fork=1:
5
temp=temp+a_mcov（k）*exp（-j*（k-1）*w（h））;
end
P_mcov（h）=x
（1）^2/（temp.*conj（temp））;
end
P_mcov=10*log10（abs（P_mcov））;
figure;
plot（nf,P_a,'r'）;
gridon;
holdon;
plot（nf,P_mcov,'--'）;
title（'功率谱比较'）;
legend（'真实功率谱','修正协方差'）;
xlabel（'数字频率/\pi'）;
ylabel（'功率谱/dB'）;
%Burg法求功率谱
[g_Burg]=burg（x,4）;
a_Burg=gtoa（g_Burg）
forh=1:
length（nf）
temp=0;
fork=1:
5
temp=temp+a_Burg（k）*exp（-j*k*w（h））;
end
P_Burg（h）=1/（temp.*conj（temp））;
end
P_Burg=10*log10（abs（P_Burg））;
figure;
plot（nf,P_a,'r'）;
gridon;
holdon;
plot（nf,P_Burg,'--'）;
title（'功率谱比较'）;
legend（'真实功率谱','Burg法'）;
xlabel（'数字频率/\pi'）;
ylabel（'功率谱/dB'）;
升级会员 