谱估计（复习大纲）

第一章 经典谱估计

经典谱估计方法是以傅里叶变换为基础的方法，主要有两类：周期图法和布莱克曼—图基法（简称BT法，又称为谱估计的自相关法）。这两类方法都与相关函数有着密切的联系，由维纳——欣钦定理可知，功率谱和相关函数之间的关系是一对傅里叶变换，因而可以从观测数据直接估计相关函数，根据估计出来的相关函数，求它的傅立叶变换，就可以得到功率谱的估计值。

一、 相关函数和功率谱

若 mx(n)?mx?常数，rxx(n1,n2)?rxx(n1?n2)即rxx(k)?E[x(n?k)x*(n)] 则称{x(n)}为广义平稳序列。

若{x(n)}和{y(n)}均为广义平稳序列，且rxy(n1,n2)?rxy(n1?n2)即

rxy(k)?E[x(n?k)y*(n)]，则称{x(n)}和{y(n)}为广义联合平稳序列。

广义平稳随机序列{x(n)}的相关函数rxx(k)和它的功率谱密度Pxx(?)之间是傅立叶变换对的关系，即

??Pxx(?)?rxx(k)?1?rxx(k)ek????j?kd? (1.6)

2?????Pxx(?)ej?kd? (1.7)

这一关系式常称为维纳——欣钦定理。

由自相关函数和功率谱密度的定义，不难得出它们的一些基本性质，主要有：

1、当{x(n)}为复序列时，rxx(?k)?rxx*(k)；若{x(n)}为实序列，则相关函数为偶函数，即rxx(?k)?rxx(k)。

2、相关函数的极大值出现在k?0处，即rxx(k)?rxx(0)。

3、若x(n)含有周期性分量，则rxx(k)也含有同一周期的周期性分量，否则，当k??时，rxx(k)?0。

4、当x(n)为实序列时，Pxx(?)为非负实对称函数，即Pxx(??)?Pxx(?)和

Pxx(?)?0。

5、平稳随机序列{x(n)}的自相关函数rxx(0)是实的且为正，而且对任一a(n)序列和任一M，自相关函数（ACF）满足：

M?12M?1M?1E{?a*(n)x(n)}?n?0??a*(m)a(n)rm?0n?0xx(m?n)?0 （1.10）

这个函数称为半正定的。

自相关函数（ACF）和互相关函数（CCF）的z变换定义为：

??Pxx(z)??rxx(k)zk??????k （1.11）

Pxy(z)??rxy(k)zk????k （1.12）

1212若令??2?f,f为归一化频率，频率区间??f?为基本周期。则式（1.6）的自功

率谱密度和式（1.8）互功率谱密度又可分别表示为

??Pxx(f)??rxx(k)ek????j2?fk （1.13）

??Pxy(f)??rxy(k)ek????j2?fk （1.14）

Pxx(f)是实的，且非负。

当一平稳随机序列{x(n)}通过一个脉冲响应为h(n)的线性非时变系统时，其输出序列

{y(n)}也是一平稳随机序列。它的自相关函数为

ryy(k)?h(k)?h(?k)?rxx(k) (1.16)

*Pyy(z)?H(z)H*(1z*)Pxx(z) （1.17） 1)?H()。令z?exp(j?)?exp(j2?f)，得到相应的功率z*z2若h(n)为实系统，则H*(谱表达

1Pyy(?)?H(?)Pxx(?) (1.26a)

或 Pyy(f)?H(f)Pxx(f) (1.26b) 上述关系对以后讨论谱估计问题是很有用的。

ryy(0)?12?2????P(?)d???1/2?1/2Pyy(f)df为输出过程{y(n)}的平均功率。

经常遇到的一种过程是离散白噪声，它的自相关函数（ACF）定义为：

rxx(k)??x?(k) （1.29）

2其中?(k)是离散冲激函数。这就是说，各样本之间彼此是不相关的。

Pxx(f)??1/2?1/2rxx(k)e?j2?f2df??x （1.30）

这表明它在各频率上是完全平坦的。换句话说，白噪声的所有频率分量均具有相同的功率。

二、 相关函数的估计

1、 自相关函数的各态历经性

一般说来，严格各态历经过程允许我们用时间平均来代替系综平均（集合平均或统计平均）。

用时间平均作为广义平稳随机过程均值的估计。

limM???x(n)?E{x(n)}??x （1.34）

2M?1n??MM1M lim1M??2M?1n??M?x*(n)x(n?k)?E[x*(n)x(n?k)]?rxx(k) （1.36）

2、 相关函数的估计

我们实际所能得到的随机序列的样本数总是有限的，由有限个样本通过某种运算求出的序列的均值和自相关函数统计特征值叫做它们的估计值。下面讨论随机序列有限个样本的相关函数的估计问题。

设{x(n),n?0,1,?,N?1}为实随机序列{x(n)}的一批样本，共有N个值。有时简称之为长度为N的随机序列x(n)。

方法一：根据假定的自相关函数的各态历经性（或遍历性），可用下式估计它的自相关函数，即

?xx(k)?r1N?k1N?k1N?kN?k?1?x(n?n?0N?k?1k)x(n) （1.38）

?xx(k)]?E[r?E[x(n?n?0N?k?1k)x(n)]

??rn?0xx(k)?rxx(k) （1.39）

?xx(k)}?0， 当N??时， var{r?xx(k)是相关函数rxx(k)的无偏估计且是渐近一致的，即当k为有限值时，r?xx(k)是因此rrxx(k)的一致估计。

方法二：有限长度序列{x(n),n?0,1,?,N?1}的相关函数rxx(k)的另一种估计方法

可表示为

?xx(k)?r1NN?k?1?x(n?n?0k)x(n) （1.42）

?xx(k)]?E[r1NN?k?1?n?0E[x(n?k)x(n)]?N?kNrxx(k) （1.43）

可见，它是相关函数rxx(k)的有偏估计。但是，当N??，估计值是渐近无偏的。

?xx(k)}?0， 当N??时，var{r?xx(k)也是rxx(k)的一致估计。 即（1.42）式的r（1.42）式所定义的相关函数取傅立叶变换求功率谱估计时，在计算上有某些方便之处，以后的讨论中，如不作特别申明，将采用这种有偏估计表示式求相关函数的估计式。

三、功率谱密度的另一个定义

可以证明，功率谱密度（PSD）的一个近似等效的定义是

Pxx(f)?limE{M??12M?1M2?x(n)exp(?j2?fn)n??M} （1.45）

上式定义的PSD与维纳一辛钦定理

??Pxx(f)??rk???xx(k)exp(?j2?fk) （1.46）

是等效的。

由（1.45）式和由（1.46）式维纳一辛钦定理给出的PSD的等效定义将作为经典谱估计方法的基础。

四、 周期图

周期图谱估计器定义为：

?(f)?PPER1NN?12?x(n)exp(?j2?fn)n?0N?1xxk??(N?1) （1.48）

可以证明，周期图等于估计出的自相关序列的傅里叶变换，或

?P(f)?PER?r?(k)exp(?j2?fk)

?xx(k)是有偏自相关函数的估计值，定义为 其中r1NN?1?k?xx(k)?r?x(n)x(n?n?0k)

周期图的期望值是

?E[P(f)]??{wB(k)rxx(k)}?PER?1/2?1/2WB(f??)Pxx(?)d? （1.49）

式中WB(f)是巴特利特（Bartlett）窗（三角形窗）的傅里叶变换。

?由式（1.49）可知周期图的均值E[P(f)]是真实PSD和巴特利特窗傅里叶变换的卷PER积，在平均意义上得到真实功率谱密度（PSD）的平滑形式。因此对有限记录数据，周期图

一般有偏的，但是当N??时，它是无偏的。

N???limE[P(f)]?Pxx(f) PER这是由于WB(f)收敛到狄拉克?函数。

22对于高斯白噪声的特殊情况，rxx(k)??x?(k),Pxx(f)??x，结果为

?(f)]?? E[PPER2x?Pxx(f) （1.50）

对于白噪声情况，即使有限记录数据，周期图也是无偏的。

?(f)]??对于白高斯过程，E[PPER2x?Pxx(f)，方差

sin2?Nf22?var{P(f)}?P(f)[1?()] （1.52） PERxxNsin2?f对任何非白过程，只要记录数据足够长，

?(f)]?P2(f)[1?(sin2?Nf)2] var[PPERxxNsin2?f对于不靠近0或?12的频率，且N??时，上式近似地退化为：

2?var[PPER(f)]?Pxx(f) （1.53）

可以看出，方差与数据长度N无关，即方差不随N的增大而减小至零。由此可得出一个重要的看法：周期图估计器是不可靠的，因为标准离差（方差的平方根）和均值一样大，因而周期图不是一致估计而其均值近似地等于要估计的量值。

上述的重要结论表明，我们不能寄希望于直接用周期图方法获得良好的谱估计，必须采用适当的修正措施减小估计方差，才能使之成为一种实用的方法。

五、周期图法改进措施

1、 数据窗

周期图法只用了N个样本，这可以看作是用一长度为N的矩形窗函数与原来无限长的序列相乘的结果，我们知道，时域中两函数相乘对应于频域中它们的傅立叶变换的卷积。由

此可以想到，当用周期图方法作谱估计时，它的谱分辨率约与N成反比，且和信号本身的特征，例如信噪比等无关。此外，如果序列是由多个正弦波信号组成的，而各分量强度不等，

则弱信号分量可能淹没在强信号谱的旁瓣中而无法发现。这种所谓信号能量（向旁瓣）泄漏现象如果不设法消除，也将妨碍周期图谱估计法的应用。

周期图改进的方法之一是将长度为N的序列x(n)乘以同一长度的数据窗w(n)。 数据加窗后的周期图谱估计值的数学期望值等于谱的真实值与窗谱函数的平方的卷积。显然它不等于功率谱的真实值，因而是有偏估计。

若序列为单频信号，则Pxx(f)为?函数，这样，数据加窗后的谱估计值的均值与窗谱函数的平方形状相同，因此选用低旁瓣的数据窗可使得杂散响应减少。但旁瓣的降低必然使主瓣加宽，而且降低了分辨率。

数据加窗后，周期图谱估计值的方差大于或近似等于谱估计值的平方，且不随数据长度的增大而减小到0。

从以上的分析可知，数据加窗用于周期图谱估计可以降低谱估计值的旁瓣，但要降低谱估计的分瓣率。而用数据加窗的办法不能减小估计方差。

2、 平均周期图

为了改进周期图的统计特性，可以近似地用对一组周期图进行平均的方法完成期望运算。假定在区间0?n?L?1上有k组独立记录数据，并且都是同一随机过程的现实。平均周期图估计器定义为

?P(f)?AVPER1K??Pm?0K?1(m)PER(f) （1.65）

?(m)(f)是第m个数据组的周期图 其中PPER?(m)(f)?PPER1LL?12m?xn?0(n)exp(?j2?fn) （1.66）

1/2?E[PAVPER(f)]???1/2WB(f??)Pxx(?)d?

但是方差将减少K倍。

1??(m)(f)] （1.68） var[P(f)]?var[PAVPERPERK谱估计器具有1/L周/样本的分辨率。

显然，为了获得最大分辨率，L将尽可能选得大一些或就选L?N?1，此时，平均周期图就成为标准周期图估计。但有时为了减小方差，应该把K?N/L选得大一些，或等效地把L取得小一些。因为这两个目的不可能同时实现，所以必须调整L在方差和偏差之间进行折衷。

减小偏差的一种很有用的技术是在谱估计之前对数据预白化。

六、 BT谱估计器

?(f)?PBTN?1?w(k)r?k??(N?1)xx(k)e?j2?fk （1.74）

其中w(k)是实序列，称为滞后窗。根据滞后窗的性质，（1.74）式可以写成：

?(f)?PBTM?w(k)r?k??Mxx(k)e?j2?fk （1.75）

上式称为布莱克曼—杜基（Blackman and Tukey 1958，BT）谱估计器，如果w(k)?1,k?M?N?1，则它等效于周期图。

BT谱估计器有时称为加权协方差估计器。对自相关函数估计器加权，将减小谱估计器的方差，但以增大偏差为代价（白噪声过程除外，对于它使用任何滞后窗偏差均为0）。有很多可利用的窗，

对于实数据过程，由（1.76）式可得BT谱估计器的均值

?(f)]?E[PBT1/2??1/2?(?)]d? W(f??)E[PPER对于长记录数据计算，可得：

?(f)]?E[PBT?1/2?1/2W(f??)Pxx(?)d? (1.77)

这时因为当N??时，周期图是无偏的。

这个均值是真实PSD的一种模糊形式，这就是说W(f)起到一个谱窗的作用。对于多数的滞后窗，BT估计器能够分辨大约1/M周/样本或相当于谱窗主辨带宽的谱的细节。和周期图一样，预白化可以用来减小这个偏差。

?(f)的方差为： PBT?(f)]?Pxx(f)var[PBTN2?1/2?1/2W(?)d??2Pxx(f)N2M?wk??M2(k) (1.79)

?(f)的均值(1.77)式和方差(1.79)式可清楚表明需要对偏差和方差折衷选择。偏差要研究PBT小，M应选得足够大，因为这样才能使谱窗起到狄拉克?函数的作用。另一方面方差要小，M就应选得尽量小。通常建议M的最大值取N/5。

七 本章小结

根据相关函数和功率谱密度函数之间为傅立叶变换对的维纳-欣钦关系，可从相关函数估计值经傅立叶变换后得到功率谱的估计值。从原理上说这种方法是直截了当的，可是由此得出的谱估计值不仅有偏，而且方差很大，它还不随数据长度的增大而减小到零。因此必须采取其他措施改善它的估计性能，才能使得这种方法有实用价值。常用的两种办法是相关函数加窗和分段平均周期图法，这两种办法都有实用价值，它们都采用了FFT算法，使得计算量大大减小。

归结起来，BT方法和周期图方法的主要优点是；（1）计算量小：（2）功率谱估计值正比于正弦波信号的功率。（3）对于长记录数据，是一种良好的实用模型。这们的主要缺点有：（1）频率分辨率（区分两个邻近频率分量的能力）不高。这是因为它们的频率分辨率约为数据长度的倒数，且与数据的特征或其信噪比无关，而实际应用中一般不可能获得很长的数

据记录；（2）经典谱估计方法在工程中都是以离散傅立叶变换为基础的，它隐含着对无限长数据序列进行加窗处理（加了一个有限宽的矩形窗）。矩形窗的频谱主瓣不是无限窄的，且有旁瓣存在，这将导致能量向旁瓣中“泄漏”，主瓣变得模糊不清。严重时，会使主瓣产生很大失真，甚至主瓣中的弱分量被旁瓣中的强泄漏所淹盖。（3）弱信号被强信号的旁瓣淹没；（4）需采用某种平滑或平均措施以改善谱估计的统计特性；（5）某些加窗的相关函数会使功率谱估计值出现负值。

第二章 谱估计的参数模型方法

以参数模型（ARMA、AR、MA）为基础的谱估计方法一般按下列三个步骤进行： 第一步：为被估计的随机过程确定或选择一个合理的模型。这有赖于对随机过程进行的理论分析和实验研究；

第二步：利用可提供的数据样本估计假设模型的参数。这涉及到对各种算法的研究。通常，模型参数的数据量比观测数据的数据量少很多，因此，为数据压缩创造了条件；

第三步：将估计出的模型参数代入模型的理论功率谱密度公式中，从而得到谱估计值。

在选择模型时要考虑模型能表示谱峰、谱谷和滚降的能力。根据极、零点对谱的不同贡献，自回归（AR）模型比较适合于具有尖峰但没有深谷的谱；滑动平均（MA）模型比较适合于具有深谷但无尖峰的谱；而自回归滑动平均（ARMA）模型对两种极端情况都能够表示。

一、有理传递函数模型（ARMA、AR、MA）

（一）定义

1、自回归滑动平均（ARMA）模型：（零、极点模型）

pq差分方程： x(n)???a(k)x(n?k)?k?1?b(k)u(n?k) (2.1 )

k?0系统函数：

q H(z)?B(z)A(z)?b(k)z?k?op?k （2.2）

?k?a(k)zk?o假定H(z)的极点在单位圆内，就可保证H(z)是一个稳定和因果的滤波器。若没有这个假定，由（2.1）式给出的x(n)就不能用来描述广义平稳随机过程。

令Puu(z)?Z[ruu(k)],Pxx(z)?Z[rxx(k)]，则在线性滤波器输出端自相关函数的z变换Pxx(z)与输入的z变换Puu(z) 间的关系如下：

Pxx(z)?H(z)H(*****1z*)Puu(z)?B(z)B(1/z)A(z)A(1/z)Puu(z) （2.3）

当(2.3)式沿单位圆计算时，z?ej2?f,?12?f?12，上式成为功率谱密度Pxx(f)。通常

假定驱动过程是均值为零和方差为?2的白噪声序列，则白噪声序列的功率谱密度为?2，即Puu(f)??2。

ARMA输出过程的功率谱密度为： PARMA(f)?Pxx(f)??其中A(f)?A(ej2?f),增益都可以并入?2。

2、滑动平均（MA）过程（全零点模型）：

q2B(f)A(f)2 （2.4）

B(f)?B(ej2?f一般假定a(0)?1,b(0)?1，因为任何滤波器)。

差分方程： x(n)? ?b(k)u(n?k) （2.5）

k?0MA输出过程的功率谱密度为：

PMA(f)??2B(f)2q2?j2?fk??2?b(k)ek?0 (2.6)

3、自回归（AR）过程（全极点模型）

p 差分方程： x(n)???a(k)x(n?k)?u(n) （2.7）

k?1AR过程的功率谱密度为： PAR(f)?Pxx(f)??22?p?1?22?j2?fk （2.8）

A(f)?a(k)ek?1（二）时间序列模型间的关系

沃尔德分解定理[wold,1954]将AR，MA和ARMA三个模型联系起来。基本定理表明，任何广义平稳随机过程都可以分解成完全随机的分量和确定性分量。

沃尔德定理的推论是说，如果功率谱是纯连续的，任何ARMA或AR过程都可以用无限阶的唯一的MA模型表示，而柯尔莫哥洛夫定理[Kolmogorov,1941]表明，任何ARMA或MA过程都可以用无限阶的AR过程表示。这些定理是很重要的，因为若我们在三种模型中选择了不正确的模型，只要选取的模型阶次足够高，仍然可以得到一个合理的近似表示。

重点掌握ARMA(1，1)过程与AR(∞)或MA(∞)模型之间的参数转换关系。

二、模型参数与自相关函数的关系

★ 本节重点要求掌握AR过程模型参数与自相关函数的关系即尤拉—沃克

（Yule-walker）方程的推导过程。

人们总希望得到模型参数和自相关函数之间明显的关系式。这些关系式可以从模型z变换的z反变换得到，也可以从时间序列模型直接导出。

1、ARMA过程模型参数与自相关函数的关系 ???? rxx(k)??????pq?kxx?a(l)rl?1p(?l)??(k?l)2?hl?0*(l)b(l?k),k?0,1,?,q （2.15）

k?q?1?a(l)rl?1xx,应当注意到，ARMA过程的参数和自相关函数之间的关系是非线性的，给出自相关函数后，

q?k必须解一组非线性方程来求出模型参数，这是因为有?h*(l)b(l?k)项的缘故。

l?02、AR过程模型参数与自相关函数的关系 ???? rxx(k)??????p?a(l)rl?1pxx(k?l),2k?1 （2.16）

xx?a(l)rl?1(?l)??,k?0方程(2.16)式称为尤拉—沃克方程，它定义了AR过程的参数和自相关函数之间的线性关系。如果给定自相关函数，就可以解一组线性方程来确定AR参数。

尤拉—沃克（Yule-walker）方程用矩阵形式表示为：

rxx(?1)?rxx(?(p?1))??a(1)??rxx(0)?rxx(1)???????rxx(1)rxx(0)?rxx(?(p?2))a(2)rxx(2)???? （2.17） ??????????????????????a(p)r(p?1)r(p?2)?r(0)r(p)?xx????xx??xx?????xx???????????RxxH 注意，自相关矩阵Rxx是埃尔米特的(Rxx?Rxx)和托布列兹的，而且矩阵是半正定的。

如果加入包括?的方程，方程（2.17）式也可变为： ?rxx(0)?rxx(1)? ????rxx(p)rxx(?1)rxx(0)?rxx(p?1)???rxx(?(p?1))??????rxx(0)??rxx(?p)2?1????????a(1)0? （2.18） ??????????????a(p)0??????2上式可以从（2.16）式得到。这个形式在以后是很有用的。

3、MA过程模型参数与自相关函数的关系

?2q?k*???b(l)b(l?k), rxx(k)??l?0?0,?k?0,1,?,qk?q?1 （2.19）

可以看出，MA过程的参数和自相关函数之间的关系是非线性的。

三、ARMA、AR和MA过程举例

在这一节中，把前面得到的结果应用于产生低阶（即p和q较小）ARMA，AR和MA过程的自相关函数和功率谱。即实AR(1)过程、实ARMA(1,1)过程、实MA(1)过程，实

MA(2)过程自相关函数和功率谱。

★ 重点要求掌握在已知实AR(1)过程的参数a(1)和白噪声的方差?2时自相关函数和功率谱的求解。

四、根据功率谱或自相关函数确定模型参数

若AR、ARMA过程其H(z)的极点位于z平面单位圆内，则H(z)是—稳定的和因果的滤波器。若H(z)是非因果稳定的则其极点必须在单位园外。不论哪种情况，稳定性十分关键，否则x(n)的方差将会无穷大。

前向AR模型（因果模型）：

p x(n)???a(k)x(n?k)?u(n) （2.41）

k?1后向AR模型（非因果模型）：

p x(n)???a(k)x(n?k)?u(n) (2.42)

k?1*对于给定的功率谱或等效地对于给定的自相关函数，可能得出两个性质不同的AR时间序列模型。若由Pxx(z)求前向模型或后向模型系数，需进行谱分解：用单位圆内的Pxx(z)的极点确定前向模型AR参数，用单位圆外的Pxx(z)的极点确定后向模型AR参数。 前向和后向模型在以后讨论各种AR估计方法时会用到。

概括地说，如果给出功率谱或自相关函数，则ARMA，AR或MA参数不是唯一的，但是也可以是唯一的，这只有通过限制该过程是因果的和/或可逆的方法来实现。为了保证因果性和可逆性，A(z)和B(z)均要求是最小相位的。如果不是这样，可用谱分解的方法将滤波器变换成最小相位形式。

第三章 AR模型谱估计

一、AR过程的线性预测

以下假定x(n)是一个AR(p)过程，线性预测的问题，是根据观测数据组

?x(n?1),x(n?2),?,,x(n?p)?（即前p个样本）来预测未观测到的样本x(n)。假定预

测器是过去样本的线性组合，即

p?(n)????kx(n?k) （3.1） xk?1?1,?2,?3??选择预测系数?p?使预测误差e(n)的功率

?(n)|2] (3.2) ??E[|e(n)|]2?E[|x(n)?x达到极小。经过推导?k?a(k),k?1,2,?,p和?min??2。这就是说，最佳线性预测系数恰好是AR模型参数，并且最小预测误差功率恰好是激励噪声方差。然而仅当AR过程的阶次与线性预测器的阶次相同时，这才是正确的。另外，预测误差e(n)正是u(n)。

u(n) H(z)?输入白噪声序列 1A(z)?11?a(1)z?1?a(2)z?2??a(p)z?p x(n)

(a) p阶AR模型

x(n) H

(b) 预测误差滤波器 图3.1 线性预测的滤波解释

p?1(z)?A(z)?1?a(z)z?1???a(p)z?p 预测误差(白噪声) e(n)=u(n) 最佳预测误差滤波器（PEF）是AR过程的逆滤波器A(z)?1??a(k)zk?1?k，它的输入

为x(n)，输出为预测误差。显然，预测误差时间序列是u(n)，因此可以把预测误差预波器看作是白化滤波器。事实上，根据以上的分析，AR（P）的过程可另外定义为

?(n)?u(n) （3.7） x(n)?xp?(n)???a(k)x(n?k)是基于前p个样本的最佳一步线性预测。 式中xk?1二、Levinson-Durbin算法（列文森-德宾算法）

如果我们希望不仅确定p阶线性预测器，而且要确定p?1,p?2,?,1阶线性预测器，一种可能是对不同假定的模型阶次去解上述尤拉—沃克方程式。结果将得出一组预测系数

{a2(1),a2(2)},?,{ap(1),ap(2),?,ap(p)}，{a1(1)}，其中aj(i)是j阶线性预测的第i个

系数。显然，对i?1,2,?,p，有ap(i)?a(i)。这个过程虽然直截了当，但计算繁琐，没有必要全部进行。另一种方法是由k?1阶预测器到k阶预测器递推地更新。

（一）、前向线性预测

?k?1(n)是x(n)的基于前k?1现在若已知x(n?1),?,x(n?(k?1))等k?1个值，令x个样本即x(n?1),?,x(n?k?1)最佳k?1阶线性预测器，或者

k?1?k?1(n)???ak?1(i)x(n?i) （3.8） xi?1?(n)的下标表示在预测中所用到的前面样本的数。ak?1(i)称为前向线性预测系数。在得到x真值x(n)后，即可得前向线性预测误差为：

k?1efk?1?k?1(n)?x(n)?(n)?x(n)?x?ai?1k?1(i)x(n?i) （3.9）

?k?1(n)预测x(n)的误差ekf?1(n)的均方值为最现在我们来求ak?1(i)的最佳值，使得由x小。即

?fk?1?E[efk?1(n)] （3.10）

2为最小。由复梯度公式，经过推导可得复信号情况的前向预测正交原理：

E[ekf?1(n)?x*(n?i)]?0,若x(n)为实信号，其前向预测正交原理为

E[ek?1(n)x(n?i)]?0,ffi?1,2,?,k?1 （3.20a）

i?1,?,k?1） （3.20b）

将ek?1(n)的表达式（3.9）代入前向预测正交原理公式（3.20a）中，得到

k?1E{[x(n)??al?1k?1(l)x(n?l)]x*(n?i)}?0i?1,2,?,k?1 （3.21）

根据自相关函数定义，由（3.21）式可得

k?1rxx(i)??al?1k?1(l)rxx(i?l)?0, i?1,2,?,k?1 （3.22）

对应于最佳预测系数的最小预测误差为

k?1 ?fk?1?rxx(0)??al?1k?1 (l)rxx(?l) （3.23）

式（3.22）和式（3.23）写成矩阵形式就是

?rxx(0)?r(1)?xx????rxx(k?1)rxx(?1)rxx(0)?rxx(k?2)????frxx(?(k?1))??1???k?1??????rxx(?(k?2))ak?1(1)0? （3.24） ??????????????????rxx(0)a(k?1)0????k?1???上述方程称为尤拉—沃克(Yule-Walker)方程。这是k个方程的方程组。已知

rxx(0),rxx(1),?rxx(k?1)即可解出ak?1(1),ak?1(2),?,ak?1(k?1),?k?1，得到k?1阶最佳

f前向预测误差滤波器和相应的最小预测误差。

（二）、后向线性预测

k?1阶前向预测是根据x(n?1),x(n?2),?,x(n?(k?1))预测x(n)。相应地，k?1阶后向预测是由x(n?(k?1)),x(n?(k?2)),?,x(n?1)预测x(n?k)。

对x(n?k)的后向线性预测可表示为：

k?1?k?1(n?k)???bk?1(k?1?i)x(n?k?i) (3.25b) xi?1相应的后向线性预测误差为：

k?1ebk?1?k?1(n?k)?(n?1)?xk?1(n?k)?x?bi?0k?1(k?1?i)x(n?k?i) (3.26)

上式中bk?1(k?1)?1。

与前向预测进行相似的推导，可得最佳后向预测系数满足正交原理：

E{ek?1(n?1)x(n?k?i)}?0,b*i?1,2,?k?1 (3.28)

经推导可得最佳后向预测系数与前向预测系数之间有下列关系：

bk?1(k?1?l)?ak?1(l),*l?1,2,?,k?1 (3.30)

将（3.30）式代入（3.26）式可得后向线性预测误差

k?1ebk?1(n?1)??ai?0*k?1(i)x(n?k?i) (3.31)

最小后向预测均方误差与前向预测误差相等，可令其为?k?1，则

k?1?bk?1??fk?1??k?1?rxx(0)??ai?1k?1(l)rxx(?l) (3.33)

（三）、Levinson-Durbin算法

fb定义k?1阶前向预测误差ek?1(n)和后向预测误差ek?1(n?1)之间的相关函数为：

?k?E[ek?1(n)(ek?1(n?1))] (3.35)

fb*对于最佳预测系数，根据正交原理可得：

k?1?k?rxx(k)??al?1k?1(l)rxx(k?l) (3.36)

求解最佳线性预测系数的Levinson-Durbin算法如下： 对于零阶递推：

e0(n)?x(n)

fk0?0

?0?E[e0(n)]?E[x(n)]?E[x(n)x(n)]?rxx(0)

f22*对于1阶递推：

e1(n)?x(n)?a1(1)x(n?1)

fa1(1)?k1??rxx(1)rxx(0)2

2?1?(1?a1(1))?0?(1?a1(1))rxx(0)

对于k?2,3,?,p，递推为

k?1ak(k)???krxx(k)???*?al?1k?1(l)rxx(k?l)?k?1?k?1 (3.43)

i?1,2,?,k?1 (3.44)

ak(i)?ak?1(i)?ak(k)ak?1(k?i)2?k?(1?ak(k))?k?1 (3.45)

反射系数由kk?ak(k)给出。

kk在线性预测中起着重要的作用。

kk???k?k?1f??E[ek?1(n)(ek?1(n?1))]E[ebfk?1fb*(n)]*2

??E[ek?1(n)(ek?1(n?1))]E[ebk?1(n?1)]fb2

??cov(ek?1(n),ek?1(n))var(efk?1(n)var(ebk?1 （3.46）

(n?1))式中cov表示协方差，由上式可知kk?1。很容易看出，反射系数是前向和后向预测误差之间相关系数的负值，所以kk又称偏相关系数（PARCOR）。

如果过程的确是一个AR（p）过程，则对于k?1,2,?,p，有ap?1(k)?ap(k)和

ap?1(p?1)?kp?1?0。一般来说，对AR（p）过程，若当k?p时ak(k)?kk?0，而

对于k?p有?k??p，这就是说，当模型的阶次等于或大于真实模型的阶次时，模型中激励噪声的方差是一个常数。因此，?k维持不变这一点似乎是正确模型的良好表示。由ak(k)?kk?1这一性质可得出?k??k?1，这意味着?k在正确模型阶次处首先达到它的

最小值。如果对于某个k值出现kk?1的情况，递推必须终止，因为?k?0，然而这种情况只有过程仅由k个正弦信号组成时才会发生。

三、AR过程的三种等效表示方法

一个AR(p)过程具有下述三种等效表示 ①{a(1),a(2),?a(p),?} ② {rxx(0),k1,k2,?kp} ③ {rxx(0),rxx(1),?,rxx(p)} 这三种表示法之间可以互相转化。

2

1、{rxx(0),rxx(1),?,rxx(p)}转化到{a(1),a(2),?,a(p),?2}或{rxx(0),k1,?kp}用Levinson-Durbin算法（式（3.43）-（3.45）及相应的初始条件）实现。其中a(k)?ap(k),k?1,2,?,p,?2??p。

2、{rxx(0),k1,k2,?kp}转换到{a(1),a(2),?,a(p),?2}或{rxx(0),rxx(1),?,rxx(p)} （1） {rxx(0),k1,k2,?,kp}转换到{a(1),a(2),?,a(p),?2}

ak(k)?kkk?1,2,?,p

*ak(i)?ak?1(i)?kkak?1(k?i),p2i?1,2,?k?1,k?1,2,?,p

?2?rxx(0)?(1?kii?1) (3.50)

a(i)?ap(i),i?1,2,?p

（2） {rxx(0),k1,k2,?,kp}转换到{rxx(0),rxx(1),rxx(p)}

?0?rxx(0),?k?(1?kk)?k?1,2k?1,2,?p (3.51a)

??k1rxx(0)?rxx(k)??k?1???ak?1(l)rxx(k?l)?kk?k?1,?l?1k?1k?2,3,?,p (3.51b)

3、{a(1),a(2),?,a(p),?2}到{rxx(0),rxx(1),?,rxx(p)}或{rxx(0),k1,k2,?,kp}的转换要求有一个步降算法（step-down），它把AR参数转变为反射系数。

ak?1(i)?ak(i)?ak(k)ak(k?i)1?ak(k)2*i?1,2,?,k?1k?p,p?1,?,2 （3.52）

上式是步降（SD）方法的第一部分。它使我们可由k阶预测误差滤波器系数求得(k?1)阶预测误差滤波器系数。在这个过程中，反射系数以及所有低阶预测误差滤波器都被求出。 即由

?p??,ap(k)?a(k)*2k?1,2,?,p

i?1,2,?k?1k?p,p?1,?2ak?1(i)?ak(i)?ak(k)ak(k?i)1?ak(k)2

?k?1??k1?ak(k)2,k?p,p?1,?,1

2可完成{a(1),a(2),?,a(p),?}到{rxx(0),k1,k2,?kp}的转化。

最后，如果{rxx(0),rxx(1),?,rxx(p)}是要得到的，可在使用SD方法之后，再利用

(3.51b)式求解。

四、格型滤波器

★ 要求掌握推导格型滤波器的关系式及画出其实现结构。

利用列文森算法推导中得出的结果，可以推导出预测误差滤波器（PEF）的另一种等效形式。所谓的格型滤波器就是把反射系数作为滤波器的系数。

k阶前向预测误差：

efk??n???n??x?n??xkx(n)??al?1k(l)x(n?l)

k阶前向预测误差递推式：

b(n?1) (3.54) ekf(n)?ekf?1(n)?kkek?1k阶后向预测误差：

ke(n)?x(n?k)?bk?ai?1*k(i)x(n?k?i) (3.56)

k阶后向预测误差递推式：

ek(n)?ek?1(n?1)?kkek?1(n) (3.55)

bb*f格型滤波器的关系式：

ffb ek(n)?ek?1(n)?kkek?1(n?1)

bb*f ek(n)?ek?1(n?1)?kkek?1(n)

fb其中eo(n)?eo(n)?x(n)。

p格型滤波器如图3.5所示，与标准的预测误差滤波器A(z)?1??ak?1k(k)z?k等效。格型滤

波器有个有用的性质：它的系数的幅度总是小于1，因为ki?1,i?1,2,?p。

格型滤波器的关系式，不仅对AR（p）过程，而且对任何过程都是正确的。但是，如果过程就是AR（p）过程，则对于i?p?1有ki?0。这是因为后向预测误差在同一时刻它们不相关。

E[ei(n)*?ej(n)]?0,i?j (3.57)

bb ebk?1(n?1)2?efk?1(n?1)2?efk?1(n)2??k?1 (3.58)

bf var(ek(n))?var(ek(n))??k (3.59)

由此可知，格型滤波器将输入数据串x(n),x(n?1),?,x(n?k)转换成相互正交的数

bbb据串e0(n),e1(n),?,ek(n)。这种转换就是熟知的Gram-Schmidt正交化转换。

图3.5 用格型滤波器实现预测误差滤波器

五、AR谱估计器的特性

（一） 隐含的自相关函数延拓

AR谱估计的高分辨特性是由于隐含着测定出的自相关函数的延拓。

?xx(k)是 隐含自相关函数r?rxx(k),??xx(k)??pr?(l)r?xx(k?l),???a?l?10?k?pk?p （3.61）

?xx(k)值，我们并没有认为它们等于零，而应该按（3.61）式它说明对于在k?p范围内的r进行外推。

（二）、AR谱估计与最大熵谱估计等效

最大熵估计（MESE）的基础是由一段已知的自相关函数样本外推出未知的自相关函数样本（Burg,1975），这样，可以减小由自相关函数截断（或加窗）所造成的功率谱估计的特征模糊问题。

如果{rxx(0),rxx(1),?,rxx(p)}已知，则问题在于如何外推求得

{rxx(p?1),rxx(p?2),?}才能保证整个自相关函数矩阵是正定的。伯格（Burg）证明了选

择这样的外推方法才是最合理的：外推后的自相关函数序列所对应的时间序列应当具有最大熵。这意味着，在具有已知的p+1个滞后的自相关函数样本值对应的所有时间序列中，该时间序列将是最随机的或最不可预测的时间序列，换句话说，它的功率谱将是最平坦或白化的。通过由这样的外推得到的自相关序列求出的谱称为最大熵谱估计。

若随机信号x(n)是由零均值与单位方差正态白噪声u(n)激励因果最小相位线性系统的输出，这时x(n)为平稳高斯过程，其功率谱为Pxx(?)或Pxx(f)，且有

???lnP?xx(?)d??? （3.67）

则熵率hx与功率谱Pxx(?)的关系为

hx?ln2?e?14????ln??Pxx(?)d? （3.68）

谱熵定义为

H[Pxx(?)]?12?1??ln?Pxx(?)d? (3.69a)

或

H[Pxx(f)]??21?2lnPxx(f)df （3.69b）

最大熵谱分析法（MESE）是针对时域无限正态随机序列或正态随机过程情况，在给定自相关函数约束条件下，使时域熵率或谱熵为最大，这时求得的功率谱Pxx(f)便是谱估计结果。

最大熵谱估计为：

PMESE(f)?Pxx(f)?1??p22?j2?fk (3.77)

?a(k)ek?1其中{a(1),a(2),?a(p),?2}是利用已知的自相关函数样本通过解尤拉—沃克方程求出的。由于具有{rxx(0),rxx(1),?,rxx(p)}的知识，因此最大熵谱估计器等价于AR谱估计器，同时也等价于预测误差滤波器法。然而这个等价性仅对高斯随机过程和已知自相关函数样本的情况成立。

（三） AR谱估计等效于最佳白化处理

AR(p)参数就是预测误差功率最小时p阶线性预测器的系数。预测误差滤波器是白化

滤波器，它去掉了AR过程的相关性。

AR谱估计得出的结果可看成是最佳白化处理的结果。最大熵谱估计等效于使谱平坦度最大。这样，最大熵谱估计在固定的{rxx(0),rxx(1),?,rxx(p)}的约束条件下具有最大的谱平坦度。

六、AR模型参数提取方法

在实际应用中，常需根据信号的有限个取样值来估计AR模型的参数，应用较多的方法有Yule-Walker法或自相关法、协方差法、修正协方差法、伯格（Burg）法。以上方法都可以用由时间平均代替集合平均的最小平方准则推导得到。对于长记录数据，它们的特性是相似的，而对于短记录数据，各种方法的特性之间仍存在差别。通常，即使是对于短记录数据，AR谱估计的质量也还是相当好的。而且AR谱估计的计算量也比较适中。对于很多方法来说，为了实现“快速”算法，都尽量利用所求解方程的结构，以便降低需要的计算量。 AR谱估计器既能以AR参数估计为基础，也能以反射系数的估计为基础，这些方法均利用

使预测误差功率极小化的准则来估计参数。

1、自相关法

假定观察到的数据是?x(0),x(1),?,x(N?1)?。自相关法，或者有时称为尤拉—沃克法。估计AR参数是利用使预测误差功率的估计值为极小的方法，该预测误差功率为：

?? ?1N?2?n???ep(n)f2?1N?p?n???x(n)??a(k)x(n?k)k?1 （3.86）

为使预测误差功率极小，这可以利用求复梯度的方法进行，得

1N??p?n???[x(n)??a(k)x(n?k)]xk?1*(n?l)?0,l?1,2,?,p （3.87）

用矩阵形式，这组方程为：

?xx(0)?r??(1)r?xx????xx(p?1)?r?xx(?1)r?xx(0)r??xx(p?2)r?????xx(?(p?1))?r??xx(?(p?2))r?????xx(0)r??xx(1)??(1)??r?a?????xx(2)?(2)ra? （3.88） ????????????????xx(p)?a(p)???r

其中

?1N?1?k*x(n)x(n?k),k?0,1,2,?,p???xx(k)??Nn?0 r （3.89）

*?rk??(p?1),?(p?2),?,?1??xx(?k),*?xx(?k)?r?xx?xx(k)是有偏自相关函数的估计器。式（3.88）中的矩阵是埃尔米特（r(k))和r托布列兹的，并且可以证明也是正定的。另外，取名尤拉—沃克法是因为自相关法与使用有偏自相关函数估计器的尤拉—沃克方程等效。可用列文森递推来解该方程，并且根据最小相位理论保证估计出的极点在单位圆内。

?min，由下式给出： 求出的白噪声方差?2的估计值是?p? ?22?xx(0)??r?(k)r??ak?1xx(?k) （3.90）

?也可以由列文森递推的最后一步求得，即p阶预测误差功率，或者写成另一种形式 ?p? ?2?xx(0)??ri?1?(1?ki2) （3.91）

其中k?i是在列文森递推中得出的第i阶反射系数的估计值。利用列文森递推的自相关法由式（3.88）和式（3.90）给出。

用自相关法得出的谱估计值分辨率比用其它将要讨论的方法差。

2、协方差法

协方差法估计使用预测误差功率

???1NN?12??Pep(n)f?1N?pN?1p2?n?px(n)??a(k)x(n?k)k?1 （3.92）

n?p极小的方法求得。注意，协方差法和自相关法之间的唯一差别是预测误差功率估计中求和的范围不同。

使用复梯度法，得出AR参数估计的一组方程

1N?N?1p*[x(n)??a(k)x(n?k)]x?pn?pk?1(n?j)?0,j?1,2,?,p （3.93）

令

cxx(j,k)?1N?N?1?xpn?p*(n?j)x(n?k)

则式（3.93）可写为：

p?(k)c?ak?1xx(j,k)??cxx(j,0),j?1,2,?,p

AR参数估计值即为上述方程的解。

白噪声方差的估计值为

p? ?2?min?cxx(0,0)????(k)c?ak?1xx(0,k) （3.95）

*式(3.94)中的矩阵是埃尔米特的(cxx(k,j)?c(j,k)和半正定的。可用乔里斯基分解法去解

这个方程组。式(3.94)中的矩阵C不是托布利兹矩阵，不能用列文森-德宾递推算法求解

AR(p)参数。利用协方差法估计出的极点不能保证在单位圆内。

3、修正协方差法

修正协方差法用使估计的前向和后向预测误差功率的平均值为极小的方法估计AR参数，也就是

?? ?12?(?f?b) （3.98） ??其中

? ?f?1N?1N?1p2?p?n?0x(n)??a(k)x(n?k)k?1p （3.99a）

2n?pN?1?p?b? ?N?px(n)??ak?1*(k)x(n?k) （3.99b）

和协方差法一样，仅对哪些用到的观测数据样本的预测误差求和。应当指出，对任何一组

a(k)值，前向和后向预测误差估计值由于求和范围的缘故将有所不同，这个步骤等效于前

向和后向AR模型的预测误差组合在一起的方法。

为了使式（3.98）极小化，利用复梯度关系，可得

pN?1*N?1?p

?(k)[?x(n?k)x?ak?1n?pN?1(n?l)??xn?0*(n?k)x(n?l)]

N?1?p* ??[?x(n)x(n?l)?n?p?xn?0*(n)x(n?l)] （3.101）

令

cxx(j,k)?12(N?p)N?1N?1?p*[?x(n?j)x(n?k)?n?p?n?0x(n?j)x(n?k)] (3.102)

*则式（3.101）可写为：

p?(k)c?ak?1xx(l,k)??cxx(l,0),l?1,2,?,p （3.103）

白噪声方差的估计值是

??2?min???12(N?p)N?1p[?(x(n)?n?p?(k)x(n?k))x?ak?1*(n)

N?1?pp??n?0(x(n)?*?(k)x?ak?1*(n?k))x(n)]

其中利用了式（3.100），最后得

p??2?cxx(0,0)??(k)c?ak?1xx(0,k) (3.105)

可以看出，除自相关估计器cxx(j,k)的定义外，修正协方差法与协方差法相同。式（3.104）

*中的矩阵是埃尔米特的(cxx(k,j)?cxx(j,k))和正定的（纯正弦信号情况除外），因此可用

乔里斯基分解（平方根法）来解该线性方程组。

修正协方差法可以得到高分辨率的统计稳定的谱估计值。对于由白噪声中的正弦信号组成的数据，该估计器无谱峰偏移及谱峰分裂现象。

4、伯格（Burg）法

伯格法又可称反射系数法。伯格法是先估计反射系数，然后利用列文森递推得到AR参数的估计值。用递推方法对不同阶的预测是通过使预测误差功率的估计值达到极小，从而得到反射系数的估计值的。

伯格提出用使前向和后向预测误差功率估计的平均值达到极小的方法来估计反射系数?的有约kk。在这里，预测误差功率估计仍然用时间平均来代替集合平均。伯格法是一个?束极小化方法，有约束的或递推的极小化方法一般不能得到全局极小。因此，为得到kk的估计，使下式达到极小

?k??12?kf???kb) (3.106) (??kf和??kb仅是反射系数k的函数。 ?总的说来，AR参数估计的伯格方法是

?xx(0)?r1NN?1?n?0x(n)

2?0?r?xx(0) (3.118) ??0f(n)?x(n)eb?0e(n)?x(n)n?1,2,?,N?1 n?0,1,?,N?2

对于k?1,2,?,p，有

N?1b*?kf?1(n)e?k?2?e(n?1)?1n?kN?1??kk

(n)2b?k?e(n?1))?12??(en?kfk?1?k???(1?kk2?k?1 （3.119） )??a??k?1(i)?k?*(k?i),i?1,2,?k?1?akk?1?k(i)??a

??kk,i?k??1(1)?k1） （如果k?1，则a?e?kf(n)?e?kf?1(n)?k?b(n?1),ekk?1?(n)?e?ebkbk?1n?k?1,k?2,?,N?1n?k,k?1,?,N?2?*e?f(n),(n?1)?kkk?1 (3.120)

?p}。 ?p(1),a?p(2),?,a?p(p),?所要求出的估计值为{a??1的缘故。 伯格法所估计出的极点在单位圆上或单位圆内。这是因为有特性kk一般说来，对真正是AR过程的数据，伯格算法能得出精确的AR谱估计。但是对于正弦信号的数据，就有一定的困难。伯格算法有谱线分裂现象，且谱峰的位置与相位密切相关。为了减小对相位的依赖关系，可以利用修正反射系数估计

N?1?w?kk?k?1(n)e?k?1(n?1)?2?wk(n)en?kN?1fb* (3.122)

2b?k?e(n?1))?12?wn?kk?kf?1(n)(n)(e这里wk(n)是适当选择的具有非负权的窗。这个修正反射系数估计最初由Burg在1975年提出，已经证实可以减小相位的影响。

反射系数的伯格估计只是大量的具有最小相位性质的佑计中的一个。另一个估计由Itakura和Saito在1971年提出，用时间平均来代替集合平均，得：

N?1?I?kk??en?kN?1fk?1?k?1(n?1)(n)e2N?1b* （3.123）

b?ke(n?1)?12?n?k?kf?1(n)e?n?k可以证明

?B?k?I （3.124） kkk?B表示由（3.117）式给出的伯格估计，还有很多其它可能的估计方法，实际上，它其中kk们在谱估计性能方面的差别似乎很小。

七、AR模型阶数的判定

正确模型阶数的判定是个很重要的问题。因为用AR模型拟合数据序列时，若模型阶数选择低时，就如同用一低阶曲线去拟合一高阶曲线，会产生平滑的结果；而当模型阶数选择高时，则如同用一高阶曲线去拟合一低阶曲线，会产生急剧变化和振荡结果。换句话说，阶次太低将会导致平滑的谱估计值，而阶次太大又会引起伪峰并且会产生一般的统计不稳定性。

在低噪声或无噪声时，AR模型选择的阶数超过某个值还会产生谱线分裂现象。 在AR模型阶数判定中，最初是用白色检验法，即调整模型的阶数，使残差或预测误差最接近为白色的，就将此阶数作为AR模型最佳阶数。后来用准则函数法，即选取一个合适的准则函数，使其达到最小作为最佳阶数。

在准则函数方面，日本学者阿凯克（Akaike）的工作非常出色。1971年他首先提出了最终预测误差（FPE）准则，简称为最小FPE法（Final Prediction Error）准则。1973年，Akaike又提出了最小信息准则，称为阿凯克信息论准则，简称AIC准则（Akaike Information Criterion），近年来，并推广为BIC准则等。这些准则是目前应用比较广泛且比较有效的定阶方法。

八、AR谱估计存在的问题及其改善措施

在实际应用中常会观察到AR谱估计还存在一些问题，例如，虚假谱峰，谱线分裂，谱峰位置受相位影响，噪声使谱估计恶化，等等。为解决这些问题，人们进行了理论与实际研究，已经取得了不少进展，并继续寻找更加完善的理论解释与处理方法。

1、 虚假谱峰

?p(i)?0（i大于p）引起的，因虚假谱峰是由于自相关函数或反射系数的估计误差造成a?p(i)?0（i大于p）相应的将会产生n?p个额外的极点。若这些额外的极点出现在单为a位圆附近，就会形成虚假的谱峰。

2、 谱线分裂与谱峰频率偏移

如果要估计的随机过程是由一个正弦信号叠加噪声所构成的，那么在实验中会观察到：AR谱估计中谱峰出现的位置与正弦信号的初相位有着很密切的关系。

在谱估计中应该存在一个谱线的地方出现了两个紧挨着的谱峰，这一现象称为谱线分裂。谱峰频率偏移就是估计的谱峰位置偏离了真实谱峰。

谱峰位置对相位的依赖性随数据记录长度的增加而减小。

3、 噪声对AR谱估计的影响

当待分析的数据序列附加了观测噪声时，将使谱峰加宽而且变得平滑，导致分辨率下降，同时估计的谱峰偏离了真实谱峰。人们观察到，对于白噪声中含有两个幅度相等的正弦信号所构成的过程，AR谱估计的分辨率随信噪比的下降而下降。在信噪比低的情况下，AR谱估计已经不再优于周期图。

分辨率下降的原因，是AR谱估计中所假设的全极点模型，在有观测噪声的情况下，已

2经不在合适了。当存在观测噪声w(n)时，由于?w的作用，AR谱变为ARMA谱，它不仅存在仍由A(z)确定的谱极点，还存在谱零点。当观测白噪声方差??u/?22w2w很大时，即当

?0时，谱零点即为原来过程的谱极点，此时谱零点与谱极点相互抵消，产生完

全平滑的谱。因此谱零点有抵消极点的作用，而使谱变得平滑，并且这种抵消或平滑作用，

222222与?u/?w的大小有关。?u/?w愈小，功率谱愈平滑；?u/?w愈大，愈接近原来的AR谱。?2w?0时?2u/?2w??，即为原来的AR谱，这是理所当然的。

改进观测噪声影响的措施

由于附加了观测噪声，使得估计变得困难，人们提出了一些准最佳估计方法，以减小观测噪声对AR谱估计造成的性能下降。这些抗噪声技术包括：

①使用一般ARMA估计器； ②将数据进行滤波以减小观测噪声；

③噪声补偿法——补偿AR参数或反射系数估计值中由于噪声引起的偏差； ④使用高阶AR模型。

第五章 白噪声中正弦信号频率的估计

估计白噪声中正弦信号的频率是一个在许多领域中都会出现的问题。对于用傅里叶方法可以很好分辨的单个或多个正弦信号，通常可以用周期图作为频率估计器，在这种情况下，周期图是最佳的。然而对于不可分辨的正弦信号，就不能用周期图。

近似极大似然技术把正弦过程看作是一个频率未知的确定性信号。而特征分析法则采用广义平稳随机过程模型，使频率作为自相关矩阵的未知参数。虽然后一种方法在很大程度上依赖于前面讲述过的谱估计器，但应当把它看作是频率估计器，而不是谱估计器。

一、皮萨论科谱分解法

对于白噪声中正弦组合过程的混合谱估计，皮萨论科（Pisarenko）谱分解法是将其作为

一特殊情况的ARMA过程，并采用特征分析技术来求解。这种方法可以较为理想地恢复正弦信号的频率和幅度信息。

1、 正弦过程的时间序列模型与退化AR过程

无论是实正弦过程，还是复正弦过程，都可化为一个退化的AR过程来研究。 由p个实正弦所组成的过程为

p x(n)??i?1Aisinn(?i??i) （5.1）

式中初相位?i为区间（??,?)上均匀分布的独立随机变量，它在一次实现中为常量。

p个实正弦所组成的过程（5.1）可用2p阶差分方程

2p x(n)???a(k)x(n?k) （5.10）

k?1来表示，它是退化的AR（2p）过程。其p个频率应该由特征多项式：

2p

?a(k)zk?o?k?0 （5.8）

来确定。其中，a(0)?1，系数a(1),a(2),?,a(2p)必须受式（5.8）根共轭成对出现与模等于1的限制，即zi?1。 由上式不难得出

a(k)?a(2p?k),0?k?p

p个复正弦组成的过程为

p x(n)??Aeii?1j(2?fin??i) （5.11）

式中符号意义与（5.1）相同，在这种情况，?x(n)?是一个退化的AR（p）过程，即有

p x(n)???a(k)x(n?k) （5.12）

k?1式（5.12）的特征多项式为

p

j2?fi?ak?0kz?k?0 （5.14）

其根zi?e(1?k?p)，模仍为1，但不共轭成对出现。

由上可见，p个实正弦过程为退化的AR（2p）过程，其参量间存在

a(k)?a(2p?k)(0?k?p)的关系，独立参量个数为p个；p个复正弦过程为退化的AR(p)过程，独立参量个数仍为p个。实正弦过程相应的退化AR过程的阶数比复正弦过

程时相应的阶数高出一倍，这是由于前者特征多项式的根是成共轭对出现的，知其一，另一即可由其共轭关系确定。因此，实、复两正弦过程不仅描述它们所需的独立参量个数相同，并且它们相应的退化AR过程所需的独立参量个数也相同。

2、 白噪声中正弦组合为一特殊ARMA过程

若白噪声中正弦组合的过程为

p y(n)?x(n)?w(n)??i?1Aisin2(?fin??i)?w(n) （5.15）

式中w(n)为白噪声。利用式（5.10），得

2p2p y(n)??a(i)y(n?i)?w(n)??a(i)w(n?i) (5.17)

i?1i?1上述为ARMA(2p,2p)模型的一种特殊结构，与一般ARMA(2p,2p)模型相比，它有下面三个不同之处：

1、AR部分与MA部分具有相同的参量，它们存在共同因子；

2p2、由于特征多项式

?a(i)zi?o?i?0 的根的模为1，故AR部分特征多项式不满足平

稳性条件，MA部分特征多项式也不满足可逆性条件；

3、AR部分的y(n)?x(n)?w(n),y(n)是含白噪声w(n)的观测值，而通常为信号

x(n)，不含白噪声。

3、 特征分析技术求解法

对上述特殊结构的ARMA过程，不能用通常的ARMA模型法求解。若给定观测值y(n)的自相关函数，可用特征分析技术来求解。

将式（5.17）写成矩阵形式，有

ya?wa （5.18） 用矢量y左乘式（5.18）两边，并取数学期望得

E[yy]a?E[yw]a （5.19） 得下列特征方程：

Ryya??wa （5.20） 其中?w为自相关矩阵Ryy的特征值，a是对应特征值?w的特征矢量。Ryy可表示为 Ryy?Rxx??wI （5.21） 式中Rxx是?x(n)?的自相关阵，可以证明，当M?2p?1时（M为Ryy的维数即Ryy为

222TTTT2，Ryy的特征值为 ： M?M维）

MAi422w

?i1??i2??i??w,2???M2Pi??w,2i?2p?1?i?1,2,?,p?? （5.22） ??2w其中Pi?Ai22为第i个信号的平均功率。因此特征方程（5.20）中?为Ryy的最小特征值，

2w2a为对应最小特征值?w的特征矢量。当Ryy的阶数超过(2p?1)?(2p?1)时，?就是其

多重最小特征值。

由式（5.15），可求得白噪声中p个正弦波组合过程的自相关函数为：

归纳起来，利用Pisarenko分解法确定阶数p（正弦波个数）和正弦波频率fi的步骤如下： 算法1 （pisarenko谐波恢复） （1）计算矩阵

??ryy(0)ryy(?1)???ryy(1)ry(0)????????ryy(m)ryy(m?1)?ry(?m)??ry(?m?1)????ryy(0)?? Ryy(m?2p)

的特征值分解。

（2）求矩阵Ryy的最小特征值。与最小特征值对应的特征向量各分量用

a(0),a(1),?,a(2p)表示。

（3）求多项式的根，有

a(0)?a(1)z???a(2p)z**2p?0

*（4）将上述根记为z1,z1,z2,z2,?,zp,zp，其中 zi?ej2?fi

由此定出fi（取正值）。因

p ryy(k)?其中Pi?即

122?Pi?1icos(2?fik)??w?(k)

2Ai 为正弦波的功率，因此，接收信号的功率谱由上式的傅里叶变换给出，

p Pyy(f)??i?1Pi2[?(f?fi)??(f?fi)]??2w （5.30）

二、奇异值分解原理

奇异值分解是一种矩阵分解方法。对方阵而言，它就是特征分析。下面从特征分解开始讨论。

（一） 特征分解

设A为N?N阶的埃尔米特矩阵(A?AH)，有N个特征值及N个相应的特征矢量，将它们分别记为?i与?i，i?1,2,?,N。为讨论方便起见，将?i按非递增顺序排列，即，有 ?1??2????N，然后用矩阵形式表示式（5.61）

AV?V? (5.65) 式中矩阵V?[v1,v2,?,vN]

??10?0???0?2??? ???????0???0?0?N??对于埃尔米特矩阵，其特征值?i,i?1,2,?,N，是实数，并且对应的特征矢量?i,i?1,2,?,N是互相正交的，即

?i?j??H?0,?1,i?ji?j (5.66)

即V是酉矩阵（VH?V?1）。矩阵A的特征分解式为：

A?V?VH (5.69)

或

N A?（二）奇异值分解

??vviii?1Hi （5.70）

若A为M?N复长方矩阵，则有下列奇异值分解式：

A?U?VH （5.71）

式中U、V分别为M?M阶和N?N阶酉矩阵，并且

??'????00?? (5.72) 0?其中

??10?0???0?2?0'? (5.73) ???????????0?0?r??方阵?'的维数r不大于M、N两数中的小者，而?i,i?1,2,?,r为一组正实数。

三、正弦组合模型自相关矩阵的特征分析

目前在谱估计中所采用的奇异值分解处理，都是基于自相关矩阵的奇异值分解处理。由于这种矩阵往往是方阵，所以只需考虑其特征分析。白噪声中正弦组合信号是谱估计领域中常用的信号模型，它在应用中对谱估计的分辨率和稳定性都有很高的要求，（稳定性通常是指谱估计的抗干扰能力，这种干扰一般来源于相关函数估计的不准确性。）而正弦信号和白噪声对其自相关矩阵的特征分解有完全不同的贡献。因此，可以通过对其自相关矩阵的特征分析来区别两者，消除噪声的影响，从而提高谱估计的稳定性。

设信号是由P个复正弦加白噪声组成，其表示式为：

px(n)??i?1Aiexpj(2?fin??i)?w(n),n?0,1,?,M?1 (6.94)

2w式中w(n)是均值为0且方差为?的高斯白噪声。正弦信号的参数

{A1,A2,?,Ap;?1,?2,?,?p;f1,f2,?,fp}依次为振幅，相位和频率。假设相位?在

[0,2?]内为均匀分布的独立随机变量。

p rxx(k)?rss(k)?rww(k)?式（5.94）写成矩阵形式有：

?Ai?12iexp(j2?fik)??w?(k) (5.95)若将

2X?EA?W (5.96)

式中 X?[x(0),x(1),?,x(M?1)]

A?[A1ej?1T,A2ej?2,?,Apej?p]

T E?[e1,e2,?,ep] ei?[1,ej2?fi(M?p) ,?,ej2?fi(M?1),ej2?fi?2]

TT W?[w(0),w(1),w(2),?,w(M?1)]

信号 x(n),n?1,2,?M?1的自相关矩阵(M?M)为：

Rxx?E[XXH]?E[(EA?W)(EA?W)]?E[(EA?W)(AEHHH?WH)]

?EE[AAH]EH?E[WW ?EPEH2H]

??wI?Rss?Rww (5.97)

式中P?E[AAH]为p个正弦分量复振幅的自相关矩阵，I为M?M阶单位矩阵。若p个正弦分量彼此不相关，则

P?diag[P1,P2,?,Pp] (5.98)

式中 Pi为第i个正弦信号的功率。且Pi?E[(Aie因P是对角阵，（5.97）式又可表示为

pj?i)?(Aiej?i)]?Ai,*2i?1,2,?p

Rxx?EPEH??I?2w?Peeiii?1Hi??wI?Rss?Rww (5.99)

2可以看出，Rxx为信号自相关矩阵Rss与噪声自相关矩阵Rww之和。M?M阶（M?p时）信号自相关矩阵Rss?EPE中。

若信号自相关矩阵Rss?EPEHH不是满秩的（秩为p），也就是说Rss必定是奇异的。可是，

2总自相关矩阵Rxx是满秩的，因为它包括?wI项，频率信息包含在信号矢量ei(i?1,2,?,p)的特征值和相应的特征矢量分别记作?i?和vi，

(i?1,2,?,M)，并设?i?按非递增顺序排列，显然，最后M?p个特征值为0，即

??0 (5.100) ?????????Mp?1p?2MRss????vviii?1Hi (5.101)

pRss?'???vviii?1Hi (5.102)

??,?p称为主特征值（p个大特征值）特征值 ?1?,?2，特征向量v1,?vp称为主特征矢量。

主特征矢量与信号向量e1,e2,?ep张成相同的信号子空间。 因为 vi是正交的特征向量，所以

MI??vvii?1Hi (5.103)

则自相关矩阵

ppiiHiMiiHiRxx??Peei?1p??I?2w???vvi?1M??2w?vvii?1Hi

??(?i?1'i??)viv2wHi???i?p?12wvivi (5.104)

H其中特征值

2???i???w, ?i??2???w,i?1,2,?,pi?p?1,p?2,?,M (5.106)

另外，因为Rss有M?p个零特征值，则 EHvi?0,则

eiHvk?0,k?p?1,?M (5.109) i?p?1,p?2,?,M (5.108)

由上面的分析可知，当信号由p个线性独立的复正弦和白噪声组成时，自相关矩阵Rxx的特征分解有以下几个特点。

21、Rxx的最小特征值为?w，有M?p个，即

?1??2????p??p?1??p?2????M??w

2???i???w, ?i??2???w,2i?1,2,?,pi?p?1,?M

特别有意义的是主特征值?1,?2,?,?p及其相应的主特征向量v1,v2,?,vp分别为信号特征值和信号特征向量。而次特征值?p?1,?p?2,?,?M及其相应的次特征向量

vp?1,vp?2,?,vm分别为噪声特征值和噪声特征向量。

2、与最小特征值对应的特征向量vp?1,vp?2,?,vM与信号向量e1,e2,?,ep正交，即

H eivk?0,k?p?1,?,M,i?1,?,p。

{vp?1,vp?2,?,vM}?{e1,e2,?,ep} (5.110) 这种正交特性可以用来确定p个正弦分量的频率。

3、从特征分解为线性子空间分解的意义来看，Rxx特征矢量vi(i?1,2,?M)所张成的线性空间，可分解为两个子空间的直和：

Vs?[v1,v2,?,vp] 主特征向量组成信号子空间 Vn?[vp?1,vp?2,?,vM] 次特征向量构成噪声子空间

span{Vs}?span{v1,v2.?,vp} (5.111)

称为信号子空间，主要含有信号的信息，而

span{Vn}?span{vp?1,?,vM} (5.112)

称作噪声子空间，只含有噪声信息，不包含信号信息。

四、噪声子空间频率估计

多重信号分类（MUSIC）算法

多重信号分类（MUltiple Signal Calssification）简称MUSIC，是由Schmidt提出的一种著名方法。

设观测信号（或接收到的信号）为：

p?Aexp(j2?fn??)?w(n)n?0,1,?,M?1

?? ??Aee?w(n) (5.138)

x(n)?pi?1iiijij2fini写成矩阵形式有

i?1 X?EA?W （5.139） 对模型(5.139)，假定以下条件：

①M?p且对应于不同fi值的向量ei是线性独立的； ②E[w(n)]?0,E[w(n)w*(n?k)]??2?(k)； ③矩阵P?E[AAH]是非奇异的且为正定的。

当以上三个条件成立时，观测向量X的协方差矩阵由：

?Rxx?E[XXH]?EPEH??I （5.140）

2给定，其中P?E[AAH]。

注意到Rxx是一对称阵，其特征值分解具有下列形式：

?Rxx?V?V2H （5.141）

2其中： ??diag[?1,?2,?,?情况下，EPEH2222M],对角元素?i??i叫做Rxx的特征值。

在条件①即M?p时，矩阵E显然是非奇异的，即rank(E)?p。从而在条件③成立

的秩也为p，因此Rxx的特征值必须满足下列关系：

?i???i??22(i?1,2,?,p) (i?p?1,p?2,?,M)

令?1,?2,?,?p分别对应的特征向量为v1,?vp，而?p?1??M对应的特征向量为

vp?1,vp?2,?,vM，定义V1?[v1,v2,?,vp]，V2?[vp?1,vp?2,?,vM]。V1和V2分别叫

做信号子空间和噪声子空间。于是特征矩阵V分为两个子矩阵：

V?[V1V2]

2下面研究噪声子空间V2和矩阵E之间的关系。一方面，由于?和V2分别是Rxx的特征值和对应的特征向量，故有特征方程：

2 RxxV2??V2 (5.142)

同时有

EV2?0 (5.146)

H不难看出，真实参数值{f1,?fp}是上式的唯一解。

当V2为估计值时，EHV2不为零。这时，可取EHV2的2?范数为最小值的f?i作第k个信号源频率的估值。因此，连续改变f值，利用EHV21eV2VHH222进行搜索，由此得到p个最小值

所对应的f就是p个信号源的频率。实际做法之一是构造如下函数： PMU(f)??1M (5.149)

2e?i?p?1H?ie?其中e为频率的向量函数。上式最大值所对应的f就是信号源频率的估计值。

MUSIC算法概括如下： ?; ①估计样本自相关矩阵Rxx?作特征值分解; ②对Rxx?的最小特征值的数目n，求出这n个最小特征值?,?,?,令 ③确定Rp?1MEExx ?2?1nE(?p?1??p?2????M)

并求出与之相对应的特征向量vp?1,?,vM，构造噪声矩阵V2?[vp?1,?,vM] ④利用式（5.149）计算PMU(f)，并用它的谱峰来估计正弦信号的频率。其中 f?(i?1)?f，网格?f可取0.001等。 ⑤计算信号协方差矩阵P。

本文来源：https://www.bwwdw.com/article/amdr.html