AR模型谱估计方法研究及其应用.docx

AR模型谱估计方法研究及其应用摘 要数字信号处理(DSP)重要的应用领域之一，是建立在周期信号和随机信号基础上的功率谱估计。在实际应用中往往不能获得具体信号的表达式,需要根据有限的数据样本来获得较好的谱估计效果，因而谱估计被广泛的应用于各种信号处理中1。本论文研究了功率谱估计的几种常用的方法，包括经典谱估计和现代谱估计的各种方法，且对每种方法的估计质量做了数学推导，并给出仿真程序及仿真图。经典法主要包括周期图法、自相关法，但这两种方法都存在缺陷，即认为观测数据之外的数据都为零，所以对经典法中的周期图法进行了加窗、平均等修正，因此提出了周期图法的改进方法；现代谱估计的方法分类比较多，AR模型法,MA模型法和ARMA模型法是现代功率谱估计中最主要的参数模型，本论文着重讨论了AR模型参数法2。同时论文将通过对经典谱估计和现代谱估计的实现方法及仿真图的比较，得出经典功率谱估计方法的方差性较差，分辨率较低，而现代谱估计的目标正是在于努力改善谱估计的分辨率，因此能得到较好的谱估计效果，为此应用更为广泛3。关键字：数字信号处理，功率谱估计，周期图法，自相关法，AR模型法 ABSTRACTDigital signal processing (DSP) important application of one of the field. Actually, we cant get the expression of a specific signal, so we need to estimate the power spectral of a signal according to some sample data sequences.So spectrum estimation which is widely used in various signal processing. In this thesis, some common methods of Power Spectral Estimation, such as classical spectral estimation and modern spectral estimation, are studied. The quality of each estimation method is derived, simulation program and simulation figure is given. Classical methods of Power Spectral Estimation mainly include the Periodogram and the BT method. But both of them have a common drawback: the data sequences, beyond the area of the observed sequences, are all presumed to zero. So the Windows and the average method are introduced to improve the quality of the Periodogram. Therefore the improvement of The Periodogram estimation method is proposed. The classification of modern spectral estimation methods are more , AR，MA, and ARMA is the most important parameters of modern spectral estimation. This thesis will focus on discussion of AR model parameters method. At the same time , It can be seen from the comparison and realization of classical spectral estimation and modern spectral estimation, classical power spectrum estimation variance is poor, low resolution .The goal of modern spectral estimation is working to improve the resolution of spectral estimation, better results of the estimation of the power spectrum can be obtained, so it is applied more widely. Keywords: digital signal processing, Power Spectrum Estimation, The Periodogram, the BT methods，AR model38目 录摘 要I1 绪论11.1功率谱简介11.2经典谱估计21.3现代谱估计31.4功率谱估计应用及用途42 谱估计简介52.1随机信号简介52.2平稳随机信号72.3估计质量的评价标准103 现代谱估计123.1平稳随机信号的参数模型123.2 AR模型的正则方程与参数计算133.3 MA模型谱估计163.4 ARMA模型谱估计173.5 AR模型功率谱估计实验193.5.1、实验内容193.5.2、实验分析203.5.3、实验结果及分析203.5.4、实验思考243.5.5、实验源代码253.6 AR模型的应用303.7 小结364论文总结36参考文献381 绪论1.1功率谱简介1.功率谱估计技术渊源流长，在过去的几十年获得了飞速的发展。功率谱估计涉及信号与系统、随机信号分析、概率统计、随机过程、矩阵代数等一系列的基础科学，广泛应用于雷达、声纳、通信、地址勘探、天文、生物医学工程等众多领域，其内容、方法不断更新，是一个具有强大生命力的研究领域4。功率谱估计(PSD)是用有限长的数据来估计信号的功率谱, 它对于认识一个随机信号或其它应用方面来讲都是极其重要的, 是数字信号处理的重要研究内容之一，在军事、生物医学、通信等领域得到了较为广泛的应用5。“谱”最早是由英国科学家牛顿提出来的，后来法国工程师傅里叶提出了著名的傅里叶谐波分析理论，该理论至今仍然是我们进行信号分析和处理的理论基础。傅里叶级数首先在观察自然界中的周期现象得到应用，但傅里叶的计算比较复杂，促使人们研制相应的机器来计算傅里叶级数。在19世纪末，Schuster提出傅里叶系数的平方，并命名为周期图，这是经典谱估计的最早提出法，至今仍被人们沿用6。后来，鉴于周期图的起伏剧烈，提出了平均周期图的概念，并提出了在对有限长数据计算傅里叶系数时所存在的边瓣问题，这就是后来我们所熟悉的窗函数的影响。周期图较差的方差性能促使人们研究另外的分析方法。Yule在1927年提出了用线性回归方程来模拟一个时间序列，从而发现隐含在该时间序列中的周期,从而发现了现代谱估计中最重要的方法参数模型法。Walker利用Yule的分析方法研究了衰减正弦时间序列，并得出了在对最小二乘分析中经常应用的Yule-Walker方程。7Yule的工作使人们重新想起了早在1795年Prony提出的指数拟合法，从而Prony方法形成了现代谱估计的又一重要内容。之后又陆续提出了Wiener-khintchine定理、谱估计自相关法BT法等。所有这些都为现代谱估计的发展打下了良好的基础7。2.功率谱估计可以分为经典谱估计(非参数估计) 和现代谱估计(参数估计)。经典谱估计的方法主要方法有自相关估计法和周期图法以及对周期图的改进方法; 现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，方法大致可分为参数模型谱估计和非参数模型谱估计，前者有AR 模型法(最大熵谱分析法)、MA模型，ARMA模型、Prony 指数模型等；后者有最小方差法，多分量的MUSIC方法等8。其中周期图法和AR 模型法是用得较多且最具代表性的方法。从信号的来源分，又可分为一维谱估计、二维谱估计及多维谱估计。从使用的统计量来分，目前大部分工作是建立在二阶矩基础上的，但由于功率谱密度是频率的实函数，缺少相位信息，因此，建立在高阶矩基础上的谱估计方法正引起人们的注意。从信号的特征来分，在这之前所说的方法都是对平稳随机信号而言，其谱分量不随时间变化，对非平稳随机信号，其谱是时变的，近20年来，以wigner分析为代表的时域分析引起了人们的广泛兴趣，形成了现代谱估计的一个新的研究领域。3.在通信系统中，往往需要研究具有统计特性的随机信号。由于随机信号是一持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信号的频谱9。然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。如果随机信号是平稳的，那么其相关函数的傅里叶变换就是它的功率谱密度函数，简称功率谱。功率谱反应了单位频带内随机信号的一个样本信号来对该随机过程的功率谱密度函数做出估计。1.2经典谱估计直接法：又称周期图法，它是把随机序列的N个观测数据视为一能量有限的序列，直接计算的离散傅立叶变换，得，然后再取其幅值的平方，并除以N，作为序列真实功率谱的估计。周期图这一概念早在1989年就提出了，但由于点数N一般比较大，该方法的计算量过大在当时无法使用，在1965年FFT出现后，此法才变成谱估计的一个常用方法。周期图法包含了下列二条假设：认为随机序列是广义平稳且各态遍历的，可以用其一个样本中的一段来估计该随机序列的功率谱,这当然必带来误差。由于对有限序列采用DFT，就默认此有限序列时域是周期的，以及该有限序列在频域是周期的。这种方法把随机序列样本看成是截得一段的有限序列的周期延拓，这也就是周期图法这个名字的来历。间接法：也叫相关法。间接法先由序列估计出自相关函数，然后对进行傅立叶变换，便得到的功率谱估计。周期图法与相关法相比，相关法在求相关函数时将有限长序列以外的数据看做是零，因此相关法认为除有限长序列外是全零序列，这种处理方法显然和周期法不一样。但是，当相关法被引入基于FFT的快速变换后，相关法和周期图法开始融合。改进的周期图法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。其又有以下几种方法：1.Bartlett法Bartlett平均周期图的方法是将N点的有限长序列分段求周期图再平均。2. Welch法Welch法对Bartlett法进行了两方面的修正:一是选择适当的窗函数，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。不同窗函数的welch谱估计在选择窗函数时，一般有如下要求：1) 窗口宽度M要远小于样本序列长度N，以排除不可靠的自相关值；2) 当平稳信号为实过程时，为保证平滑周期图和真实功率谱也是实偶函数，平滑窗函必须是实偶对称的；3) 平滑窗函数应当在是峰值，并且m随绝对值增加而单调下降，使可靠的自相关值有较大的权值；4) 功率谱是频率的非负函数且周期图是非负的，因而要求窗函数的fourier变换是非负的。在经典谱估计中，无论是周期图法还是其改进方法，都存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，这其实是假设了窗以外的数据或自相关函数全为零，这种假设是不符合实际的，正是由于这些不符合实际的假设造成了经典谱估计分辨率较差。另外，经典谱估计的功率谱定义中既无求均值运算又无求极限运算，因而使得谱估计的方差性能较差，当数据很短时，这个问题更为突出，如何选取最佳窗函数、提高频率分辨率，如何在数据情况下提高信号谱估计质量，还需要进一步研究10。1.3现代谱估计现代谱估计与经典谱估计的主要区别就在于，现代谱估计一般采用信号模型法，信号模型法将原始信号视为白噪声通过一系统的输出信号，通过对输出信号的观测，按照一定的准则，求出相应的系统函数，这样再由输入白噪声和以求得的系统函数就很容易得到输出信号的功率谱。由已知白噪声和系统函数求得的输出序列，实际上是对原始观测到的输出信号的两端进行了估计或延拓。数据长度加宽以后，频谱分辨率会得到改善！因此现代谱估计优于经典谱估计。1.4功率谱估计应用及用途功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛的应用于各种信号处理中。在信号处理的许多场所，要求预先知道信号的功率谱密度。例如，在最佳线性过滤问题中，要设计一个维纳滤波器就首先要求知道信号与噪声的功率谱密度，根据信号与噪声的功率谱才能设计出能够尽量不失真的重现信号，而把噪声最大限度抑制的维纳滤波器常常利用功率谱估计来得到线性系统的参数估计。例如，当我们要了解某一系统的幅频特性时，可用一白色噪声通过该系统，再从该系统的输出样本估计功率谱密度，故通过估计输出信号的PSD，可以估计出系统的频率特性。从宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所谓谱估计的分辨率可以粗略的定义为能够分辨出的二个分立的谱分量间的最小频率间隙，提高谱估计的分辨率已成为目前谱估计研究中的一个重要方向。功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波、信号识别、信号分离、系统辨识等。谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计、高阶谱估计等11。2 谱估计简介2.1随机信号简介2.1.1 随机变量随机变量（random variable）表示随机现象（在一定条件下，并不总是出现相同结果的现象称为随机现象）各种结果的变量（一切可能的样本点）。例如某一时间内公共汽车站等车乘客人数，电话交换台在一定时间内收到的呼叫次数等等，都是随机变量的实例。随机变量在不同的条件下由于偶然因素影响，其可能取各种不同的值，具有不确定性和随机性，但这些取值落在某个范围的概率是一定的，此种变量称为随机变量。随机变量可以是离散型的，也可以是连续型的。如分析测试中的测定值就是一个以概率取值的随机变量，被测定量的取值可能在某一范围内随机变化，具体取什么值在测定之前是无法确定的，但测定的结果是确定的，多次重复测定所得到的测定值具有统计规律性。随机变量与模糊变量的不确定性的本质差别在于，后者的测定结果仍具有不确定性，即模糊性。按照随机变量可能取得的值，可以把它们分为两种基本类型：离散型随机变量，即在一定区间内变量取值为有限个，或数值可以一一列举出来。例如某地区某年人口的出生数、死亡数，某药治疗某病病人的有效数、无效数等。连续型随机变量，即在一定区间内变量取值有无限个,或数值无法一一列举出来。例如某地区男性健康成人的身长值、体重值，一批传染性肝炎患者的血清转氨酶测定值等。1.随机变量的分布函数设X是随机变量，对任意实数，事件X<x的概率称为随机变量X的分布函数。记为，即，易知，对任意实数a,b,。分布函数的性质(1) 单调不减性：若, 则；(2) 归一性：对任意实数，且 ，(3) 左连续性：对任意实数x，2.数学期望、方差、标准差定义： ,为的数学期望值，或简称为均值。；以上分别称为X的标准差和方差。若为离散型随机变量，则上述的求均值运算将有积分改为求和。例如，式中的是取值为时的概率12。3.随机向量在某些实际问题中，往往需要同时用两个或两个以上的随机变量来描述试验的结果。设E是一个随机试验, 样本空间是，设和是定义在上的随机变量， 由它们构成的一个向量叫做二维随机向量或二维随机变量。（注: 二维随机向量 的性质不仅与和有关,而且还依赖于这两个随机变量的相互关系。）4.概率密度函数概率密度函数是为了表示瞬时数据落在指定幅值范围的概率。其定义为：。瞬时值小于或等于某值x的概率定义为概率分布函数或累计概率分布函数5.相关函数表征了一个随机过程自身在不同时刻的状态间，或者两个随机过程在某个时刻状态间线性依从关系的数字特征。相关函数是两随机变量之积的数学期望，称为相关性。统计学中用相关系数xy来描述变量x，y之间的相关性，函数的相关系数，简称相关函数：2.1.2随机信号的特征随机信号具有不重复性、不确定性，通常用概率与统计方法研究其中是否存在某些重复、确定的成分。随机过程在某一时刻的均值（一阶矩）可将总体中各样本函数在的瞬时值相加，然后除以样本函数的个数而得到。自相关函数即为随机过程两不同时刻之值的相关性，又称二阶矩。用和两时刻瞬时值乘积的总体平均值得到。自相关函数的性质：(1) 自相关函数是 的偶函数 ；(2) 当 时，自相关函数具有最大值，；(3) 周期信号的自相关函数仍然是同频率的周期信号，但不保留原信号的相位信息。(4) 当随机信号中含有周期信号时，中也必定有周期性分量，且周期相同。(5) 对变化迅速的信号（宽带随机过程），相关的程度在很小时就完全丧失。2.2平稳随机信号2.2.1平稳随机信号的定义平稳信号分严平稳和宽平稳，严平稳的条件在信号处理中太严格，不实用，一般所说的平稳是指宽平稳，满足三个条件：(1) 均值为与时间无关的常数；(2) 均方有界；(3) 自相关函数与信号时间的起始点无关，只和时间差有关（宽平稳信号的方差和均方也是与时间无关的）。(1) 平稳随机过程的定义:如果对于任意n和以及有则称为严平稳随机过程，或称狭义平稳随机过程。(2) 平稳随机过程的数字特征：1) ,平稳随机过程的数学期望与时间无关；2) ，平稳随机过程的方差与时间无关；3) 其中：；4) 。平稳随机过程的数学期望及方差与无关，它的自相关函数和协方差函数只与时间间隔有关；随机过程的这种“平稳”数字特征，有时就直接用来判断随机过程是否平稳。即若一个随机过程的数学期望及方差与时间无关，而其相关函数仅与有关，即我们就称这个随机过程是广义平稳的。(3) 宽平稳随机过程（广义平稳）:若的数学期望为常数，且自相关函数只与有关，则称为宽平稳随机过程，或称广义平稳随机过程。不难看出，严平稳过程一定是宽平稳过程，反之不一定。但对于正态随机过程两者是等价的。本论文若不加特别说明，平稳过程均指宽平稳过程13。2.2.2 平稳随机信号的自相关函数实随机信号的自相关函数定义：，由于平稳随机信号的统计特性与时间的起点无关，设, 则有。所以，平稳随机信号的自相关函数是时间间隔的函数，记为。平稳随机信号自相关函数的性质：设为平稳随机过程，其自相关函数为，自协方差函数，则它们有如下性质：（1）时的自相关函数等于均方差，自协方差函数等于方差， 即 (2) 当平稳随机信号是实函数时，其相关函数是偶函数，即： (3) =0时的自相关函数、自协方差函数取最大值，即 (4) 若X(t)=X(t+T)，则其自相关函数也是周期为T的周期函数，即 (5) 若均值，当时，与相互独立，有 ，即对于零均值的平稳随机信号，当时间间隔很大时，与相互独立，互不相关。2.2.3 平稳随机信号的功率谱1. 功率谱密度定义：设,<<是均方连续的随机过程，称为的平均功率。称为的功率谱密度，简称谱密度。2. 功率谱密度的性质（1）若 ，则 是 的傅里叶变换；(2) SX()是的非负实函数；(3) 实平稳过程的谱密度是偶函数；当 是的有理函数时，其形式必为，其中,为常数，且, ，分母无实根。3.随机序列的功率谱随机序列，它的相关函数满足其功率谱密度 具有如下式子：2.3估计质量的评价标准1.无偏性对于待估参数，不同的样本值就会得到不同的估计值，这样，要确定一个估计量的好坏，就不能仅仅依据某次抽样的结果来衡量，而必须由大量抽样的结果来衡量对此，一个自然而基本的衡量标准是要求估计量无系统偏差。也就是说，尽管在一次抽样中得到的估计值不一定恰好等于待估参数的真值，但在大量重复抽样时，所得到的估计值平均起来应与待估参数的真值相同换句话说，我们希望估计量的均值（数学期望）应等于未知参数的真值，这就是所谓无偏性(Unbiasedness)的要求。定义:设来自总体X的一个样本，是总体参数 的一个估计量，若 ,则称是 的无偏估计量(Unbiased Estimator)。一个估计量如果不是无偏的就称它是有偏估计量。称为估计量的偏差。无偏估计的实际意义就是无系统偏差，估计量是否无偏是评价估计量好坏的一个重要标准，若 ，但有，则称是的渐近无偏估计。 2.有效性比较两个无偏估计量优劣的一个重要标准就是观察它们哪一个取值更集中于待估参数的真值附近，即哪一个估计量的方差更小，这就是下面给出的有效性(Effectiveness)概念。 定义:设与都是总体参数的无偏估计，若 ,则称比 更有效。在的所有无偏估计量中，如果存在一个估计量，它的方差最小，则此估计量应当最好，并称此估计量为 的最小方差无偏估计，也称其为最有效的3.相合性估计量的无偏性和有效性都是在样本容量n固定的情况下讨论的。由于估计量和样本容量n有关，我们自然希望当很大时，一次抽样得出的的值能以很大的概率充分接近被估参数，这就提出了相合性(Consistency)（一致性）的要求。定义：设 是总体参数的估计量，如果对任意都有, 则称是的相合估计量（或一致估计量）。是的相合估计就意味着依概率收敛于.根据大数定律，无论总体X服从什么分布，只要其阶原点矩 存在，则对任意 都有，所以样本的阶原点矩始终是总体阶原点矩的相合估计。 进一步地, 可以证明：只要相应的总体矩存在，矩估计必定是相合估计。特别地， 总是 的相合估计, 样本方差 和样本的二阶中心矩都是总体方差 的相合估计和又都是的相合估计。由相合性定义可以看出，若是的相合估计，当样本容量很大时，一次抽样得到的值便可作为的较好近似值14。3 现代谱估计3.1平稳随机信号的参数模型由上一章讨论可知，经典功率谱估计方法的方差性较差，分辨率较低。方差性能差的原因是无法实现功率谱密度原始定义中的求均值和求极限的运算。分辨率低的原因，对周期图法是假定了数据窗以外的数据全为零，对自相关法是假定了在延迟窗以外的自相关函数全为零。当然，这种假定是不符合实际的，正是由于这些不符合实际的假定产生了经典谱估计较差的分辨率。在第一章已经简洁的介绍了现代谱估计的基本方法，这些方法技术的目标在于努力改善谱估计的分辨率。参数模型法是现代谱估计的主要内容，也是本章讨论的重点，参数模型法的思路如下：（1）假定所研究的过程是由一个输入序列激励一个线性系统的输出。（2）由已知的，或其自相关函数来估计的参数。（3）由的参数来估计的功率谱。是一个因果的线性移不变离散时间系统，当然，它应该是稳定的，其单位抽样响应是确定的。输出序列可以是平稳的随机序列，也可以是确定性的时间序列。若是确定性的，那么是一个冲激序列，若是随机序列，那么应是一个白噪声序列。工程实际中所遇到的功率谱大体分为三种，一种是“平滑”，即白噪声的谱，另一种是“线谱”，这是由一个或多个纯正弦信号所组成的信号的功率谱，这两种是极端的情况；介于二者之间的是既有峰值又有谷值，这种谱称为ARMA谱。显然，由于ARMA模型是一个零极点模型，它易于反映功率谱中的峰值和谷值。不难想象，AR模型易于反映谱中的峰值，而MA模型易于反映谱中的谷值。AR,MA和ARMA是功率谱估计中最主要的参数模型。本章将会详细的讨论AR模型参数的计算、谱的性能及其他算法（如线性预测、最大熵谱估计等）的关系，最后简要给出MA模型及ARMA模型谱估计算法23。3.2 AR模型的正则方程与参数计算3.2.1 正则方程的求导参数模型法功率谱估计的主要思想是：将广义平稳的过程表示成一个输入序列激励线性系统的输出；由已知的或其自相关函数来估计的参数；由的参数估计的功率谱。AR模型又称为自回归模型，它是一个全极点模型，其当前输出是现在输入和过去输入的加权和，表示如下(其中为白噪声序列；p为AR模型的阶数): (4-1) (4-2)由随机信号通过线性系统理论知输出序列的功率谱 (4-3)其中为白噪声序列的方差，因此进行功率谱估计，必需求得AR模型的参数(k=l，2p)及。假定、都是平稳的随机信号，为白噪声，方差为，现在，我们希望建立AR模型的参数和的自相关函数的关系，也即AR模型的正则方程24。将方程（4-1）两边同乘以，并求得均值，最后得到 (4-4)又因为 (4-5)由Z变换的定义，在（4-2）式中，当时，有，综合（4-4）和（4-5）两式，有 (4-6)在上面的推导中，应用了自相关函数的偶对称性。上式写成矩阵形式，即上述两式即是AR模型的正则方程，又称Yule-Walker方程25。3.2.2 AR模型参数求解的典型算法用线性方程组的常用解法(例如高斯消元法)解YuleWalker方程，需要的运算量数量级为，但若利用系数矩阵的对称性和Toeplitz性质，则可构成一些高效算法，LevinsonDurbin算法是其中最著名、应用最广泛的一种，这种算法的运算量数量级为。这是一种按阶次进行递推的算法，即首先以和模型参数作为初始条件，计算模型参数；然后根据这些参数计算模型参数等，一直到计算出模型参数为止，当整个迭代计算结束后，不仅求得了所需要的P阶AR模型参数，而且还得到了所有各低阶模型的参数。根据线性预测理论知：一个P阶AR模型的个参数同样可用来构成P阶的最佳线性预测器，其预测的最小均方误差等于AR模型激励白噪声的能量，即AR模型是在最小方差意义上对数据的拟合。“前向预测”是利用n之前的P个值对性预测，如公式(4-8)、(4-9)、(4-10)所示；与之对应的“后向预测”公式为(4-11)、(4-12)、(4-13)，其中为为预测误差，P预测误差功率，f表示前向预测，b表示后向预测。模型参数算法就是基于上述最小均方误差时由模型参数估计信号功率的方法，主要有以下几种经典算法：自相关法(BT法)。用自相关法进行功率谱估计，但估计时令前向预测误差功率最小，即对前后都加窗构成，WienerHopf方程系数为Toepli tz矩阵，使用LevinsonDurbin算法可方便快速的求解AR系数。因此自相关法也是已知所有AR系数求解方法中简单的一种，但谱分辨率相对较差26。Burg算法。用Burg算法进行功率谱估计时令前后向预测误差功率之和最小，即对、前后都不加窗，使用LevinsonDurbin递推可快速的求解AR系数。Burg算法是建立在数据基础之上的，避免了先计算自相关函数从而提高计算速度；是较知为通用的方法，计算不太复杂，且分辨率优于自相关个法，但对于白噪声加正弦信号有时会出现谱线分裂现象。计算步骤如下：由初始条件 ，再由式 (4-7)求出；由 得时的参数,求出 和，再由估计；依照Levinson递推关系，求时的参数 及；重复上述过程，直到，求出所有阶次时的AR参数。 (1) 改进协方差算法。同Burg算法一样，改进协方差算法进行功率谱估计时令前后向预测误差功率之和最小，即对、前后都不加窗，但得到的协方差矩阵不是Toeplitz矩阵，因此正则方程不能用Levinson递推算法求解。Marple于1980年提出实现协方差方程求解的快速算法，大大提高了谱估计的性能27。3.3 MA模型谱估计给出模型的三个方程 由（4-1）得 将上式两边同乘以，并求均值，得 （4-8）式中。因为 (4-9)对模型，由式(2) 式得 ; 所以，可以求出模型的正则方程，即有 (4-10)的功率谱为 (4-11)等效于经典谱估计中的自相关法，即MA谱估计等效为信号长度为的自相关法谱估计。 3.4 ARMA模型谱估计ARMA(p,q)模型的差分方程 (4-12)式中。类似地，可导出其正则方程如下： (4-13)式中是系数和的函数，前个方程是高度非线性的。从第个方程开始是线性的，可以解出AR部分的系数，将上式中的第二个方程写成如下展开形式：上式虽然可解出AR部分的系数，但存在以下两个问题：由于式中的真实自相关函数是未知的，因此只能使用估计值来代替，且要用到大延迟的估计值（最大延迟是），而对于给定的信号长度，这将造成估计很不准确。因而，也就不能得到AR部分系数的准确估计。式中阶次和都是未知的，需要事先指定。实际上是式中自相关阵的维数，和决定了的选用范围。因此和的不正确指定有可能导致自相关阵出现奇异。因此，在实际应用中，对自相关阵采用更一般的形式，即取L个方程，这里 ，即，式中，由此得到的最小二乘解为 求得ARMA(p,q)模型中的AR参数，余下的任务就是求解MA部分的参数30。利用求得的AR系数先得到一个FIR系统为序列经此FIR系统滤波，得到一个输出序列，ARMA(p,q)模型与FIR系统级联，近似于模型。因此，可以利用输出序列估计自相关序列并按MA(q)模型谱估计公式来得到MA谱，即，得到MA谱估计后，利用下式即可求得ARMA谱估计：3.5 AR模型功率谱估计实验3.5.1、实验内容AR过程的线性建模与功率谱估计。考虑AR过程：是单位方差白噪声。(a) 取b(0)=1, a(1)=2.7607, a(2)=-3.8106, a(3)=2.6535, a(4)=-0.9238，产生x(n)的N=64个样点。(b) 计算其自相关序列的估计，并与真实的自相关序列值相比较。(c) 将的DTFT作为x(n)的功率谱估计，即：。(d) 利用所估计的自相关值和Yule-Walker法(自相关法)，估计和的值，并讨论估计的精度。(e) 用(d)中所估计的和来估计功率谱为：。(f) 将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。(g) 重复上面的(d)(f)，只是估计AR参数分别采用如下方法：(1) 协方差法；(2) Burg方法；(3) 修正协方差法。试比较它们的功率谱估计精度。3.5.2、实验分析1、 计算真实的自相关值时，采用逆Levinson-Durbin递归方法，由a、b参数得到，其中为滤波器的阶数，再采用公式外推得到的自相关值；2、 实际功率谱，可调用Matlab中的FFT算法得到；3、 自相关序列的估计值采用公式得到；4、 采用各种功率谱估计方法对功率谱进行估计。3.5.3、实验结果及分析仿真参数设置：采样点数为64，频域采样点数为128自相关序列的估计与真实自相关序列值的比较见图1，由图可知估计值与真实值存在一定的误差，但整体变化趋势相差不大。图1 自相关序列题目（c）中功率谱的估计方法实际为周期图法，周期图法估计的功率谱与自相关法估计的功率谱的比较见图2，由图可知，周期图能辨认出两个峰值，而自相关法不能，说明周期图的分辨率大于自相关法。图2 周期图法和自相关法得到的功率谱图3图7的（a）部分分别为采用周期图法、自相关法、协方差法、Burg方法、修正协方差法进行功率谱50次估计的交叠图，（b）部分给出了其整体平均及真实的功率谱。由这些图可以看出，对于这一AR（4）过程，除自相关法外，所有估计都能分辨出两个峰值，且峰值的位置大致相似。此外，周期图法的方差大于其它估计方法。（a）(b)图3 周期图法估计AR(4)过程的功率谱（a）(b)图4 自相关法估计AR(4)过程的功率谱（a）(b)图5 协方差法估计AR(4)过程的功率谱（a）(b)图6 Burg法估计AR(4)过程的功率谱（a）(b)图7 修正协方差法估计AR(4)过程的功率谱表1为采用自相关法、协方差法、Burg方法、修正协方差法得到的a参数和b参数，表中平方误差和的计算公式为，从表中可以看出，自相关法估计的参数与真实值相比相差较大，协方差法、Burg方法、修正协方差法参数估计的性能相当。a(1)a(2)a(3)a(4)b(0)平方误差和真实值2.7607-3.81062.6535-0.923810自相关法1.5371-1.32120.36348-0.134374.728913.5619协方差法2.7324-3.73732.5794-0.88852设为10.0129Burg方法2.7083-3.66992.5098-0.8569设为10.0477修正协方差法2.7071-3.66752.5069-0.85566设为10.0495表1 各种方法估计的a参数和b参数3.5.4、实验思考当观测数据存在观测噪声时，即，其中是单位方差的白噪声，与不相关，考虑观测噪声对各种谱估计方法的影响。仿真参数设置：采样点数为64，频域采样点数为128。图8图11为存在观测噪声时，各种方法对AR（4）过程的功率谱估计，由图可知，观测噪声对各种谱估计方法的估计性能具有一定的影响，在上述仿真参数设置下，周期图法和协方差法可辨认出两个峰值，而Burg方法和修正协方差法不能。不过，当增加频域采样点数时，Burg方法和修正协方差法也能辨认出两个峰值，但频域采样点数的增加意味着计算量的增加。图8 观测噪声下周期图的功率谱估计图9 观测噪声下协方差法的功率谱估计图10 观测噪声下Burg方法的功率谱估计图11 观测噪声下修正协方差法的功率谱估计3.5.5、实验源代码clc;clear;clf;%a参数a1=2.7607;a2=-3.8106;a3=2.6535;a4=-0.9238;p=4;%A=1 -a1 -a2 -a3 -a4;b0=1;%b参数