现代功率谱估计.doc

淮北师范大学2012届毕业论文现代功率谱估计现代功率谱估计淮北师范大学物理与电子信息学院 235000摘要 功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。它是随机信号处理的重要内容，广泛应用于人民的日常生活及军事、工业、农业活动中。其实现方法主要可分为经典谱估计和现代谱估计。经典谱估计方法由于其种种缺点，迫使人们大力研究现代谱估计方法。现代谱估计法是以参数模型为基础的方法，大致可以分为参数模型谱估计和非参数模型谱估计，前者有AR模型、MA模型、ARMA模型、PRONY模型等；后者有最小方差方法、多分量的MUSIC方法等。本文将着眼于现代谱估计的各种方法，首先简要介绍随机信号功率谱估计的相关基础知识，然后从经典法入手，探讨现代谱估计的理论基础，分析各种方法的优劣性及适用范围，并且给出对应的Matlab仿真结果，从而深刻理解各种方法的特点，从而在实际工作中做出合理的选择。关键词 功率谱 估计 现代 信号处理 Matlab引言功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理论基础。 傅立叶级数提出后， 19世纪末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为“周期图”（periodogram）。这是经典谱估计的最早提法，这种提法至今仍然被沿用。周期图较差的方差性能促使人们研究另外的分析方法。1927年，Yule提出用线性回归方程来模拟一个时间序列。Yule的工作实际上成了现代谱估计中最重要的方法参数模型法谱估计的基础。Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程，可以说，Yule和Walker都是开拓自回归模型的先锋。 1948年，Bartlett首次提出了用自回归模型系数计算功率谱。自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker的快速计算方法。这些工作为现代谱估计的发展打下了良好的理论基础。1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计的迅速发展。现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方差性能不好的问题。1967 年，Burg 提出的最大熵谱估计，即是朝着高分辨率谱估计所作的最有意义的努力。由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。本文将着眼于现代谱估计的各种方法，首先简要介绍随机信号功率谱估计的相关基础知识，然后从经典法入手，探讨现代谱估计的理论基础，分析各种方法的优劣性及适用范围，并且给出对应的Matlab仿真结果，从而深刻理解各种方法的特点，从而在实际工作中做出合理的选择。1 从经典谱估计到现代谱估计首先，给出功率谱的两个最基本的定义如下：Sxxej=k=-Rxxke-jk (1)Pxej=limME12M+1n=-MMxne-jn2 (2)可以证明，这两个定义是等效的。无论是建立在第一个还是第二个公式上的定义，在实际中都几乎是不可能实现的（除非x(n)可以用解析法精确的表示），因此，只能用所得的有限次记录（往往仅一次）的有限长数据来予以估计，这就产生了功率谱估计这一极其活跃，同时也极其重要的研究领域。1 经典谱估计之周期图法在实际应用中，通常观测到的是信号的有限个(例如N个)取样值，用yN(n)表示。可以认为它是分段平稳随机信号的一段，也可将它看成是从平稳随机信号中截取的一段数据。对于平稳随机信号，无论从何时开始任取一段长为N的数据，所计算出来的均值或自相关函数都是相同的。信号yN(n)可以看成是用一个宽为N的数据窗w(n)从平稳随机信号yN(n)中截取出来的，即yNn=ynwn (3)根据遍历性，用时间平均代替集合平均，若已知N个数据为yNn=y0,y1,yN-1 (4)则用时间平均来近似计算的自相关函数为Ryyk=1Nn=0N-1-|k|yn+kyk , kN-1 (5)称之为取样自相关。它可以看成是有限长序列yNn和yN-n的卷积运算结果除以N，即Ryyk=1NyNn*yN-n (6) 取样自相关函数的双边Z变换叫做周期图，它是功率谱的一种估计，用 Syyz 表示，Syyz=k=-(N-1)N-1Ryykz-k (7)联系式(9)，由上式得到Syyz=1NYzYz-1 (11)这里Yz是yNn的Z变换。式(10)和(11)是计算周期图的两种基本方法，前者称为间接法，后者成为直接法。令z=ej,由式(11)得到Syyz=1NY()2=1Nn=0N-1yN(n)e-jn2 (12)该式很适合用FFT计算。 改进周期图的有4种办法：修正周期图法，平均周期图法，加床平滑法，Welch法。2 经典谱估计之自相关法根据Wiener-Khintchine定理，平稳离散随机信号x(n)的自相关函数Rxxm=Ex*nxm+n (3)与功率谱Sxx之间构成一对傅里叶变换关系，即Sxx=m=-Rxx(m)e-jm (4)Rxxm=12-Sxx()ejmd (5)这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。它是1958年由Blackman 和Tukey提出。这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N的有限长序列列第二步：由N长序列求（2M-1）点的自相关函数序列。即 (6)这里，m=-(M-1),-1,0,1,M-1，MN，是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,M-1的傅里叶变换，另一半也就知道了。第三步：由相关函数的傅式变换求功率谱。即 (7)以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的代表估值。一般取M0 自相关函数的头p+1个值是R0,R1,Rp，因此，式_表示成下列矩阵形式：R0R1R2 RpR1R0R1 Rp-1R2 R1 R0 Rp-2 RpRp-1Rp-2 R01a1a2ap = 2000这就是AR(p)模型的Yule-Walker方程。只要已知或估计出p+1个自相关函数值，即可由该方程解出p+1个模型参数a1,a2,2。1 使用Levinson-Durbin算法对该方程求解。这是一种递推方法，当整个计算结束后，不仅求得了所需要的p阶AR模型的参数，而且还得到了所有各低阶模型的参数。它的关键是推导出由AR(k)模型的参数计算AR(k+1)模型的参数的迭代计算公式，这些公式如下：k2=R0+i=1kak,iRi Dk=i=1kak,iRk+1-i, ak,0=1k+1=Dkk2k+12=(1-k+12) k2ak+1,i= ak,i-k+1ak,k+1-i , i=1,2,k将模型参数带入，即可计算功率谱估计值：Sxxej=p21+i=1pap,ie-ji2递推公式提供了一种确定模型阶数的实验方法，如模型的阶数不知道，由低阶开始递推，当递推到M阶时，预测误差满足允许的值，停止递推，选AR模型的阶数为M。2 Burg算法Burg算法是Burg于1975年提出的求解AR参数的有效方法，其特点是在Levinson算法的基础上，不对自相关函数进行估计，而是利用前、后向线性预测系数之间的递推关系，直接求出反射系数。Burg算法与自相关法不同，它是使序列x(n)的前后向预测误差功率之和：PPfb=1N-Pn=pN-1ePf(n)2+ePb(n)2最小。当阶次由1至p时，ePf(n)和ePb(n)有以下的递推关系:emfn=em-1fn+kmem-1bn-1embn=em-1bn-1+kmem-1fnm=1,2,p初始条件为：e0fn=e0bn=x(n)将大括号带入PPfb式，可知PPfb仅为km的函数。令PPfbkm=0，可得km=-2n=mN-1em-1fnem-1bn-1n=mN-1em-1fn2+em-1bn-12 , m=1,2,p再利用Levinson-Durbin递推算法可得AR模型系数：amm=kmamk=am-1k+kmam-1m-k , k=1,2,pPmfb=Pm-1fb1-km2Burg算法是建立在数据基础之上的，避免了先计算自相关函数从而提高计算速度，是较为通用的方法。3 协方差法及其改进方法不作观测数据区间之外的数据为0的假设，在均方误差意义下使得数据的前向预测误差最小。由此估计的自相关矩阵式半正定的，且不具有Toeplitz性，得到的AR模型可能不稳定。同Burg算法一样，改进协方差算法进行功率谱估计时令前后向预测误差功率之和最小，即对前后向预测误差都不加窗，但得到的协方差矩阵不是Toeplitz矩阵，因此正则方程不能用Levinson方法求解。Maprle于1980年提出了一种快速求解方法，具体请参考文献*。3 MA模型如果除a0=1外所有其他的AR系数都等于零，这种模型称为q阶滑动平均模型，简称为MA(q)模型。其传输函数为HMAz=Bz=k=0qbkz-k模型输出功率谱为Sxxz=2BzBz-1或Sxxej=2Bej2=2k=0qbke-jk2这是一个全零点模型。3现代谱估计方法的特点根据Wold分解定理，可得推论：如果功率谱连续，那么任何ARMA过程或AR过程可以用一个无限阶的MA过程表示。Kolmogorov提出的一个定理有着类似的结论：任何ARMA或MA过程可以用一个无限阶的AR过程表示。估计ARMA或MA模型参数一般需要解一组非线性方程，而估计AR模型参数通常只需解一组非线性方程，因此AR模型得到了深入的研究和广泛的应用。如果被估计过程是ARMA或MA过程，或者是高于p阶的AR过程，那么用AR作为它们的模型估计时，虽然不可能很精确，但可以尽可能的逼近之，关键是选择足够高的阶。MA谱估计以全零点模型为基础，将其用于估计窄带谱时得不到高分辨率，但用于MA随机过程时，由于MA随机过程的功率谱本身具有宽峰窄谱的特点，故能得到精确估计。当采用AR谱估计，特别是采用Burg法时，能得到可靠的高分辨率估计。但当噪声污染了数据时，只有采用ARMA模型才能获得良好的谱估计。采用ARMA模型，以较少的模型参数就能改善AR谱估计的性能，因而关于ARMA模型方法的研究也受到了人们的普遍重视。常用的几种AR模型参数提取方法有：1)自相关法假定观测数据区间之外的数据为0，在均方误差意义下使得数据的前向预测误差最小。由此估计的自相关矩阵式正定的，且具有Toeplitz性，可以用Levison-Durbin算法求解。自相关法的计算效率高，且能保证预测误差滤波器是最小相位的，但数据两端要附加零取样值，实际上等效于数据加窗，这将使参数估计的精度下降。特别是当数据段很短时，加窗效应就更为严重。2)Burg法Burg法一方面希望利用已知数据段两端以外的未知数据，另一方面又总是设法保证使预测误差滤波器是最小相位的。它不直接估计AR参数，而是先估计反射系数，然后利用Levinson递推算法由反射系数来求得AR参数。在约束AR模型的参数满足Levison-Durbin递归条件的前提下，在均方误差意义下使得数据的前向预测误