对功率谱估计常用方法的探讨及相关matlab实现

对功率谱估计常用方法的探讨及应用分析DSP课程设计对功率谱估计常用方法的探讨及应用分析邹亮[e-mail]对功率谱估计常用方法的探讨及应用分析进行傅里叶变换在频域中研究信号，是研究确定性信号最简单且有效的手段，但在现代信号分析中，对于常见的随机信号，不可能用清楚的数学关系式来描述，其傅里叶变换更不存在，转而可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度。功率谱估计是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计和现代功率谱估计。本文介绍了各种经典功率谱估计方法，不仅从理论上对各种方法的谱估计质量进行了分析比较，而且通过Matlab进行了仿真。在对经典谱估计进行讨论之后，还分析了现代谱估计即参数谱估计方法，通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。现代谱估计的内容极其丰富，设计的学科及应用的领域都相当广泛，至今每年都有大量的科研成果出来。在本文的最后利用现代谱估计的方法讨论了功率谱方法在噪声源信号识别中的应用。文章还给出了常见谱估计方法的比较，便于深刻理解各种方法的特点，从而在实际工作中做出合理的选择。1.功率谱方法的发展功率谱估计是随机信号处理的重要内容，其技术渊源很长，而且在过去的40余年中获得了飞速的发展。涉及到信号与系统、随机信号分析、概率统计、矩阵代数等一系列的基础学科，广泛应用于人民的日常生活及军事、工业、农业活动中，是一个具有强大生命力的研究领域。本文将简要回顾一下功率谱估计的发展历程，对常用的一些方法进行总结。功率谱的估计方法有很多，主要有经典谱估计和现代谱估计。经典谱估计又可以分成两种：一种是BT法，也叫间接法；另一种是直接法又称周期图法。现代谱估计的方法又大致可分为参数模型谱估计和非参数模型谱估计，前者有AR模型、MA模型、ARMA模型、PRONY模型等，后者有最小方差方法、多分量的MUSIC方法等。1.1功率谱研究的发展过程功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。下面对谱估计的发展过程做简要回顾:英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法国工程师傅立叶提出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理论基础。傅立叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为对功率谱估计常用方法的探讨及应用分析“周期图”（periodogram）。这是经典谱估计的最早提法，这种提法至今仍然被沿用，只不过现在是用快速傅立叶变换（FFT）来计算离散傅立叶变换（DFT），用DFT的幅度平方作为信号中功率的度量。1958年，R,Blackman和J.Tukey首先提出BT法，并命名为布莱克曼-杜基谱估计器（简称BT谱估计器）。这种方法是先按照有限个观测数据估计自相关函数，再对其求傅里叶变换得到功率谱。在1965年FFT未出现以前，BT法一直是最常用的方法。周期图较差的方差性能促使人们研究另外的分析方法。1927年，Yule提出用线性回归方程来模拟一个时间序列。Yule的工作实际上成了现代谱估计中最重要的方法——参数模型法谱估计的基础。Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程，可以说，Yule和Walker都是开拓自回归模型的先锋。1930年，著名控制理论专家Wiener在他的著作中首次精确定义了一个随机过程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即，“功率谱密度是随机过程二阶统计量自相关函数的傅立叶变换”，这就是Wiener—Khintchine定理。该定理把功率谱密度定义为频率的连续函数，而不再像以前定义为离散的谐波频率的函数。1949年，Tukey根据Wiener—Khintchine定理提出了对有限长数据进行谱估计的自相关法，即利用有限长数据估计自相关函数，再对该自相关函数球傅立叶变换，从而得到谱的估计。1958年，Blackman和Tukey在出版的有关经典谱估计的专著中讨论了自相关谱估计法，所以自相关法又叫BT法。周期图法和自相关法都可用快速傅立叶变换算法来实现，且物理概念明确，因而仍是目前较常用的谱估计方法。1948年，Bartlett首次提出了用自回归模型系数计算功率谱。自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于1947年提出了解Yule-Walker的快速计算方法。这些工作为现代谱估计的发展打下了良好的理论基础。1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计的迅速发展。现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方差性能不好的问题。1967年，Burg提出的最大嫡谱估计，即是朝着高分辨率谱估计所作的最有意义的努力。虽然，Bartlett在1948年，Parzem于1957年都曾经建议用自回归模型做谱估计，但在Burg的论文发表之前，都没有引起注意。现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，至今，每年对功率谱估计常用方法的探讨及应用分析都有大量的论文出现。目前尚难对现代谱估计的方法作出准确的分类。从现代谱估计的方法上，大致可以分为参数模型谱估计和非参数模型谱估计，前者有AR模型、MA模型、ARMA模型、PRONY模型等；后者有最小方差方法、多分量的MUSIC方法等。非参数模型谱估计的特点是其模型不是用有限参数来描述，而直接由相关函数序列得到，这种方法能提高低信噪比时的谱分辨率。参数模型谱估计是先根据过程的先验信息或者一些假定，建立一个数学模型来表示所给定采样数据的过程，或者选择一个较好的近似实际模型，而后利用采样数据序列或者自相关序列，估计该模型的参数，最后把参数代入到该模型对应的理论功率谱表达式，得到所需要的谱估计。目前大量的论文集中在模型参数的求解上，以求得到速度更快、更稳健、统计性能更好的算法。1.2功率谱估计方法提出在通信系统中，往往需要研究具有目中统计特性的随机信号。由于随机信号是一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。1.3功率谱估计应用及用途功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数)。例如，在最佳线性过滤问题中，要设计一个维纳滤波器就首先要求知道(或估计出)信号与噪声的功率谱密度(或自相关函数)。根据信号与噪声的功率谱(或)(mxx)才能设计出能够尽量不失真的重现信号，而把噪声最大限度抑制的维纳滤波器。常常利用功率谱估计来得到线性系统的参数估计。例如，当我们要了解某一系统的幅频特性)(H时，可用一白色噪声n通过该系统。再从该系统的输出样本y(n)估计功率谱密度)(yyP。由于白色噪声的PSD(用)(P表示)为一常数即2)(P，于是有：22)()(HPyy(1-1)故通过估计输出信号的PSD，可以估计出系统的频率特性)(H(模特性)。从宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这对功率谱估计常用方法的探讨及应用分析要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频率间隙(距)。提高谱估计的分辨率已成为目前谱估计研究中的一个重要方向功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。维纳滤波、卡尔曼滤波，可用于自适应滤波，信号波形预测等（火控系统中的飞机航迹预判）。2经典谱估计2.1周期图法周期图法又称直接法。它是从随机信号x(n)中截取N长的一段，把它视为能量有限x(n)真实功率谱)(jwxeS的估计)(jwxeS的抽样.周期图这一概念早在1899年就提出了，但由于点数Ｎ一般比较大，该方法的计算量过大而在当时无法使用。只是1965年FFT出现后，此法才变成谱估计的一个常用方法。周期图法包含了下列二条假设：1．认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(nxN来估计该随机序列的功率谱。这当然必然带来误差。2．由于对)(nxN采用DFT，就默认)(nxN在时域是周期的，以及)(kxN在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(nxN的周期延拓，这也就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数)(mRx时将)(nxN以外是数据全都看成零，因此相关法认为除)(nxN外x(n)是全零序列，这种处理方法显然与周期图法不一样。但是，当相关法被引入基于FFT的快速相关后，相关法和周期图法开始融合。比我们发现：如果相关法中M=N，不加延迟窗，那么就和充（N-1）个零的周期图法一样了。简单地可以这样说：周期图法是M=N时相关法的特例。因此相关法和周期图法可结合使用。2.2相关法谱估计（BT）这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。它是1958年由Blackman和Tukey提出。这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N的有限长序列列)(nxN对功率谱估计常用方法的探讨及应用分析第二步：由N长序列)(nxN求（2M-1）点的自相关函数)(mRx序列。即)()(1)(10mnxnxNmRNnNNx(2-1)这里，m=-(M-1)…,-1,0,1…,M-1，MN，)(mRx是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。第三步：由相关函数的傅式变换求功率谱。即jwmMMmXjwxemReS)()(1)1((2-2)以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是将想x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jwxeS代表估值。一般取MN，因为只有当M较小时，序列傅式变换的点数才较小，功率谱的计算量才不至于大到难以实现，而且谱估计质量也较好。因此，在FFT问世之前，相关法是最常用的谱估计方法。当FFT问世后，情况有所变化。因为截断后的)(nxN可视作能量信号，由相关卷积定理可得)(mRx)]()([1mxmxNNN(2-3)这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对上式两边取（2N-1）点DFT，则有2121212)(1)()([1)(KXNKxkxNmRNNNx(2-4)于是将时域卷积变为频域乘积，用快速相关求)(mRx的完整方案如下：1.对N长)(nxN的充（N-1）个零，成为（2N-1）长的。2.求（2N-1）点的FFT，得220121212)()(NNmkNNNWnxKX。3.求212)(1KXNN。由DFT性质，)(12nxN是纯实的，)(12kxN满足共轭偶对称，而212)(1KXNN一定是实偶的，且以（2N-1）为周期。4.求（2N-1）点的IFFT：mkNNNkNxWKXNNmR121)1(212)(1121)((2-5)对功率谱估计常