功率谱估计方法与实现的研究_崔桂华(1)

姓名：崔桂华    学号：SC10023028   Email：cuigh@mail.ustc.edu.cn

一、背景介绍

1.1科学意义

在许多科学领域和生活实际中，由于各种干扰因素的存在，我们需要处理随机信号。常见的随机信号有各种无线电系统及电子装置中的噪声与干扰，建筑物所承受的风载，船舶航行时所受到的波浪冲击，包括心电图、脑电图和肌电图等在内的许多生物医学信号等等。随机信号和确定信号不同，它不能通过一个确定的数学公式来描述，也不能准确地预测。因此，对随机信号一般只能在统计意义上研究。

一个随机信号包含了无穷多个样本，每一个样本都是一个无穷长的时间序列，因此，要得到一个完整的随机信号是很困难的。即使可以得到，由于在计算机处理信号时，信号总是有限长，因此存在信号的截短问题。另外，信号从产生到传送再到采集，总要受到噪声干扰。对随机信号处理的一个重要任务就是由有限长且受到干扰的信号中得到信号的某些特征，如信号的均值、方差、自相关函数及功率谱等，或恢复出没有被干扰的信号。基于随机信号的上述特点，信号特征的提取或信号的恢复都要通过估计来实现，必须涉及到估计的理论和方法。

功率谱估计在随机信号处理中占有很重要的地位，这一技术广泛应用于雷达、声纳、通信、地质勘探、天文和生物医学工程等众多领域。例如在雷达系统，如何从被目标反射回来时受到严重干扰的微弱的回波信号中，提取发送的脉冲信号并通过计算回波到达时间和频率偏移进而确定目标的方位和运动速度，这就是一个典型的估计问题。又如在控制系统，如何从接收到的被噪声干扰的信号中分离出有用的控制信号。如果噪声是加性的且与有用信号各占一个频带，则可方便地采用低通、高通或带通滤波器，将干扰滤除，把有用信号提取出来。如果噪声与有用信号的频谱互相交叠在一起，则人们就很难利用上述频率选择滤波器把有用信号提取出来。为此，必须建立在随机过程理论基础上，从统计观点出发，通过对有用信号和噪声统计特性的分析，采用更复杂的滤波方法提取特定信息。这种过滤随机信号过程，实质上也是一种估计，所以估计器也是一种滤波器，它是统计信号处理的重要内容，其理论和技术在工程中得到广泛应用。如在地质勘探中，人们通过量测人工地震信号来预测地层结构，从而估计石油贮存的位置与数量；在化学工程中，通过温度与压力的测量来确定特殊液体的密度等等[1]。 

1.2 国内外相关研究现状分析

“谱”的概念最早由英国科学家牛顿提出[2]。后来，1822年，法国工程师傅里叶提出了著名的傅里叶谐波分析理论。该理论至今仍然是我们进行信号分析和处理的理论基础。傅里叶级数首先在观察自然界中的周期现象时得到应用，如在研究声音、天气、太阳黑子的活动、潮汐等现象时用于测定其发生的周期。由于傅里叶系数的计算是一个困难的工作，所以又促使人们研制相应的机器，如英国物理学家Thomson发明了第一个谐波分析仪用来计算傅里叶系数，这些机器也可用新得到的傅里叶系数预测时间波形。利用该机器画出某一港湾一年的潮汐曲线约4小时，这些都是人们最早从事谱分析的尝试。

1.2.1 经典谱估计方法

19世纪末，Schuster提出用傅里叶系数的幅平方作为函数功率的测量，并命名为周期图[3]，这是经典功率谱估计的最早提法，至今仍然被沿用，只不过我们现在是通过FFT来计算离散傅里叶变换。由于周期图起伏剧烈，Schuster提出了平均周期图的概念，并指出了在对有限长数据计算傅里叶系数时所存在的“边瓣”问题，这就是我们后来所熟知的窗函数的影响。Schuster用周期图来计算太阳黑子活动的周期，以1749年至14年每月太阳的黑子数为基本数据，得出黑子的活动周期是11.125年，而在天文学文献中太阳黑子的活动周期为11年。

1930年，著名控制论专家Wiener出版了他的经典著作Generalized Harmonic Analysis[4]。他在该书中首次定义了一个随机过程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即功率谱密度是随机过程二阶统计量自相关函数的傅里叶变换。这就是Wiener-Khintchine定理。该定理将功率谱密度定义为频率的连续函数，而不再是以前离散的谐波频率的函数。1949年，Tukey根据Wiener-Khintchine定理提出了对有限长数据做谱估计的自相关法，即利用有限长数据估计自相关函数，再用该自相关函数做傅里叶变换，从而得到估计谱[5]。Blackman和Tukey在1958年出版的有关经典谱估计的专著中讨论了自相关谱估计法[6]，因此后人又把经典谱估计的自相关法成为BT法。

周期图法和自相关法是经典谱估计的两个基本方法。它们都存在着以下两个难以克服的固有缺点：

(1) 频率分辨率(区分两个相邻频率分量的能力)不高。这是因为它们的频率分辨率(以赫兹计)反比于数据记录长度(以秒计)，而实际应用中一般不可能获得很长的数据记录。

(2) 经典谱估计方法在工程中都是以离散傅立叶变换(DFT)为基础的，它隐含着对无限长数据序列进行加窗处理(加了一个有限宽的矩形窗)。矩形窗的频谱主瓣不是无限窄的，且有旁瓣存在，这将导致能量向旁瓣中“泄漏”，主瓣变得模糊不清。严重时，会使主瓣产生很大失真，甚至主瓣中的弱分量被旁瓣中的强泄漏所掩盖。

为了克服以上缺点，人们曾做过长期努力，提出了平均、加窗平滑等办法，在一定程度上改善了经典谱估计的性能。实践证明，对于长数据记录来说，以傅立叶变换为基础的经典谱估计方法，的确是比较实用的。但是经典方法始终无法根本解决频率分辨率和谱估计稳定性之间的矛盾，特别是在数据记录很短的情况下，这一矛盾尤为突出。

1.2.2 现代谱估计方法

针对经典谱估计方法的缺点，现代谱估计方法提出并得到了迅速发展。目前，现代谱估计的主要方法分为参数模型方法和非参数模型方法，前者有AR模型、MA模型、ARMA模型、PRONY指数模型等；后者有最小方差方法、多分量的MUSIC方法等。根据Wold分解定理可以推导出无穷阶的MA模型可用一个有限阶的AR模型来等效，也就是说一个有限阶的AR模型可近似为一个高阶的MA模型；一个有限阶的ARMA模型可用一个阶次足够高的AR模型或MA模型来近似[7]。

AR模型参数的求解方法有自相关法、Burg法、协方差方法、改进的协方差方法和最大似然估计等方法。

1.2.2.1 自相关法[8] 

AR模型可表示成

将其代入自相关函数的表示式，得

由于u(n)是方差为σ2的白噪声，有

由Z变换的定义，当z→∞时，有h(0)=1，有

上式可写成矩阵形式

上述两式即是AR模型的Yule-walker方程。系数矩阵不但是对称的，而且沿着和主对角线平行的任一条对角线上的元素都相等，这样的矩阵称为Toeplitz矩阵。用线性方程组的常用算法(例如高斯消元法)求解上式，需要的运算量数量级为O(p3)。但若利用系数矩阵的对称性和Toeplitz性质，则可构成一些高效算法，Levinson-Durbin算法是其中最著名、应用最广泛的一种。这种算法的运算量数量级为O(p2)。这是一种按阶次进行递推的算法，即首先以AR(0)和AR(1)模型参数作为初始条件，计算AR(2)模型参数；然后根据这些参数计算AR(3)模型的参数，等等，一直到计算出AR(p)模型参数为止。这样，当整个迭代计算结束后，不仅求得了所需要的p阶AR模型的参数，而且还得到了所有各低阶模型的参数。定义am(k)为p阶AR模型在阶次为小时的第k个系数，k=1，2，…，m，m=1，2，…，p，以为m阶时的前向预测的最小误差功率。

当m=1时，有

解出

a1(1)=-rx(1)/ rx(0)

ρ1= rx(0)[1- a21(1)]

定义初始条件

ρ0= rx(0)

由Toeplitz矩阵的性质，可得到如下Levinson-Durbin递推算法：

Levinson-Durbin算法从低阶开始递推，直到阶次p，给出了在每一个阶次时的所有参数，这一特点特别有利于选择AR模型的合适的阶次。

由于线性预测的最小均方误差总是大于零的，必有|km|＜1,如果|km|=1，那么递推应该停止。由反射系数的这一特点，我们可以得出预测误差功率的重要性质

ρp<ρp-1<…<ρ1<ρ0

实际中，我们往往并不能精确地知道x(n)的前p+1个自相关函数，而知道的仅仅是N点数据，即xN(n),n=0,2,…,N-1，为此，可以：

(1)首先由xN(n)估计x(n)的自相关函数，得r’x(m)；

(2)用r’x(m)代替上述递推算法中的rx(m)，重新求解Yule-Walker方程，这时求出的AR模型参数是真实参数的估计值，即a(1),a(2),…,a(p)：

(3)将这些参数代入功率谱估计公式，得到x(n)功率谱的估计。

1.2.2.2 修正协方差算法

确定AR模型参数的Burg方法基本上是具有预测系数满足Levinson递推的附加约束的最小二乘格型算法。由于这一约束的结果，增加AR模型的阶仅需要在各级最优化一个参数。与这种方法相比我们可以利用无约束改进协方差法(最小二乘算法)确定AR参数。

该算法的主要思路是为了摆脱因采用递推运算对确定预测系数的约束，一旦求出ρp，则ap(p)= ρp，而其余的预测系数就由低一阶模型的系数ap-1(k)来确定，不能灵活地选取)，让每一预测系数(模型参数)的确定直接与前、后向预测的总的平方误差最小(最小二乘法)联系起来。即总的平方误差：

它与Burg方法有相同的性能指标。但是，我们并不为了求AR参数而把Levinson递推加入上式。可得到一组线性方程：

其中，自相关rx(l，k)定义为

所得到的剩余最小二乘误差为

所以，可得改进协方差功率谱估计为

元素为rx(l，k)的相关矩阵不是Toeplitz阵，因此，不能利用Levinson算法求解。但是，相关矩阵的结构足以设计复杂性与p2成正比的有效计算方法。Marple[9] [10]设计了这样的算法，他的算法具有格型结构，利用了Levinson类型的阶递推和另外的时间递推。这种形式的修正协方差(无约束最小二乘法)方法也一直被称为数据未加窗最小二乘方法。修正协方差方法基本上克服了谱线、谱峰偏移和出现伪峰等缺点，在此意义上优于Burg方法。但此法不能保证AR模型稳定，并且所需的运算量也偏大。

1.2.2.3 AR模型阶次的选择

模型定阶问题是参数谱估计的一个难题。阶次选得太低，功率谱受到的平滑太厉害，平滑后的谱可能会分辨不出真实谱中的两个相邻很近谱峰。阶选得太高，固然会提高谱估计的分辨率，但同时会产生虚假的谱峰或谱的细节，都不能反映模型的真实情况。既要估计AR模型参数，又要估计模型的阶，对于这样复杂的情况，如何评价各种谱估计的性能，目前尚未见到有人详细研究。己经提出几种不同的误差准则作为确定模型阶的依据。以下是其中主要的三种：

(1)最终预测误差(FPE)准则[11] [12] 

AR(k)过程的最终预测误差定为

它是过程中不可预测部分的功率与AR参数估计不精确产生的误差功率之和。上式中N是数据样点数目。括号内的数值随着k的增大(趋近于N)而增加，反映出预测误差功率的估计的不精确性增加。由于盯：随阶的增加而减小，所以FPE将有一个最小值。FPE的最小值对应的阶便是最后确定的阶。该准则的实际应用表明，虽然对于AR过程来说效果很好，但在处理地球物理数据时，一般都认为这一准则确定的阶偏低。

(2)Akaike信息准则(AIC)法

对AR过程，AIC定义为

这里i是假设的AR模型的阶，显然，对于上式，AIC有一个最小值，它所对应的阶就是要选择的阶。可以证明，当N→∞时FPE和AIC等效的，即

(3)判别自回归传输函数(CAT)法[13] 

这一准则是：使实际预测误差滤波器和估计滤波器的均方误差之间的差的估计的最小所对应阶作为最佳阶。

使得上式最小的p值就是最终确定的阶。

1.2.2.4 Prony复极点模型[14] 

基于线性预测理论，用一个指数衰减线性组合模型来拟合一个等间隔的离散

信号其表达式为：

其中xn为时间序列取样值，Si=ai+jωi为复极点，ai为衰减因子，ωi为自然频率，Ai为留数，Δt取样间隔，M为取样总数。上式为M个非线性方程式，有2N个(N个Ai，N个Si)未知数。当M≥N时，Mittra与BlaricuIII引用古老的Prony算法求解了这个非线性方程组，具体是将此非线性方程组转化为一个高次代数方程和两个线性方程组来求解，归纳为如下三个方程组：

上式包括M-N个方程的线性递推方程组，含N个未知数ai，xi为已知数，因此当M=2N时，利用一般矩阵求逆法，可以求得ai值；如M>2N，则可求得最小二乘意义下的近似值。将求得的ai值代入，解此高次方程，可得N个根Zi，于是被测信号的复极点Si=lnZi/Δt即可求得。将各Zi代入，得线性方程组，利用其中任意N个方程可求得留数Ai。最后，作傅立叶变换，求模，可得到功率谱密度。

其中：

从物理概念上说，Prony技术是一种自适应处理过程，它适时地调整了指数衰减线性组合模型的各个参数与被测信号拟合。与傅立叶周期图截然不同的是，周期图是以固定频率、固定数目的不衰减正弦波的线性组合来拟合。它所存在的缺点与自回归相同，即必须正确估计模型阶数，对噪声也比较敏感，不适合在低信噪比下工作。

1.2.2.5最大似然谱估计法[15] 

由Capon提出的最大似然谱估计器实质上是在最小均方意义下的有限时间脉冲响应滤波器组，它开始用于地震研究中，称为阵列——波数分析，又称Capon谱估计法。

根据信号的谱分布来自适应地更新有限时间脉冲响应滤波器的权系数A，以使在待求频率点上滤波器的频率响应为1(无失真输出)和输出方差为最小，即在

约束条件EA=1之下，使为最小。最后导出结果为：

其中；Rxx为数据xn的

协方差矩阵；T为转置。

以上两式分别为最佳权系数列向量和最大似然谱估计密度估值，上式中只要能求出自相关逆矩阵，就能计算出谱密度估值。最大似然谱估计值与自回归谱估计值存在有如下的关系式：

上式中m为自回归谱估值的阶数。上式说明：m=1，2，…，P各阶白回归谱倒数的平均等于最大似然的倒数，这意味着(1)最大似然谱的分辨率低于最大熵谱的分辨率，这是由于具有低分辨率的低阶最大熵谱和具有高分辨率的高阶最大熵谱的联合平均造成的；(2)最大似然谱在统计上具有更好的稳定性。此外Banggeroer证实，在长数据段情况下，最大似然估计要比自回归谱估计具有较小的方差值。

1.3 参考文献

[1] 刘志刚. 基于现代谱估计理论的信噪分离方法及其应用研究: [学位论文].成都理工大学，2006

[2] Robinson E A. A Historical Perspective of Spectrum Estimation. Proc. IEEE, 1982,    70(Sept):885~907

[3] Schuster A. On the Investigation of Hidden Periodicities with Application to a Supposed 26 Day Period of Meteorological Phenomena, Terr. Mag., 18, 3(1): 13~41

[4] Wiener N. Generalized Harmonic Analysis. Acta Math., 1930, 55: 117~258

[5] Tukey J W. The Sampling Theory of Power Spectrum Estimates, J Cycle Res., 1957, 6: 31~52

[6] Blackman R B, Tukey J W. The Measurements of Power Spectra. New York: Dover Publication, 1959

[7] Kay S M. Modern Spectral Estimation: Theory and Application. Englewood Cliffs. NJ:  Prentice-Hall, 1988

[8] 胡广书. 数字信号处理 理论、算法与实现. 北京：清华大学出版社，2003

[9] Marple S L. A New Autoregressive Spectrum Analysis Algorithm. IEEE Trans. on ASSP, 1980, 28: 441~454

[10] Marple S L. Digital Spectral Analysis with Applications. Englewood Cliffs. NJ: Prentice-Hall, 1987

[11] H. Akaike. Fitting Autoregressive Models for Prediction. Ann. Inst. Statist. Math., 1969, 21: 243~247

[12] H. Akaike, Statistical Predictor Identification. Ann. Inst. Statist. Math., 1970, 22: 203~217

[13] E Parzen. Some Recent Advances in Time Series Modeling. IEEE Trans. Automat. Contr., 1974, AC-19: 723-730

[14] C W Chuang, D L Moffatt. Natural Resonances of Radar Targets via Prony’s Method and Target Discrimination, IEEE Trans. Aerospace Electron. Syst., 1976, AES12: 583~5

[15] Capon J. High-resolution Frequency-wavenumber Spectrum Analysis. Proc. IEEE. 1969, 57(Aug): 1408~1418

[16] Burg J P. Maximum Entropy Spectral Analysis. Ph.D. dissertation Dept. of Geophysics, Stanford Univ., CA, 1975

[17] P F Fougere, E J Zawalick, H R Radoski. Spontaneous Line Splitting in Maximum Entropy Power Spectrum Analysis. Physics Earth Planetary Interiors. 1976, 8(12): 201~207

[18] D N Swingler. A Modified Burg Algorithm for Maximum Entropy Spectral Analysis. Proceedings of the IEEE, 1979, 67(9): 1368~1369

[19] M Kaveh, G A Lippert. An Optimum Tapered Burg Algorithm for Linear Prediction and Spectral Analysis. IEEE Trans. Acoust. , Speech, Signal Processing, 1983, 31

[20] B I Helme, C L Nikias. Improved Spectrum Performance via A Data-adaptive Weighted Burg Technique. IEEE Trans. Acoust., Speech, Signal Processing, 1985, 33

[21] 田坦, 李延. 伯格谱估计算法的一种改进. 数据采集与处理. 2002, 17(3): 276~278

[22] W Campbell, D N Swingler. Frequency Estimation Performance of Several Weighted Burg Algorithms. IEEE Transactions on Signal Processing, 1993, 41(3): 1237~1247

二、研究内容、研究目标以及拟解决的关键问题

    本文研究了求解AR模型参数的Burg算法。Burg算法是一种高分辨率的谱估计算法，因其计算量较小而得到了广泛的应用。然而，大量计算表明，当数据为实数样本时，谱估计的性能不尽人意。特别是在某些数据长度S下，频率估计偏差受信号初相位的影响较大，有时可达到1/(16T)的量级。因此寻求改进算法以减小频率估计偏差仍是近年来研究的课题。

本文的研究目标在于分析传统Burg算法在谱估计时产生偏差的原因，提出提高Burg算法谱估计的准确度的方法，尤其是在实数据、短序列样本的情况下，并且减小信号的初始相位对功率谱谱峰的影响。尽管修正的协方差算法等方法在谱峰偏移方面优于Burg算法，但是运算量偏大，这是制约修正的协方差算法应用的一个重大因素。同理，如果对Burg算法进行改进后，所需的运算量偏大，那么Burg算法也就失去了一个重大优点，因此在对Burg算法进行改进时，必须考虑到是否会显著增加运算量。

本文拟解决的关键问题如下：

（1）分析产生谱偏差和谱线的原因。

（2）改进Burg算法，减小谱偏和减弱谱线。

三、拟采取的研究方案及可行性分析

3.1 有关方法、技术路线

3.1.1 Burg算法

为了克服L-D算法中因估计相关函数给功率谱带来的影响，Burg提出一种新的算法[16]，其基本思想是直接从观测的数据利用线性预测器得前向和后向预测的总均方误差之和为最小的准则来估计预测系数，进而通过L-D算法的递推公式求出AR模型优化的参数。

前、后向预测误差的定义式分别为

将am(k) )分别代入以上两式，得

按前后向预测误差均方(平方)误差的总和，有

并对反射系数求偏导，且之为零，求得总平均方误差为最小时得反射系数等于

上式分母中的两项分别是前向与后向最小方误差，可以证明总的最小平方误差为

按L-D递推公式可求得

综上所述，Burg算法步骤为

(1)根据前、后向预测误差的定义知道

这时均方差等于预测器平均输出功率，故有

(2)按递推关系，估计p阶预测系数(AR模型参数)及最小预测误差的均方值，即

(3)递推高一阶前、后向预测误差，即

(4)求得的AR模型参数估值，得到功率谱估计

Burg算法的优点在于

（1）频率分辨率高

（2）所得的AR模型稳定

（3）计算效率高

然而，Burg方法也有公认的不足。以信号xn=cos(2πfn/fs)+un为例来说明。其中f=100Hz为信号频率，fs=1000Hz为抽样频率，n=1,2,…,N, N=18，即样本数为18。un为均值为0的高斯白噪声，其功率待定。AR模型的阶数待定。

（1）在高信噪比时，它的谱线呈现出[17]。将信噪比设为SNR=80dB，AR模型的阶数设为5。使用原始Burg算法得到的功率谱如图所示。由图可以看到在100Hz频率谱线的右侧又出一条谱线。

谱线在下列情况下最有可能发生：

a) 高信噪比;

b) 正弦分量初始相位为45°的奇数倍；

c) 数据序列长度为正弦分量1/4周期的奇数倍;

d)估计的AR参数的数目与数据的个数相比占有较大百分比;

e) 虚假谱峰常伴随产生谱线现象，这与数据记录短有关，   因此增加数据记录长度这种现象随之消失。

（2）对于高阶模型来讲，Burg方法也可能引入假谱峰。将信噪比设为SNR=10dB，AR模型的阶数设为16。使用原始Burg算法得到的功率谱如图所示。由图可以看出除了100Hz的正确谱线外，还有多处虚假的谱峰。

（3）对于噪声中的正弦信号，Burg方法对正弦信号的初始相位呈现出敏感性。特别是对短的数据记录更是如此。

3.1.2 窗函数

传统Burg算法利用如下公式求预测误差

Swingler [18]指出，采用加权的Burg算法可以减小谱偏。即在上述求解预测误差的公式中，进行简单的加窗处理，如下公式所示

上式中wn,p为窗函数，如三角窗、汉明衰减窗等。

Swingler提出使用的窗函数为

Kaveh 和 Lippert[19]提出使用优化窗函数

Helme和Nikias[20]提出使用数据自适应窗函数

这些窗函数随着递推的进行，即阶数p的增大，窗函数中心也移动，保证有效的利用数据。对比原始Burg算法的预测误差求解公式可以看出，原始Burg算法相当于使用了一个所以权值为1的窗函数。

与原始Burg算法一样，对加权预测误差求偏导并经过推导可得

可以看出加权后算法得到的反射系数与原始反射系数比，分子分母带有都有窗函数。

3.1.3 一阶反射系数[21] [22] 

伯格算法首先从低阶开始，根据数据估计反射系数，然后利用莱文森递归公式计算预测误差滤波器的系数，从而避免了直接计算相关值和矩阵求逆。而Burg算法中一阶反射系数在计算功率谱密度中起着关键的作用。这是因为运用递推公式时，各阶PEF系数均与其有关。然而，Burg算法中一阶反射系数的计算本身就有误差，结果使谱估计的峰值产生偏差，特别是当信号长度为1／4信号周期的奇数倍且初相位为1／4(或其奇数倍)时最为严重。

设观察数据是形如

的正弦波，其中，Δt为采样间隔，φ为初始相位。该信号的自相关函数为

由一阶最大熵方程

可得到系数的精确值为a1,1 = -cosΩ。

但当利用原始伯格算法计算时，由反射系数的求解公式可得

将信号的表达式代入得

由上式可知，这样求出的一阶反射系数a1,1一般情况下不等于精确值-cosΩ。只在下列两种情况上式中第二项为零：

（1）信号长度为半周期的整数倍

（2）信号长度为1/4周期的整数倍，且φ=0或π/4

 而当信号长度为1/4周期的奇数倍且初始相位为π/4的奇数倍时，误差最大。由于计算误差的存在，利用伯格算法所得到的系数均带来误差，结果引起频率估计误差。

改进算法不直接计算一阶反射系数，而是通过二阶PEF系数再求一阶反射系数，结果使得到的一阶反射系数与初始相位无关，从而使谱峰偏差减小，改善了谱估计的质量。

由于一阶反射系数的误差是引起谱估计偏差的主要原因之一，因此从改进一阶反射系数的计算入手，将能改善谱估计的性能。改进算法的基本思想是不用原始的反射系数计算公式直接计算a1,1，而是从二阶误差总功率着手，先计算二阶PEF系数a2,1，a2,2，再计算a1,1。

二阶PEF输出误差功率为

上式中的前、后向预测误差分别为

对a2,1和a2,2分别求偏导得

其中，定义为

解二元一次方程组得

由二阶AR模型系数可以通过递推得到一阶反射系数

对实数据来说，上式可以简化为

将代入简化后的反射系数公式得

可以求得一阶反射系数的值为-cosΩ，且与初始相位无关。由于这种方法求得的一阶反射系数比原始Burg算法求得的值更精确，所以经过递推后所求得的p阶模型参数值也更精确，从而减小谱偏。

3.2 实验手段

    本文中涉及的实验全部都在MATLAB中实现，通过编写并运行相应的实验代码，对算法性能进行评价。

3.3 关键技术

分别使用窗函数和一阶反射系数都能减小谱偏，本文综合了二者，即在求预测误差时，采用加权的Burg算法，对其进行加窗操作；在求解一阶反射系数时，先求二阶模型系数，再通过递推求解一阶反射系数。

由于使用了窗函数，在求解一阶反射系数时，二阶预测误差公式变为

这样，一阶反射系数的表达式形式虽然不变，但是变为如下形式

最终的改进算法计算步骤可归纳如下：

（1）当阶数m=1时，根据下式之一计算一阶反射系数

信号为复数据时

信号为实数据时

其中

（2）当阶数m≥2时，利用步骤（1）得到的a1,1和本文3.1.2节得到的加权反射系数公式

以及Levinson-Durbin递推算法计算其余阶数的AR模型参数。

3.4 可行性分析

    本文利用了MATLAB软件来进行实验，实验结果证明该方法具有一定的不可见性和鲁棒性，而且整个实验过程运行时间不到10秒。因此，该方法在运算成本、运算实时性和适用性方面都是可行的。

3.5 实验结果

    本实验使用四种方法：原始Burg算法、加权法、改进计算一阶反射系数法、综合加权和改进一阶反射系数后的Burg算法，通过分析和对比这四种方法得到的谱估计结果得到综合后的Burg算法的优点和不足。

3.5.1 谱偏问题

（1）信号为无噪声的单频正弦信号

    表达式为

xn=Acos(2πfn/fs + φ)，n=1,2,…,N

其中A=为信号振幅，信号的功率可由A2/2计算得到。f=100Hz为信号频率，fs=1000Hz为抽样频率，φ为初始相位，N为样本数。令阶数p=3。

设φ=45°，N=35，这是最坏的情况，即初始相位为45°的奇数倍，且数据长度为正弦分量1/4周期的奇数倍。分别使用原始Burg算法、加权法、改进计算一阶反射系数法、综合加权和改进一阶反射系数后的Burg算法估计功率谱。算法中的窗函数取汉明衰减窗

这种窗函数的形状如图所示

无噪声时的谱估计结果如图所示

图a, b, c, d分别为原始Burg算法、加权法、改进计算一阶反射系数法、综合加权和改进一阶反射系数后的改进的Burg算法得到的结果。频率偏差分别为2.0996，0.3906，-0.0977，-0.0977。可以看出后三种都比原始算法有改进，改进一阶反射系数法和综合改进算法效果相同，都大大优于原始Burg算法。

（2）信号为有噪声的单频正弦信号

    表达式为

xn=Acos(2πfn/fs + φ)+un，n=1,2,…,N

un是均值为0，功率为1的高斯白噪声，则信噪比SNR=10dB。阶数p=10，其他参数与无噪声时相同。进行30次谱估计的结果分别如下图所示。

图a, b, c, d分别为原始Burg算法、加权法、改进计算一阶反射系数法、综合加权和改进一阶反射系数后的改进的Burg算法得到的结果。计算30次的平均谱偏和方差，得到下表所示数据

（3）谱偏与数据长度N的关系

当信号不含噪声时，令初始相位为0，阶数为3，观察样本数为15~45时，谱偏Δf相对样本数的图像。

由图像可以看出，原始Burg算法和加权后的算法谱偏与数据之间存在周期关系，周期即是样本的周期。但是对比二者的变化幅度可以看出加权后的谱偏幅度明显变小，说明了加权对谱偏的改进作用。而改进一阶反射系数和综合后的Burg算法的谱偏则与数据长度无关，且大大减小。

当采用含有噪声的信号时，令信噪比为SNR=10dB，阶数p=5，得到入下所示的图像。

 由图像可以看出在噪声存在的情况下，四种算法的谱偏都随N的变化而变化，但综合后的Burg算法变化幅度较小。

（4）谱偏与初始相位的关系

当信号不含噪声时，阶数为3，样本数N=40，初始相位的变化范围为0~2π，谱偏Δf相对初始相位的图像。

由图像可以看出，原始Burg算法和加权后的算法谱偏与初始相位之间存在周期关系，周期即是样本的周期。但是对比二者的变化幅度可以看出加权后的谱偏幅度明显变小，说明了加权对谱偏的改进作用。而改进一阶反射系数和综合后的Burg算法的谱偏则与初始相位无关，且大大减小。

当采用含有噪声的信号时，令信噪比为SNR=10dB，阶数p=9，得到入下所示的图像。

    由图像可以看出在噪声存在的情况下，四种算法的谱偏都随初始相位的变化而变化，后三种Burg算法变化幅度比原始Burg算法小，而综合后的算法是最小的。

3.5.2 谱线问题

    由前文所述可知，原始Burg算法在高信噪比时会导致谱线。前述信号在SNR = 80dB时四种方法的谱估计结果如下所示。

（1）阶数p=3时的结果如下图所示，从图中可以看出低阶时都获得好的谱线，无现象。

（2）p=5的结果如图所示，由图可以看出，原始的Burg算法已经出现了谱线，但是另外三种经过改进的算法仍然有好的谱线。

（3）p=10的结果如图所示，由图可以看出，原始的Burg算法得到的谱线得更明显，同时只进行加权改进的Burg算法也出现了，而改进一阶反射系数和综合的Burg算法的谱线仍然没有。

以上的三次模拟说明，加权和改进一阶反射系数都可以减弱谱线，但是相比之下，改进计算一阶反射系数能更好的减弱谱线现象。综合后的Burg算法由于结合了改进一阶反射系数，因此也可以减弱谱线。

3.5.3 虚假谱峰问题

原始Burg算法在高阶时会导致谱线。前述信号在SNR = 10dB，阶数p=30时四种方法的谱估计结果如下所示。由图可以看出，改进后的三种Burg算法同样具有高阶时的虚假谱峰问题。其中，加权算法产生的虚假谱峰幅度较小，但实质上改善不大。

3.5.4 频率分辨能力

    若输入信号为带噪声的两个同幅度的正弦信号，频率分别为100Hz和110Hz，表达式为

xn=Acos(2πf1n/fs)+ Acos(2πf2n/fs)+ un，n=1,2,…,N

其中A=为信号振幅，信号的功率可由A2/2计算得到。f1=100Hz, f2=110Hz为信号频率，fs=1000Hz为抽样频率， N=35为样本数，令阶数p=10。un是均值为0，功率为1的高斯白噪声，则信噪比SNR=30dB。谱估计的结果如下所示。

    由图可以看出，不论是原始的Burg算法，还是三种改进方法，都没有提高分辨能力。

3.5.5 其他窗函数

若算法采用下式所示窗函数

其形状如图所示

固定阶数p时的形状如图所示

若算法采用汉明窗函数，如下式所示

其形状如图所示

固定阶数p时的形状如图所示

用上述窗函数分别加权法和综合算法的汉明衰减窗，再进行谱估计，三种窗函数得到的六个结果如图所示，其30次估计的平均谱偏和方差如下表所示。

3.6 结语

本文提出了综合加权和改进一阶反射系数后的Burg算法。针对原始Burg算法性能方面存在的问题，本实验将其与三种算法：原始Burg算法、加权法、改进计算一阶反射系数法进行对比，表明了本文方法的在原方法上的改进之处。

（1）谱偏问题

对无噪信号，综合法和改进一阶反射系数法性能相同，都优于原始方法和加权法。对有噪信号，综合法在平均值和方差方面都大大优于另外三种方法。

（2）频谱问题

综合法和改进一阶反射系数法性能相同，都优于原始方法和加权法。

（3）虚假谱峰问题

综合法在此问题上没有改善。

（4）频率分辨率

综合法在原始算法基础上没有提高。

    由以上分析可以看出，综合了加权法和改进计算一阶反射系数的Burg算法在谱偏的均值和方差方面有较大改进，在其他方面没有性能的下降，具有应用价值。今后的工作可以考虑使用其他窗函数进行谱估计，并多方面对比性能。