现代谱估计法
（殷恒刚 107010254）
1. 现代谱估计简介
经典谱估计法可以利用FFT 计算，因而有计算效率高的优点，在谱分辨力要求不是太高的地方常用这种方法。
但频率分辨率地是经典谱估计的一个无法回避的缺点。
如周期图法在计算中把观测到的有限长的N 个数据以外的数据认为是零，而BT 法仅利用N 个有限的观测数据作自相关函数估计，实质上也就是假设除已知数据外的自相关函数全为零，这些显然都是与事实不符的。
为了克服以上缺点，人们提出了平均，加窗平滑等方法，在一定程度上改善了经典谱估计的性能。
但是，经典谱估计，始终无法解决，频率分辨率与谱估计稳定性之间的矛盾，特别是在数据记录长度比较短时，这一矛盾尤其突出。
现代谱估计理论也就是在这种背景下产生的，以1967年Burg 提出的最大熵谱分析法为代表的现代谱估计法，不认为在观察到的N 个数据以外的数据全为零。
因此克服了经典法的这个缺点，提高了谱估计的分辨率。
后来发现线性预测自回归模型法(简称AR 模型法)与Burg 的最大熵谱分析法是等价的，它们都可归结为通过Yule-Walker 方程求解自回归模型的系数问题。
目前常用的求自回归模型系数的算法有三种：①为Levinson 递推算法；②为Burg 递推算法；③为正反向线性预测最小二乘算法。
2.现代谱估计的三种模型
由信号与系统相关知识可知，任何具有有理功率谱密度的随机信号都可以看成是由一白噪声激励一物理网络所形成。
如图一所示。
我们可以先假设一个模型，然后根据已记录数据估计参数值，这样就不用假设N 以外的所有数据全为零，这就克服了经典谱估计的缺点。
图1
一个系统的Z 域传递函数的一般形式如下：
00
()()
b
a n j
j
j n i i
i b
z
Y z X z a z
-=-==
∑∑ (1.1)
参数建模的任务也就是如何确定阶数a n 和b n 以及系统数组(1,,)i a a i n = 和
(1,,)j b b i n = 。
一般情况下，阶数a n 和b
n 已经给定，所以主要任务是求出传递
函数的系数。
对于线性系统的模型，可以分为三种参数模型，即AR 模型（自回归模型）,MA 模型（滑动平均模型）和ARMA 模型（自回归滑移平均模型）。
这三个模型是平稳随机信号的三种标准线性模型。
现分别介绍如下。
（1） AR 模型 该模型除b 0=1外，
(1,,)
j b b i n = 的取值均为零。
其系统函数为
1
1
()1p
k
k
k H z a
z
-==
+
∑
其差分方程为
1
()()()
p
k k x n w n a x n k ==--∑，可知其功率谱密度为
)
()()(1
2
-=
Φz A z A z w
x σ
所以可推出其功率谱为
2
2)
(1)(ω
σ
ωj w
x e
A P =
AR 模型又称为p 阶自回归模型，简称AR 模型。
其系统函数只有极点，没有零点，因此也称为全极点模型。
（2） MA 模型 该模型除a
=1外，
(1,,)
i a
a i n = 的取值均为零。
其差分方程为
()()
q
k
k x n b
w n k ==
-∑，其系统函数为
()()()
q
k
k
k X z H z b
z W z -==
=
∑
故其功率谱密度为
)
()()(1
2
-=Φz
B z B z w x σ，所以可得到其功率谱为：
2
2
)
()(ω
σ
ωj w
x e B P =
该模型称为q 阶滑动平均模型，简称MA 模型。
其系统函数只有零点，没有极点，因此也成为全零点模型。
（3） ARMA 模型
ARMA 是AR 与MA 模型的组合。
该模型中，a 0=1，b 0=1，(1,,)i a a i n = 和
(1,,)
j b b i n = 不全为零。
其系统函数为
1()1q
k
k k p
k
k
k b z
H z a
z
-=-==
+
∑
∑
其差分方程为
1
()()()
q
p
k
k k k x n b
w n k a x n k ===
---∑∑，这里)(n w 是零均值，功率为2w
σ
的白噪声。
a k
就是自回归参数，b k
叫做动平均参数。
可得
)
()()()()(1
12
--=Φz
A z A z
B z B z w
x σ
，故
(1.5.2)
其功率谱为2
2)
()()(ω
ωσ
ωj j w
x e
A e
B P =。
由其系统函数可知，ARMA 模型的系统函数同时
存在极点和零点。
基于模型的功率谱估计方法大体上可按以下几个步骤进行： （1）模型的选择；
（2）根据已知数据进行模型参数的估计；
（3）将模型参数代入功率谱的计算公式得到功率谱估值。
对于模型的选择， Wold 分解定理告诉我们，如果功率谱是连续的，则任何MA 过程或ARMA 过程也可以用一个无限阶次的AR 过程表示，也就是说，如果选择了一个不合适的模型，只要模型的阶数足够高，仍然能够比较好地逼近实际的随机过程。
但这个定理并不是说明模型的选择不重要，例如，如果实际过程是AR 或MA 过程，而我们选用的是ARMA 模型，则我们所要的估计参数会增多。
实际上，由于可用的数据总是有限的，不论采用何种估计方法，待估计的参数越多，估计的精度就越差。
以上三种参数模型中，AR 模型得到了广泛的应用，其原因为：
一.任意ARMA 或MA 信号模型都可以用无限阶或阶数足够大的AR 模型来表示。
二.AR 模型的参数计算是线性方程，比较简便。
同时对窄的频谱有良好的适应性。
三.AR 模型在作谱估计时，由于具有递推特性，所以所需的数据较短，很适合表
示窄的功率谱；而MA模型表示窄谱时，一般来说，需要的参数较多；ARMA 模型表示时，虽然需要的参数较少，但参数估计的算法是非线性方程组，其运算远比AR模型复杂。
四.实际的物理系统往往是全极点系统。
所以研究有理分式传递函数的模型，主
要研究AR模型。
3. 运用matlab实现谱估计
本文采用Brug法进行功率谱估计。
Matlab源程序（my_zuoye3.m）如下：clear;
clc;
N=1024;
Nfft=1024;
n=[0:N-1];
randn('state',0);
wn=randn(1,N);
xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(20)*sin(2*pi*0.3*n)+wn;
subplot(2,1,1)
plot(n,xn)
xlabel('时间(s)')
ylabel('幅度')
[Pxx1,f]=pburg(xn,15,Nfft,1);%用Brug法进行功率谱估计，阶数为15，点数为1024
Pxx1=10*log10(Pxx1);
hold on;
subplot(2,1,2);plot(f,Pxx1,'g');
xlabel('Frequency');
ylabel('Power Spectrum (dB)');
title('Burg法 order=15,N=1024');
grid on;
N=1024， Nfft=1024时谱估计的结果如下图所示。
200
400
6008001000
1200
-20-10010
20时间(s)
幅度
0.05
0.1
0.15
0.2
0.250.30.35
0.4
0.45
0.5
-2002040
60Frequency
P o w e r S p e c t r u m (d B )
Burg 法 order=15,N=1024
Burg 递推算法因为不需要计算自相关，而是用使前向与后向预测误差能量之和最小的方法求出模型的参数，避免了由有限个数据估计自相关函数的计算及矩阵求逆，并且对于短的时间序列x(n)仍能得到较正确的估计，因此得到普遍应用。