三、功率谱分析字体[大][中][小]周期信号的功率谱为其双边幅值频谱的平方|cn|2;非周期信号的功率谱为其幅值谱密度的平方|X(ω)|2=X(ω)X*(ω)。随机信号属于时域无限信号,其频率、幅值和相位为随机变量。因而,采用具有统计特性的功率谱估计进行谱分析(一)自功率谱密度及其估计各态历经随机信号的功率谱密度Sx(ω)与自相关函数Rx(τ)为傅里叶变换偶对,即为了方便,也可用在非负频率范围内(ω0)定义的单边功率谱密度Gx(ω)代替双边功率谱密度Sx(ω),两者之间的关系为自功率谱估计可分为线性估计法与非线性估计法。前者以快速变换为基础,应用较早,也称为经典谱分析法;后者是与时序模型结合的一种新方法,又称为现代谱分析方法。1.周期图各态历经随机信号的均方值ψx2为信号能量的时域描述。巴什瓦定理表明,信号能量的时域计算与频域计算相等,即由此定义自功率谱密度及其估计为:式中表12-45典型信号的自相关、频谱、概率密度(续)X(ω)为测试数据x(t)的傅里叶变换,X(k)为N个数据x(n)的离散傅里叶变换,由FFT直接求出。由于X(k)具有周期函数的性质,所以称由此获得的自功率谱估计为周期图。自相关估计x′(r)的快速傅里叶变换可作为自功率谱估计的另一计算公式以上两种估计都是自功率谱Sx(ω)的有偏估计,只是偏差大小不同。两种估计在时域对数据或对自相关估计进行截断,相当于加窗处理,致使谱估计成为真实功率谱(或称为真功率谱)与窗谱W(ω)的卷积,即Ŝx(ω)=Sx(ω)*W(ω)窗谱旁瓣的泄漏效应和卷积的作用使真功率谱的尖峰数值变化,邻近点的数值变大,造成谱估计的模糊与失真以上两种估计的方差较大;相距2π/N的各点估计值互不相关,故数据点数N越大,这些点的估计值的随机起伏越严重。为改善谱估计的估计质量,在增大数据点数的同时,采用平均化处理和窗处理方法减小谱估计的方差。2.修正周期图平均化和加窗处理可使改善的周期图估计质量提高,而且还可保留周期图便于应用FFT计算速度快的优点。(1)平均周期图把N点的长数据序列分成k段,每段数据点数为M=N/k;求得各段周期图Ŝxi(ω)后再用平均法求得平均周期图Ŝxav(ω)。平均处理使谱估计的方差减小为当N一定时,段数K大则各段数据点数M小,谱估计的偏差大、方差小、谱平滑但频率分辨率低;若K小、M大,则偏差小、方差大。(2)加窗谱估计选择适当的窗函数ω(r),对自相关估计序列x′(r)作加窗处理,然后求谱估计加窗谱估计平滑、方差小,故又称为平滑周期图。常用窗函数有矩形窗、三角窗、汉宁窗、海明窗和高斯窗等。图12-7表示了分别用矩形窗、三角窗和汉宁窗分析同一数据所得结果。图a是被分析信号的真谱;图b和c是用两种时宽的矩形窗分析的结果。可以看出,时域窗口窄,则频域分辨率低,两相邻谱线无法区分。图d和e表明三角窗和汉宁窗频率分辨率低,但旁瓣衰减快,谱的分布区域窄而边沿清晰。显然,应按信号的性质和处理要求适当选择窗函数。图12-7窗函数效果比较(3)平滑平均周期图平滑谱估计计算要求所选用的窗函数保证求出的谱估计非负,但有的窗函数不满足此要求。如果先将数据x(n)分段,求出分段数据的加窗谱估计,然后再将各分段谱估计作平均化处理,即可满足谱估计非负的要求,又可减小估计偏差和估计方差,使谱估计质量提高。按此法计算的谱估计称为平滑平均周期图。N个数据x(n)分成k段,k=N/M。按照快速傅里叶变换的要求,将段内数据补零使M=2γ(γ为正整数)。求取各段的加窗谱估计然后用平均法求平滑平均周期图由于Pxi(k)计算用的数据经过加窗修正,为使谱估计为真谱的渐近无偏估计,在计算平滑平均谱估计Ŝx(k)时必须作相应的反修正,上式中U就是反修正因子——归一化因子(二)互功率谱密度对于零均值各态历经随机信号x(n)、y(n)的时域和频域描述,与自相关估计、自功率谱估计类似。但是,由于互相关函数并非偶函数,所以互相关函数的估计量分别给作互功率谱密度简称互谱,与互相关函数为傅里叶变换关系。互谱估计为互相关估计的离散傅里叶变换,可通快速傅里叶变换FFT求得。当有限长数据x(n)、y(n)的数据点数为N时,则估计方法和估计质量与自功率谱估计类似。在互谱分析中,常引用相干函数γy(ω)分析两个各态历经随机信号x(t)、y(t)之间以及它们与系统特性的相互关系。若在某频率上,γxy2(ω)=0,则表示x(t)与y(t)在此频率上不相干。如果x(t)与y(t)是统计独立的,则在所有频率上都有γxy2(ω)=0。若在所有频率上都有γxy2(ω)=1,则x(t)与y(t)完全相干。当γxy2(ω)1时,则说明:1)测量中有不通过系统的外来干扰;2)联系x(t)和y(t)的系统不是线性系统;3)y(t)是x(t)和其他输入的综合输出。相干函数的估计为其中的自功率谱和互谱的估计都应当是经过频率平均或总体平均后的估计,而且还应保证x(t)、y(t)的测试数据是同时测量、具有相同分析参数(如记录长度T等),以免分析结果出现错误。(三)极大熵谱估计伯格(J.P.Burg)提出的极大熵谱分析方法(MaximumEntropyMethod,简称MEM)是一种新的功率谱估计方法。伯格在最大熵谱估计准则的提出和具体算法上有所创新。由此演变出来的算法有多种,被统称为“现代谱分析”。1.极大熵准则传统谱估计方法实际上都是把无限长序列加窗截断后,由有限长序列获得功率谱估计。不论是对原始数据加窗还是对自相关函数加窗处理,其目的都是减小谱估计的方差、提高频率分辨率。然而窗处理不可避免地产生频域的“泄漏”,使功率谱失真。尽管在窗函数和处理方法上进行了许多研究,使得以周期图为基础的谱估计方法广泛应用,但谱估计的频率分辨率并不能令人满意。此外,传统谱估计方法通过增大数据点数来获得较高的谱估计精度。这样,不仅数据处理工作量大,而且对短记录数据或缓变信号等显然是无能为力的。在谱估计计算中,对原始数据或自相关函数加窗处理,假设窗外的数据为零,而且还对窗内部分进行某些修正。这些人为措施增加了确定性因素,使原来具有的不确定性减少,在一定程度上歪曲了观测得到的信息。极大熵谱分析方法在谱估计计算中不作加窗处理,而是采用适当的方法把由有限长数据求得的自相关估计外推。外推的原则是使相应的数据序列在外推点上取值的可能性具有最大的不确定性,亦即不对结果人为地增加任何附加信息。在数理统计学中,“熵”表征了各种随机试验的统计特性,是随机总体的平均不确定性的量度。极大熵谱分析法把熵的概念引入谱分析。上述外推原则就是要求在使随机过程的熵达到最大的条件下,确定未知的自相关函数值,外推原则的最大熵描述就是谱估计的极大熵准则。2.极大熵谱估计的基本原理随机试验a有有限个不相容的结果Ai,相应的概率为P(Ai),且满足则随机试验的熵H(a)为其不确定性的量度。对于连续随机变量,则熵为式中的对数可以取10或取e为底,在比较熵的大小时并没有影响正态随机变量x的概率密度为则熵为高斯随机过程的样本函数x(n),其数据点数为N,则其熵为式中|R|为N×N阶自相关矩阵[R]的行列式,也可用det[R]表示。[R]为托布列兹矩阵。方差为σ的自噪声的熵为Hε=-E[ln(p(ε1)p(ε2)…p(εN))]=Nlnσε由于H随N增长将要发散,故用熵率h代替,熵率定义为由于功率谱密度与自相关函数之间为傅里叶变换关系,则自相关矩阵[R]的行列式|R|的极限与功率谱密度Sx(f)的关系是:式中fN=2fc,为Nyquist采样频率;fc——信号的频限,即带宽,由此可以得到熵率为:当由有限长随机信号测试数据计算自相关估计并外推时,为使每步估计的熵达到最大,则应有(k≥m+1,m——数据点数)由此得到:式中,功率谱密度Sx(f)应受到熵率导数为零的条件约束,1/Sx(f)可用截短的傅里叶级数表示成为2m+1项有限级数之和,于是有和从而求得极大熵谱式中△t——采样时间间隔;ak——自回归模型参数,ak*与ak共轭;Pm——外推误差的方差MEM谱与参数ak和Pm有关。典型的算法有Burg、Levinson、Akaike和Marple算法。3.极大熵谱的特点与经典谱估计方法相比,极大熵谱估计方法克服了窗处理带来的一系列缺陷。而且由于是连续谱,谱线光滑、谱峰陡峭,频率分辨力远高于经典谱估计;而其谱分辨力与样本函数长度无关,特别适用于短数据样本、缓变过程的谱估计。图12-8表示短数据分析时FFT谱与MEM谱的对比。被分析信号是1Hz正弦波上叠加10%的白噪声。图a样本长度为1s时,两种谱均显示出主峰频率为1Hz;图b表示样本长度为0.57s时,MEM谱主峰频率清晰准确,而FFT谱却无法识别。图12-9表示缓变信号的MEM谱与FFT谱的比较。被分析的信号是0.6Hz的正弦波,由图可见,MEM谱较为精确。图12-8短数据的MEM谱与FFT谱比较图12-9缓变数据的MEM谱与FFT谱谱分析与谱估计在生产实践和科学研究中获得了日益广泛的应用。例如,在声纳系统中,为了寻找海洋水面舰艇或潜艇,需要对噪声信号进行谱分析,以提取有用信息,判断舰艇速度、方位、大小等;对各种机电产品主体或部件作实际运行的谱分析,可检验设计的效果并提出改进设计的依据,或者进行在线监测与故障诊断,以便及时排除潜在的故障因素,保证安全运行。