功率谱估计

功率谱估计

引言：

对信号和系统进行的分析研究、处理有两类方法：一类是在时域内进行，维纳滤波、卡尔曼滤波以及自适应滤波等都属于时域处理方法；另一类方法是频域研究方法。对于确定性信号，傅里叶变换是在频率分析研究的理论基础，但是在实际生活中大多数信号是随机信号，而随机信号的傅里叶变换是不存在的，在实际应用中，通常通过采集和观测平稳随机过程的一个抽样序列的一段（有限个）数据，根据这有限个已知的数据来估计随机过程的功率谱问题来对随机信号进行分析，这即是频率谱估计。功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内通过用某种有效的方法来估计出其功率谱密度，从而得出信号、噪声及干扰的一些性质来，提取被淹没在噪声中的有用信号。功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。??谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。????

按照Weiner—Khintchine定理，随机信号的功率谱和其自相关函数服从傅里叶变换关系，可以得出功率谱的一个定义，如公式（1）所示：

Pxxe?jw???rm????xx(m)e?jwm 公式（1）

对于平稳随机信号，服从各态历经性，集合平均可以用时间平均来代替，可

以推出功率谱的另一定义。如公式（2）所示：

N1jwPxx?e??limE[x(n)e?jwn] 公式（2） ?2N?1n??NN??2

频率谱估计主要分为经典谱估计和现代谱估计，经典谱估计是将数据工作区外的

未知数据假设为零,相当于数据加窗,主要方法有相关法和周期图法;现代谱估计是通过观测

数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱,主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的,应用最广的是AR参数模型。谱估计方法的

分类如下图所示:下面来介绍其中的几种方法：周期图法、BT法、Bartlett法、Welch法以及参数模型法。

??周期图法????相关图法（BT图法）经典谱估计（线性谱估计）???Bartlett法???频率谱估计方法?Welch法??参数模型法（AR型、MA型和ARMA型）???现代谱估计（非线性谱估计）?特征分解法??最大熵谱法???

频谱估计方法的分类

经典谱估计

1、周期图法

周期图法又称直接法。它是从随机信号x(n)中截取N长的一段，把它视为能

量有限x(n)真实功率谱

Sx(ejw)的估计

Sx(ejw)的抽样.周期图这一概念早在1899

年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。只是1965年FFT出现后，此法才变成谱估计的一个常用方法。周期图法包含了下列二条假设：

1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一x(n)段N来估计该随机序列的功率谱。这当然必然带来误差。

x(k)在时域是周期的，以及N在频

x(n)域是周期的。这种方法把随机序列样本x(n)看成是截得一段N的周期延拓，

2．由于对

xN(n)采用DFT，就默认

xN(n)这也就是周期图法这个名字的来历。

随着样本长度 N 的增长,估计值渐进无偏,但是随着样本长度 N 的增长,周期图法并没有变得更平滑.可见,基本的周期图估计方法估计方差仍然比较大,难以满足一致性估计的条件.周期图法得到的功率谱分辨率高，但是方差性能差，功率谱起伏剧烈，容易出现虚假谱峰；周期图法得出的估计谱方差特性不好：当数据长度N太大时，扑线的起伏加剧；N太小时谱的分辨率又不好。采取平均

x(n)的方法对周期图法改进，平均就是将截取的数据段N再分成L个小段，分别计算功率谱后取功率谱的平均，增加分段数可以进一步减低估计的协方差，然而若每段中的数据点数太少，就会使估计的频率分辨率下降很多。在样本信号序列总点数一定的条件下，可以采用使分段相互重叠的方法来增加分段数，从而保持每段信号点数不变，这样就在保证频率分辨率的前提下进一步降低估计协方差。

2、相关图法（BT）

相关图法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法是先估计自相关函数，再进行傅里叶变换得到功率谱。它是1958年由Blackman 和Tukey提出。这种方法的具体步骤是：

第一步：从无限长随机序列x(n)中截取长度N的有限长序列列xN(n)

第二步：由N长序列xN(n)求（2M-1）点的自相关函数Rx(m)序列。即

1Rx(m)?N???xn?0N?1N(n)xN(n?m)

这里，m=-(M-1)…,-1,0,1…,M-1，MN，Rx(m)是双边序列，但是由自相关

函数的偶对称性式，只要求出m=0,.....,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即Sx(e)??jwM?1?Xm??(M?1)?R(m)e?jwm

以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将想x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的Sx(ejw)代表估值。一般取M<

?1可得： Rx(m)?[xN(m)?xN(?m)]

N 这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对

?112上式两边取（2N-1）点DFT，则有 Rx(m)?[x2N?1(k)?x2N?1(K)?X2N?1(K)

NN 随着自相关函数的时延取值M增大，分辨率提升，方差变大。当M=N-1时,BT法与周期图法估计出的功率谱是一样的;当 M < N-1 时,BT 法的偏差大于周期图法,在窗函数满足一定条件时是渐进无偏估计;方差小于周期图的方差;功率谱比周期图法要平滑；分辨率比周期图法低,与窗函数的选择有关。BT法的缺点在于当 M →N 时, R(m)的方差很大,使谱估计质量下降;由R(m) 得到的S(w)不一定为正值,从而可能失去功率谱的物理意义。

相关法在求相关函数Rx(m)时将xN(n)以外是数据全都看成零，因此相关法认为除xN(n)外x(n)是全零序列，这种处理方法显然与周期图法不一样。但是，当相关法被引入基于FFT的快速相关后，相关法和周期图法开始融合。我们发现：如果相关法中M=N，不加延迟窗，那么就和充（N-1）个零的周期图法一样了。简单地可以这样说：周期图法是M=N时相关法的特例。因此相关法和周期图法可结合使用。

对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。改进的方法有：平均周期图法和窗函数法,下面介绍的Bartlett法和Welch法实际上就是对周期图法的改进。

???3、Bartlett法:

Bartlett平均周期图法的思想： 对一个随机变量进行观测、得到L组独立的记录数据，用每一组数据求其均值，然后将L个均值加起来求平均。这样得到的均值，其方差将是用一组数据得到的均值的方差的1/L。

假设随机信号x(n)的观测数据区间为：O?n?M?1，共进行L次独立观测，

得到L组记录数据，表示为xi(n)，第i段的周期图为

1 Ii(?)?MM?1n?02?j?ni?x(n)e

将得到的L个周期图进行平均，作为信号x(n)的功率谱估计，有

L1?(e)??I(?) PxxiLi?1jw为分析偏移，对上式求统计平均，得

1Ljw?EPxx(e)??E?Ii(?)??E?Ii(?)? Li?1??周期图的统计平均值为：

?(ejw)?1EPxx2???????Pxx(ej(???))WB(ej?)d?

21?sin(Nw/2)??WB(ejw)?FT[WB(m)]???N?sin(w/2)??

可见，平均周期图法仍是有偏估计，偏移和每一段数据的个数M有关；由于

M

求平均周期图的方差：进行L次独立观测，L个周期图相互独立，因此平均周

?(ejw)?1Var?I(?)?。平均周期图的方差是周期图方差的期图的方差为：VarPxxiL1/L。显然分辨率的降低换取了估计方差和估计的均方误差减少。 实际中很难得到独立的数据组，一般常用的方法是将长度为N的数据分成L段，每段有M个数据，N=LM，第i段数据表示为：

?? xi(n)?x(n?iM?M) 0?n?M?1, 1?i?L

如果各段数据互不相关，利用平均周期图法估计功率谱的公式、偏移和方差的计算公式与上面结果一样。但实际中数据是连续的，数据段之间不能认为是完全不相关的，估计方差的减少小于1/L，或者认为近似为原来的1/L。可以根据实际情况适当选择I和M，如果对分辨率要求不高，可以取L大一些。最好是将N取得大一些，分辨率个估计误差都可以满足要求。

平均周期图的方法特别适合于应用FFT算法。因此在FFT出现以后这个方法比下面将要讨论的利用窗口函数的处理法用得更多。而在FFT出现以前主要是用窗口函数处理法平滑周期图。

2.4 Welch法：

Welch法又叫加权交叠平均法，这种方法以加窗（加权）求取平滑，以分段重

叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免带有两者的缺点，因此归根到底是一种折中。Welch提出了对Bartlett法的修正使更适合于用FFT进行计算。他主要提出二方面的修正，其一是选择适当的窗函数?(n)，并在周期图计算前直接加进法，这样得到的每一段的周期图为

1(i)IM(?)?MU1U?M这里

?x(n)?(n)ein?0N?12?j?n

?Wn?0N?12(n)为归一化因子，而Bartlett法每段的周期图为

N?1n?02i?j?n1(i)IM(?)?M?x(n)e

这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段时，可使各段之间有重迭，这样将会使方差减小(当N与M一定时)。 Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。

经典谱估计几种方法的比较与分析：

( 1) 周期图法和BT法的特点是离散性大，曲线粗糙，方差较大，但是分辨率较高；

(2) Bartlett法和Welch法的收敛性较好，曲线平滑，方差较小，但是功率谱主瓣较宽，分辨率低，这是由于对随机序列加窗截断所引起的Gibbs效应造成的； (3) 与Bartlett法相比，Welch法的估计曲线比较粗糙，但是分辨率较好，原因是 Welch 法中对数据进行截断时加的是Hanning窗，而在Bartlett法中使用的是矩形窗，相对于矩形窗，Hanning窗的主瓣包含更多的能量，因而使功率谱的主瓣较窄，分辨率较高。

在经典谱估计中，无论是周期图法还是其改进的方法，都存在着频率分辨率低、 方差性能不好的问题，原因是谱估计时需要对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，这其实是假定了窗以外的数据或自相关函数全为零，这种假定是不符合实际的，正是由于这些不符合实际的假设造成了经典谱估计分辨率较差。另外，经典谱估计的功率谱定义中既无求均值运算又无求极限运算，因而使得谱估计的方差性能较差，当数据很短时，这个问题更为突出。由上述理论分析及仿真实验可知，Welch法采用加窗交叠求功率谱，

可以有效减小方差和偏差，一般情况下能接近一致估计的要求，因而得到广泛应用。同时还可以发现，对信号加不同的窗函数，谱估计的质量是不同的。

参数模型法：

参数模型是将信号的随机性和确定性结合起来，激励白噪声体现过程的随机性，确定性模型则反映过程的可预测性。因此参数模型能较好地概括随机信号的性质。参数模型法为随机信号建立参数模型是研究随机信号的一种基本方法，其含义是认为随机信号是由白噪声激励某一确定系统的响应。只要白噪声的功率和系统的参数已知，随机信号的研究就可以转化为对模型参数和性质的研究：

w(n) h (n ) x(n) W(z) X(z) 最简单的参数模型是线性模型，即输入和输出的关系用线性差分方程描述，

H(z) 此时模型的传递函数是有理分式形式。根据传递函数形式，模型分为三种：滑动平均模型（MA）、自回归模型（AR）、自回归移动平均模型（ARMA）。

MA模型

信号由当前激励w(n)以及若干次以往激励w(n?k)的线性组合产生

x(n)??bkw(n?k)k?0qk?0 q 称为系统阶数。系统函数只有零点，没有极点，称为全零点系统，记为MA( q )。MA模型是一个全零点模型，容易反映谱的谷值，具有深谷而无锋锐的谱。由于MA (q)模型是一个全零点模型，共有q个零点，因此，其功率谱不易体现出信号中的峰值，因此分辨率较低。

2信号的功率谱为： qqH(z)??bkz?k

Sx(ej?)??21??b(k)e?j?kk?1 MA模型的相关函数满足以下方程：

q?m ?2q2b(k)b(k?m)??b(k)b(k?m),m?0,1,?q?? Rx(m)??k?mk?0?0,m?q

?

??MA参数模型求解方法是：

1.由N点数据x(n)，n=0,1,…,N-1，建立一个p阶的AR模型p远大于q，可求出p阶AR系数 a^(k),k?0,1,?pp

2.利用

ap(k),k?0,1,....p^，建立线性预测方程等效出一个q阶的AR模型，再一

次利用AR系数的求解方法，得到 b(k),k?0,1,.....p；

q3.计算功率谱： jw2S(e)??1??b(k)e?j?k k?1

从谱估计的角度来看，MA模型谱估计等效于经典谱估计中的自相关法，谱估计的分辨率较低。

AR模型

信号由当前激励w(n)以及若干次以往激励w(n?k)的线性组合产生

x(n)?w(n)??akx(n?k)k?1pH(z)?11??akz?kpk?1

p称为系统阶数。系统函数只有极点，没有零点，称为全极点系统，记为AR( p)。AR模型是一个全极点模型，容易反映谱的峰值，具有锋锐而无深谷的谱。

AR模型的自相关函数满足如下Yule-Walker方程：

?p ????iRx(m?i),m?0?i?1Rx(m)??p

?? ?iRx(m?i)??2,m?0???i?1

取m=0,1.....p,可得如下的矩阵方程：

Rx(1)?Rx(p)??1???2??Rx(0) ??R(1)?????R(0)?R(p?1)01xxx??????? ?????????????????? Rx(0)???p??0??Rx(p)Rx(p?1)?????在实际计算中，已知长度为N的序列x(n)，可以估计其自相关函数 利用矩阵方程，直接求出参数

Rx(m)^，再

?1,?2,...?p 及?，于是可以求出功率谱的估计值。

2综上所述，基本的AR谱估计的计算方法为：

（1）根据待估计信号的观测数据，估计出该信号的自相关函数。 （2）对信号的AR模型，选择恰当的模型阶数p；

（3）利用Yule-Walker线性方程组求解该信号的AR模型参数

?1,?2,?3.....?p及

?2

(4）根据下式计算信号的功率谱： 1Sx(w)?p

1???ie?j?n i?1

AR模型阶数的选择有以下三种准则： 1.最终预测误差准则： FPE ( k ) ? N ? k ?

kN?k

(k)?Nln?k?2k2.AIC信息论准则：AIC

k113.CAT准则： CAT(k)?1??Ni?1?i?k

AR谱估计的频率分辨率较高。AR模型的求解方法有自相关法、Burg算法和改进的协方差方法。它们共同的特点是由信号观测数据求模型系数时，均使用信号预测误差最小的原则。对于长记录数据，这些方法的估计质量相似；但对于短记录数据，不同的方法之间存在差异。自相关计算方法最为简单，但谱的分辨率相对较差；Burg方法计算不太复杂，可以由信号观测数据直接计算AR模型的参数，谱估计质量较好。改进的协方差算法给出了较好的谱估计性能，但计算较繁，编程较困难；

在AR谱估计中有以下的不足:1.与信号的信噪比关系较大，信噪比低，则方差大，分辨率低。2.如果信号x (n)是含噪声的正弦信号，其谱峰易受相位的影响，并且可能出现“谱线分裂”的现象。建议：选择合适的算法，AR谱估计前去除直流和线性倾斜成分等。3.谱的质量受p的影响，p取值过小，则过于平滑，精度不够，p太大，则可能会产生虚假的谱峰。建议：模型的阶不能过高，最高不超过数据长度的一半。

ARMA模型

ARMA 模型是AR模型与MA模型的结合：

x(n)??bkw(n?k)??akx(n?k)k?0k?1qpH(z)??bzkk?opk?1q?k

系统函数既有极点，也有零点，称为极零点模型，记为ARMA(p ,q )。ARMA模型能同时反映谱中的谷、峰值。

ARMA（p, q）模型的Yule-Walker方程为：

q?m?p2

???a(k)Rx(m?k)???G(k)b(m?k),k?0,1,?qk?0 R(m)??k?1?px ??a(k)Rx(m?k),m?q?? ?k?1

由于式中格林函数G (k)是ARMA的模型系数a (k)和b (k)的函数，所以上式的

1??akz?k第一个方程是非线性方程。但是当m大于q时有：

Rx(q)Rx(q?1)?Rx(q?p?1)??a(1)??Rx(q?1)? ???Rx(q?2)?Rx(q)?Rx(q?p?2)??a(2)? Rx(q?1)?????????????? ?????? ?Rx(q)???a(p)???Rx(q?p?1)Rx(q?p?2)????Rx(q?p)??

这是一个线性方程组，共有p个方程，可用来首先计算AR部分的系数，一旦a(1),a(2),…,a (p)求出，将它们代入第一个方程，再设法求解MA部分的系数，这是一种分开求解两个部分系数的方法。

模型选择考虑的原则

模型选择考虑的原则是模型能够表示谱峰、谱谷能力。对于具有尖峰的谱，应该选用具有极点的模型，如AR和ARMA模型；对于具有平坦的谱峰和深谷的信号，可以选用MA模型；既有极点也有零点的谱应该选用ARMA模型。对于没有模型准确反映其谱特性的信号，可以选用高阶的AR模型近似表示。选择模型的另一个考虑是尽量减少模型参数。当然这和选择模型是否合适有关系，前面介绍的三种信号模型均有普遍应用价值，但模型选择合适，功率谱估计和实际的谱拟合得好，如果不合适，只有提高阶数才能得到较为近似的谱，这样需要估计的参数增多，同样也会降低谱估计的质量。因此在选择模型合适的基础上，应尽量减少模型的参数。

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