对功率谱估计常用方法的探讨及应用

《对功率谱估计常用方法的探讨及应用》由会员分享，可在线阅读，更多相关《对功率谱估计常用方法的探讨及应用（37页珍藏版）》请在装配图网上搜索。

1、DSP课程的设计对功率谱估计常用方法的探讨及应用分析进行傅里叶变换在频域中研究信号，是研究确定性信号最简单且有效的手段， 但在现代信号分析中，对于常见的随机信号，不可能用清楚的数学关系式来描述， 其傅里叶变换更不存在，转而可以利用给定的N个样本数据估计一个平稳随机信号 的功率谱密度。功率谱估计是数字信号处理的重要研究内容之一。功率谱估计可以 分为经典功率谱估计和现代功率谱估计。本文介绍了各种经典功率谱估计方法，不 仅从理论上对各种方法的谱估计质量进行了分析比较，而且通过Matlab进行了仿 真。在对经典谱估计进行讨论之后，还分析了现代谱估计即参数谱估计方法，通过 观测数据估计参数模型再按照求参

2、数模型输出功率的方法估计信号功率谱。现代谱 估计的内容极其丰富，设计的学科及应用的领域都相当广泛，至今每年都有大量的 科研成果出来。在本文的最后利用现代谱估计的方法讨论了功率谱方法在噪声源信 号识别中的应用。文章还给出了常见谱估计方法的比较，便于深刻理解各种方法的 特点，从而在实际工作中做出合理的选择。1.功率谱方法的发展功率谱估计是随机信号处理的重要内容，其技术渊源很长，而且在过去的40 余年中获得了飞速的发展。涉及到信号与系统、随机信号分析、概率统计、矩阵代 数等一系列的基础学科，广泛应用于人民的日常生活及军事、工业、农业活动中， 是一个具有强大生命力的研究领域。本文将简要回顾一下功率谱估

3、计的发展历程， 对常用的一些方法进行总结。功率谱的估计方法有很多，主要有经典谱估计和现代 谱估计。经典谱估计又可以分成两种：一种是BT法，也叫间接法；另一种是直接 法又称周期图法。现代谱估计的方法又大致可分为参数模型谱估计和非参数模型谱 估计，前者有AR模型、MA模型、ARMA模型、PRONY模型等，后者有最小方 差方法、多分量的MUSIC方法等。1.1功率谱研究的发展过程功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特 征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。下面对谱估计 的发展过程做简要回顾：英国科学家牛顿最早给出了“谱”的概念。后来，1822年，法

4、国工程师傅立叶提 出了著名的傅立叶谐波分析理论。该理论至今依然是进行信号分析和信号处理的理 论基础。傅立叶级数提出后，首先在人们观测自然界中的周期现象时得到应用。19世纪 末，Schuster提出用傅立叶级数的幅度平方作为函数中功率的度量，并将其命名为 “周期图”(periodogram)。这是经典谱估计的最早提法，这种提法至今仍然被沿用， 只不过现在是用快速傅立叶变换（FFT）来计算离散傅立叶变换（DFT）,用DFT 的幅度平方作为信号中功率的度量。1958年，R,Blackman和J.Tukey首先提出BT法，并命名为布莱克曼-杜基谱 估计器（简称BT谱估计器）。这种方法是先按照有限个观测

5、数据估计自相关函数， 再对其求傅里叶变换得到功率谱。在1965年FFT未出现以前，BT法一直是最常用 的方法。周期图较差的方差性能促使人们研究另外的分析方法。1927 年, Yule提出用线 性回归方程来模拟一个时间序列。Yule的工作实际上成了现代谱估计中最重要的方 法参数模型法谱估计的基础。Walker利用Yule的分析方法研究了衰减正弦时间序列，得出Yule-Walker方程， 可以说，Yule和Walker都是开拓自回归模型的先锋。1930年，著名控制理论专家Wiener在他的著作中首次精确定义了一个随机过 程的自相关函数及功率谱密度，并把谱分析建立在随机过程统计特征的基础上，即, “

6、功率谱密度是随机过程二阶统计量自相关函数的傅立叶变换”，这就是 Wiener-Khintchine定理。该定理把功率谱密度定义为频率的连续函数，而不再像 以前定义为离散的谐波频率的函数。1949年，Tukey根据WienerKhintchine定理提出了对有限长数据进行谱估计 的自相关法，即利用有限长数据估计自相关函数，再对该自相关函数球傅立叶变换， 从而得到谱的估计。1958年，Blackman和Tukey在出版的有关经典谱估计的专著 中讨论了自相关谱估计法，所以自相关法又叫BT法。周期图法和自相关法都可用 快速傅立叶变换算法来实现，且物理概念明确，因而仍是目前较常用的谱估计方法。1948

7、年, Bartlett首次提出了用自回归模型系数计算功率谱。自回归模型和线 性预测都用到了 1911年提出的Toeplitz矩阵结构，Levinson曾根据该矩阵的特点于 1947年提出了解Yule-Walker的快速计算方法。这些工作为现代谱估计的发展打下 了良好的理论基础。1965年，Cooley和Tukey提出的FFT算法，也促进了谱估计 的迅速发展。现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方 差性能不好的问题。1967年，Burg提出的最大嫡谱估计，即是朝着高分辨率谱估 计所作的最有意义的努力。虽然，Bartlett在1948年，Parzem于1957年都曾经

8、 建议用自回归模型做谱估计，但在Burg的论文发表之前，都没有引起注意。现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，至今，每年 都有大量的论文出现。目前尚难对现代谱估计的方法作出准确的分类。从现代谱估 计的方法上，大致可以分为参数模型谱估计和非参数模型谱估计，前者有 AR模 型、MA模型、ARMA模型、PRONY模型等；后者有最小方差方法、多分量的 MUSIC方法等。非参数模型谱估计的特点是其模型不是用有限参数来描述，而直 接由相关函数序列得到，这种方法能提高低信噪比时的谱分辨率。参数模型谱估计 是先根据过程的先验信息或者一些假定，建立一个数学模型来表示所给定采样数据 的过程，或

9、者选择一个较好的近似实际模型，而后利用采样数据序列或者自相关序 列，估计该模型的参数，最后把参数代入到该模型对应的理论功率谱表达式，得到 所需要的谱估计。目前大量的论文集中在模型参数的求解上，以求得到速度更快、 更稳健、统计性能更好的算法。1.2功率谱估计方法提出在通信系统中，往往需要研究具有目中统计特性的随机信号。由于随机信号是 一类持续时间无限长，具有无限大能量的功率信号，它不满足傅里叶变换条件，而 且也不存在解析表达式，因此就不能够应用确定信号的频谱计算方法去分析随机信 号的频谱。然而，虽然随机信号的频谱不存在，但其相关函数是可以确定的。如果 随机信号是平稳的，那么其相关函数的傅里叶变换

10、就是它的功率谱密度函数，简称 功率谱。功率谱反映了单位频带内随机信号的一个样本信号来对该随机过程的功率 谱密度函数做出估计。1.3功率谱估计应用及用途功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的 功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所，要求预先 知道信号的功率谱密度(或自相关函数)。例如，在最佳线性过滤问题中，要设计一 个维纳滤波器就首先要求知道(或估计出)信号与噪声的功率谱密度(或自相关函数)。 根据信号与噪声的功率谱(或0 (m)才能设计出能够尽量不失真的重现信号，而把XX噪声最大限度抑制的维纳滤波器。常常利用功率谱估计来得到线性系统的参数估

11、计。例如，当我们要了解某一系 统的幅频特性H (w )时，可用一白色噪声(n )通过该系统。再从该系统的输出样本 y(n)估计功率谱密度P (w )。由于白色噪声的PSD(用P (w )表示)为一常数即 yywwP (w)2，于是有：wwwP (w) = 2 |H (w)|2(1-1)yyw故通过估计输出信号的PSD，可以估计出系统的频率特性H(w)|(模特性)。从 宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这 要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所 谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频 率间隙(距

12、)。提高谱估计的分辨率已成为目前谱估计研究中的一个重要方向功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关 系，实际用途有滤波，信号识别(分析出信号的频率)，信号分离，系统辨识等。 谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。 维纳滤波、卡尔曼滤波，可用于自适应滤波，信号波形预测等(火控系统中的飞机 航迹预判)。2经典谱估计2.1周期图法周期图法又称直接法。它是从随机信号x(n)中截取N长的一段，把它视为能量 有限x(n)真实功率谱S (ejw)的估计S (ejw)的抽样.xx周期图这一概念早在1899年就提出了，但由于点数N般比较大，该方法的 计

13、算量过大而在当时无法使用。只是1965年FFT出现后，此法才变成谱估计的一 个常用方法。周期图法包含了下列二条假设：1. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段 x (n)来估计该随机序列的功率谱。这当然必然带来误差。N2. 由于对x (n)采用DFT，就默认x (n)在时域是周期的，以及x (k)在频域NNN是周期的。这种方法把随机序列样本x(n)看成是截得一段x (n)的周期延拓，这也N就是周期图法这个名字的来历。与相关法相比，相关法在求相关函数 R (m)时将xx (n)以外是数据全都看成零，因此相关法认为除x (n)外x(n)是全零序列，这种处 NN理方法

14、显然与周期图法不一样。但是，当相关法被引入基于FFT的快速相关后，相关法和周期图法开始融合。 比我们发现：如果相关法中M=N，不加延迟窗，那么就和充(N-1)个零的周期图 法一样了。简单地可以这样说：周期图法是M=N时相关法的特例。因此相关法和 周期图法可结合使用。2.2相关法谱估计(BT)这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。它是1958年由 Blackman和Tukey提出。这种方法的具体步骤是：第一步：从无限长随机序列x(n)中截取长度N的有限长序列列x (n)N第二步：由N长序列x (n)求(2M-1 )点的自相关函数R (m)序列。即NxR (m)=xNn=0x (n

15、) x (n + m)N N(2-1)这里，m=-(M-1).,-1,0,1 .,M-1, MN，Rx(m)是双边序列，但是由自相关函 数的偶对称性式，只要求出m=0,。，M-1的傅里叶变换，另一半也就知道了。第三步：由相关函数的傅式变换求功率谱。即S (ejw) =R (m)e-jwm(2-2)xXm=-(M -1)以上过程中经历了两次截断，一次是将x(n)截成N长，称为加数据窗，一次是 将想x(n)截成(2M-1 )长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱 估计，式中的S (ejw)代表估值。一般取MN，因为只有当M较小时，序列傅式x变换的点数才较小，功率谱的计算量才不至于大到

16、难以实现，而且谱估计质量也较 好。因此，在FFT问世之前，相关法是最常用的谱估计方法。当FFT问世后，情况有所变化。因为截断后的x (n)可视作能量信号，由相关N(一 m)卷积定理可得(2-3)这就将相关化为线性卷积，而线性卷积又可以用快速卷积来实现。我们可对上 式两边取(2N-1 )点DFT，则有1 1(2-4)Rx ( m) = N x2 N-1( k ) * x2 N-1( K )=討 X 2 N-1( K )2于是将时域卷积变为频域乘积，用快速相关求&(m)的完整方案如下:1对N长x (n)的充(N-1)个零，成为(2N-1)长的。N2 求(2N-1)点的 FFT，得 X(K)=逻?x

17、(n)W -mk。2 N-12 N-12 N-1X X1 一 N 1 一 N求而N=0(K) 2。由DFT性质，x (n)是纯实的，x (k)满足共轭偶对称,2 N-12 N-12 N-1(K)|2 定是实偶的，且以(2N-1 )为周期。2 N-14.求(2N-1)点的 IFFT：(2-5)Rx (m)二召 迟棵X2N-1(K)叫薦k=-(N-1)这里丄X(K)|2是实偶的,m=-(N-l).O.N-l。本来IFFT求和范围是0至2N-2,N 2 N-1由于丄X(K)|2的实偶性与周期性，求和范围改为-(N-1)至(N-1)不影响计N 2 N-1算结果。同理可将m的范围改为-(N-1)至(N-

18、1)。上述的快速相关中，充零的目的是为了能用圆周卷积代替线性卷积，以便进一 步采用快速卷积算法。快速相关输出是-(N-1 )至(N-1 )的2N-1点，加W (m)窗M后截取的是-(M-1 )至(M-1)的频段，最后作(2M-1 )点FFT，得S”(k)。我们 注意到：如果数据点数与自相关序列点数相同即”=“，则(2N-1 )点的IFFT后紧 跟一个(2N-1 )点的FFT，利用&(m)的对称性，FT运算框的计算式变为S (K)二 * k (m) W -mk X(2_6)XX2 N-1 2 N-1m=-(N-1)由于N=M并假设窗形状是矩形的，第二次W (m)的截断就不需要了。比较式M(3-5

19、)和式(3-6)，S (k)二 FFTR (m) , R (m)二 IFFT丄 X(K)|2正反傅氏xxxN 2 N-1变换可以抵消，直接得S(k) = X(K)2(2-7)xN 2 N-1 丿为了实行基2FFT，也可将(2N-1)点换成2N点，这样做不影响结果的正确性。2.3巴特利特(Bartlett)平均周期图的方法首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当 N Y 时趋于零，达到一致估计的目的。如果X1，X2,XL是不相关的随机变量， 每一个具有期望值卩，方差b 2,则可以证明它们的数学平均x = (x + x +. + x )/L的12l期望值等于卩，数学平均

20、的方 差等 于b 2/L，即:E lx =丄 E lxL 1+ x + + x = Lp =Vark = E1 x - E (x )2=E x 2 - ( ex )22L L= Eyx + x + + x )2L212L所以Var为工ELx LeL e!x i jj=1 i=1j=1i丰jii=1i丰j=Lp (L -1) |H = L |H 2 - m 2x =丄艺丄艺Lpi=1(Ex )(Ex )(2-8)i=1=丄仏L2 - (Ex )2 + LL2 -L21C 2=Lc 2 =-L2L由(2-8)可见，L个平均的方差比每个随机变量的单独方差小L倍。当 LTaVa JxT0，可达到一致谱

21、估计的目的。因而降低估计量的方差的一种有效 方法是将若干个独立估计值进行平均。把这种方法应用于谱估计通常归功于Bartlett。Bartlett平均周期图的方法是将序列x(n) (0 n N 1)分段求周期图再平均。 设将x(n)分成L段，每段有M个样本，因而N = LM，第i段样本序列可写成xi (n) = x(n + iM M)0 n M 1, 1 i M, 0(m)很小，则可假定各段的周期图Ii (o)是互相独立的。对功率xxM谱密度的概念的讨论，谱估计可定义为L段周期图的平均，即(2-9)P (o)=丄Eli (o)xxLMi=1于是它的期望值为P (o)L丄为Exx LME J (o

22、)Mi=1b()丄丄 F p (6 )W (ejg-6)d6xx2 兀xxBsinK _6)M2sin K-6)/2 I2d0(2-10)这里M二N/L，因此Bartlett估计的期望值是真实谱P (o)与三角窗函数的卷XX积。由于三角窗函数不等于5函数，所以Bartlett估计也是有偏估计即Bias丰0，但 当 N T x 时，Bias T 0。由于我们假定各段周期图是相互独立的，所以可按式(3-8)得到下式: (o)11 Vartl (o)L mVar*PxxP 2(o ) 1 +L xx/ s i nOM)、.M s i o 丿VarP (o )h下降的，当 L Ts 时，Varxx(2

23、-11) (o Jt 0。因xx由此可见，随着L的增加此Bartlett估计是一致估计。比较式(2-10)的eP (o)!与式(2-1)的El o)可见在二种情况的估计量的期望值xxNW (evo)的卷积形式，但后者将前者WB中N改为M，主瓣的宽度增大。由于主瓣的宽度愈窄愈接近5函数,B今式(2-10)中W (ejo)的主瓣宽度大于后式中的主瓣宽度，因而 Biasb (o)，而主瓣愈宽显然使分辨率愈差。因此Bias可用来说明谱的 N都是真值P (o)与窗口函数xxBM 二 N/L N。因而使w (evo)B 偏倚BiasPxx分辨率，Bias愈大说明谱分辨率愈差。一个固定的记录长度N,周期图分

24、段的数目L愈大将使方差愈小，但M也愈小，因而使Bias愈大，谱分辨率变得愈差。因此 Bartlett方法中Bias或谱分辨率和估计量的方差间是有互换关系的。M和N的选择 一般是由对所研究的信号的预先了解来指导的。例如，如果我们知道谱有一个窄峰, 同时如果分辨出这个峰是重要的，那么我们必须选择M足够大。又从方差的表达 式我们可以确定谱估计的可接受的方差所要求的记录长度N=(LM)。由此可见Bartlett法使谱估计的方差减小是用增加Bias以及降低谱分辨率的代 价换来的。实际上，当N 定时，Bias与Vr的互换性是谱估计的一个固有特性。为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均

25、值高斯分 布的随机过程，其功率谱密度为1(1 一 az _1)(1 一 az)a = 1 一 0.1这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其H=1/（1-az-1）的滤波器形成的。为了简单设选用N = 8的矩形窗函数。其周期图的 期望值（用EI （w）表示）与真值P （w）均表示在图2.4中。它说明N = 8的周期图可以8xx得到的Bias的情况。图2-2表示M=8分4段与16段二种情况平均后的周期图。显 然L=4的方差比L=16的大。将L=16的曲线与图2.4的曲线比较可见在这种估计中 误差的大部分起因于Bias而不是方差（因为二种情况均为M = 8，Bias相仿，误差也

26、 相仿）。M=16, L=2及8的周期图表示在图2-3。图2-3中L不同造成的影响也是明 显的。但是这二个曲线的起伏都很大，因此有理由认为误差主要起因于方差。比较 M = 8与M=16的周期图可见，M愈大（即2.3中的N愈大），将使周期图的起伏愈 增快的结论，在这里也同样成立。比较图2-2与图2-3发现在这个例子中最好的选择是应用L=16, M=8的估计而不是 L=8、M=16的估计，即宁可减小方差，牺牲Bias。在实际中，当然功率谱密度的 真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预画先知道的。通过改变M和L以及利用预先知道的情况，通常可以很好地进行 选择。平均周期图

27、的方法特别适合于应用FFT算法。因此在FFT出现以后这个方法 比下面将要讨论的利用窗口函数的处理法用得更多。而在FFT出现以前主要是用窗 口函数处理法平滑周期图。图2-1P ( )与EI (o )的特性x图2-2平滑后的周期图（每段取8个数据）2.4 Welch 法Welch提出了对Bartlett法的修正使更适合于用FFT进行计算。他主要提出二 方面的修正，其一是选择适当的窗函数(n)，并在周期图计算前直接加进法，这样 得到的每一段的周期图为1MUn=0xi (n)e(n)e -阿(2-12)这里U =丄七1W2(n)为归一化因子，而Bartlett法每段的周期图为Mn=0xi(n)e-jn

28、=0(2-13)这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。其二是在分段 时，可使各段之间有重迭，这样将会使方差减小(当N与M 一定时)。他建议可以重 迭 50%。3经典谱估计方法编程与分析3.1直接法直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的 序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以 N，作为序列x(n)真实功率谱的估计。Matlab代码示例：clear;Fs=1000; %采样频率n=0:1/Fs:1;%产生含有噪声的序列改变n的取值范围观察图形的变换xn=cos(2*pi*40*n)+3*cos(2*pi*

29、100*n)+randn(size(n);window=boxcar(length(xn); %矩形窗nfft=1024;Pxx,f=periodogram(xn,window,nfft,Fs); %直接法 plot(f,10*log10(Pxx);图3-1周期图法中N=1000功率谱图3-2周期图中N=100000功率谱-50_gQ I|050100150200250300350400450500图3-3周期图当N=100功率谱随着采样点数的增加，该股集是渐进无骗的。从图中可以看出，采用周期突发 估计得出的功率谱很不平滑，相应的估计协方差比较大。而且采用增加采样点的办 法也不能吃周期图变得更

30、加平滑，这是周期图法的缺点。周期图法得出的估计谱方 差特性不好：当数据长度N太大时，扑线的起伏加剧；N太小时谱的分辨率又不好。 对其改进的主要方法有二种，即平均和平滑，平均就是将截取的数据段x (n)再分N成L个小段，分别计算功率谱后取功率谱的平均，这种方法使估计的方差减少，但 偏差加大，分辨率下降。平滑是用一个适当的窗函数W(ejw)与算出的功率谱 Sx (ejw)进行卷积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小，但X分辨率下降。3.2间接法间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便 得到x(n)的功率谱估计。Matlab代码示例： clea

31、r;Fs=1000; %采样频率n=0:l/Fs:l;%产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n);nfft=1024;cxn=xcorr(xn,unbiased); %计算序列的自相关函数CXk=fft(cxn,nfft);Pxx=abs(CXk);index=0:round(nfft/2-l);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1);plot(k,plot_Pxx);图3-4直接发功率谱图当M=N-1时,BT法与周期图法估计出的功率谱是一样的;当M 02 m =

32、 0wm 0(4-4)xx可见， 阶AR模型输出的相关函数具有递推的性质，因而选用AR模型进行谱 估计只需较少的观测数据将式(4 )写成矩阵形式得_ R(0)R (1)R (1)R(0)R (p) R( p1)Ta.1=Q 2w0(4-5)_ R( p)R( p 1)R(0) _aL p0上式就是著名的Yule-Walker( YW)方程.它表明，只要已知观测数据的自 相关数，就能求出AR模型参数 3 和 w，进而按式(3 )求得信号功率谱的估 值。另外，从AR模型的差分方程式可知，该模型的现在输出值是它本身过去 值的回归，这与预测器存在着一定的相似性，它们之间有着非常密切的关系，即它 们的系

33、统函数互为倒数，也就是说预测误差滤波器A (z)就是AR信号模型H(z) p的逆滤波器.因此通过预测误差滤波器优化设计使预测均方误差最小就可求得AR 信号模型的最优参数，即P阶线性预测器的预测系数 a (k) 等于p阶A R模型p的系数a，其最小均方预测误差E等于自噪声方差Q w。p因此，根据上述的Y-W方程以及AR模型与预测误差滤波器之间的关系，就可提取A R模型参数.目前主要有三种：Levinson-Durbin算法、Burg算法和Marple算 法。这三种方法都可以用时间平均代替集合平均的最小平方准则推导得到。在本文 中，我们主要采用TBurg法来估计信号的模型参数，因此主要介绍一下Bu

34、rg法。 Levinson-Durbin 算法a (k) = a 和bm1L -D递推算法是在满足前向预测均方误差最小的前提下，先求得观测数据的自相关 函数，然后利用Y - W方程的递推性质求得模型参数，进而根据式(3 )求得功率谱 的估值.它是模型阶次逐次加大的一种算法，即先计算阶次m=l时的预测系数按此依次计算2，再计算m= 2时的a (1)， a和b 2w122w 2a( p)及b 2当bpwp2满足精度要求时即可停到阶次m=p时的a，a，pp止递推.递推公式为：R (m) +刃 a(k )R (m - k)a (m)=mk=1mkEm1a (k) = a(k) + a (m) a(m

35、- k)mm1mm1其中k=1,2,m-1(4-6)(4-7)=bwm1 a (m)m2 E= r(o)闻 h |a (k)|2m1Lkk=1(4-8)Burg算法Burg算法的基本思想是直接从观测的数据利用线性预测器的前向和后向预测 的总均方误差之和为最小的准则来估计反射系数，进而通过L -D算法的递推公式 求出AR模型优化参数.设观察到的N个数据为x( 0 ),x，x ( N1 )，则 具体算法如下：(1)取 m=1初始化：ef (n)= e (n) =x(n), n=0,1N-10 0/c、1 N1/、b 2 = R(0) = y x2(n)w 0Nn=0(2)计算反射系数:(3)计算滤

36、波器系数及预测误差功率：a (m) = Pmma (k) = a(k) + pmm1(m k)m1k_l,2,m-1E =b 2 =(i-P2)Emwmmm-1递推高一阶前、后向预测误差：ef (n)二ef(n)+ P eb(“ -1)mm-1m m-1eb (n)二 eb (n -1)+p ef()mm-1m m-1把m更新为m+1，重复直到b 2满足要求。pMarple算法Marple算法又称为不受约束的最小二乘法，它的主要思路是为了摆脱因采用递推运 算对确定预测系数的约束，让每一预测系数(模型参数)的确定直接与前、后向预 测的总的平方误差最小(最小二乘法)联系起来.即令总的平方误差pm-

37、1m-1n=pY a (k) x (n - k)2+YY a (k)x(n + k 一 p)2a (0) = 1n=ppLk =0n=ppLk =0p最小.由上式可见，总的平方误差是系数a(k)的函数.若把对各预测系数pppa (k)而非单一地对a(p)=p求导数，并令其为零，就可以得到一组线性方ppp程.解此方程组所得的a(k)就是在最小平方误差准则下的最优预测系数.但由于p方程组系数矩阵不是Toeplitz型，所以不能利用L D算法求解.为了减少运算 量，Ma r p l e提出一种格型结构的高效递推算法。4.1.1 AR模型的稳定性及阶的确定AR(p)模型稳定的充分必要条件是H(z)的极

38、点(即A(z)的根)都在单位圆内。稳 定的AR(p)模型将具有以下的性质：(1) H(z)的全部极点或A(z)的所有根都在单位圆内。(2) 自相关矩阵是正定的。(3) 激励信号的方差(能量)随阶次的增加而递。(4) 反射系数的模恒小于1。但是在实际应用中，levinson算法的已知数据(自相关值)是由 来估计的，有 限字长效应有可能造成大的误差，致使估计出来的AR(p)参数所构成的A(z)的根跑 到单圆上或单位圆外，从而使模型失去稳定。在递推计算的过程中如果出现这种情 况，将致，或|即停止递推计算。通常事先并不知道AR模型的阶。阶选得太低，功率谱受到的平滑太厉害，平 滑后的谱已经分辨不出真实谱

39、中的两个峰了。阶选的太高，固然会提高谱估计的分 辨率，但同时会产生虚假的谱峰或谱的细节。一种简单而直观的确定AR模型的阶的方法，是不断增加模型的阶，同时观察 预测误差功率，当其下降到最小时，对应的阶便可选定为模型的阶。但是预测误差 功率(或AR模型激励源的方差)是随着阶次的增加而单调下降的，因此，很难确定 降到什么程度才合适。另一方面，应该注意到，随着阶次增加，模型参数的数目亦 增多了，谱估计的方差会变大(表现在虚假谱峰的出现)。因此，不能简单地依靠观 察预测误差功率的下降来确定模型的阶。与此相应的另一种简单方法是观察各阶模 型预测误差序列的周期图，当它很接近于平坦(白色谱)时即对应最佳的阶。

40、4.2 MA模型及功率谱估计MA模型是一个全零点模型，其功率谱不易体现信号中的峰值，即分辨率较低。 从谱估计的角度来看，MA模型谱估计等效于经典谱估计中的自相关法，若单纯为了 对一段有限长数据做谱估计，就没必要求解MA模型了。但在系统分析于辨识以及 ARMA谱估计中都要用到MA模型，因此仍有必要讨论MA系数的求解方法。目前提出 的方法主要有：1，谱分解法；2用高阶的AR模型来近似MA模型；3，最大似然法。 以第二种方法为例讨论MA模型参数的提取。i有N点数据x (n)建立一个p阶的AR模型，pq，可用AR模型参数的计 算方法求出p阶的AR系数a (k)，k=l,2.p；pii利用a ，k=1,

41、2p进行线性预测，等效于一个q阶的AR模型，再 一次利用AR系数的求解方法得到b (k)，k=1,2，q。进而通过两次求AR系数可得MA模型的系数。4.3ARMA 模型理论上可以证明，在一定条件下，A R MA和M A模型都可以用一个阶次足够大的 A R模型来近似，所以在实际使用中我们用一个阶次比较高的A R模型代替A R MA 模型来进行功率谱估计，避免了求解非线性方程组带来的困难。5现代谱估计的应用一基于AR模型的噪声源识别5.11研究背景在噪声治理实践中，常常需要鉴别多源噪声中各声源对总体噪声的贡献情况， 以便有针对性地进行治理。但是，在大多数情况下。不允许或无法使噪声源设备单 独运行。

42、例如，油炉燃烧时必须供水，水泵噪声将与油炉燃烧噪声并存；发电厂中 所有设备都不得随意停机等等。因此，迫切需要能够进行在线分析各声源发声情况 的声源鉴别系统。本文就是针对上述问题提出的。使用现代信号分析方法处理采集 到的仪器设备在运行过程中产生的噪声信号，实现对噪声源的识别。噪声已被公认为是与水质污染，大气污染并列的三大污染之一。在我们生活的 环境中每个人都会受到各种噪声的干扰和影响。特别是在交通，工业生产和日常生 活中，噪声更是无处不在。据全国统计，在反映环境污染的投诉中，关于噪声污染的人民来信和来访的件 数逐年增加，已从1991年的2. 78万件增加至1995年的3. 90万件，增加了 40

43、% 以上；而反映噪声污染问题的投诉占环境污染投诉的信访比例则从1991年25%增 加到1995年的35. 6%，五年中增加10个百分点。这一比例高居各类污染投诉的 首位。由于环境噪声污染影响范围较大，近年来因噪声扰民引起的纠纷不断出现， 其中以反映商业、饮食服务业和建筑施工场所噪声扰民居多。5.1.2噪声源识别与控制方法现状噪声控制的目的就是根据实际需要和可能性，用最经济的办法把噪声限制在某 种合适的范围内。为了有效地控制噪声，一般根据噪声传播的具体情况，分别在噪 声源部位，噪声传播途径或噪声接受部位采取措施。目前经常采用的措施是对噪声 传播途径的控制。即在噪声传播的过程中利用吸声材料或吸收结

44、构来吸收声能，用 隔声设旌隔离噪声，或者用消声设备降低空气中声的传播。实际上我们对噪声传播 途径的控制也是基于对主次声源识别的基础上的。识别了主要噪声源，我们才能有 针对性地治理噪声。在目前对噪声源识别的研究中，主要有以下几种方法：5.1.2.1消去法所谓消去法，就是首先把附带全部装备的试验对象（如汽车、发动机等）在一定 条件下测定其总和的工作噪声，然后去除其中的一部分装备或控制这部分噪声传 出，再按同样试验条件测定去除部分装备后的试验对象的工作噪声，则去除这部分 装备前后试验对象噪声交化值，即被视作被拆卸或控制传出装备噪声的方法。这一 方法试验难度大，耗时耗力，由于对部分设备要求停机而不能进

45、行在线测量，因而 不能获得实验对象噪声分布的整体概念.5.1.2.2声强测量法声强测量法是利用双点声压梯度的积分来近似空气质点的振动速度，并利用 FFT来实现声强的实时测量。采用声强测量法分析噪声可获得丰富的噪声信息，不 但可测出噪声的大小，而且还能测出声能的流向。该方法对环境要求低。现场识别 能力强.但是很强的背景声会对声强测试产生不良的影响，产生较大的误差，并且 其测试精度较难控制。5.1.2.3相干函数法相干函数法是根据相关理论，用互相关函数描述被测部件的振动信号与外界噪 声信号之间的关系.根据结果判断噪声是否由该部件振动引起，从而判断出各个声 源.在实际噪声测量中，常常是复杂的多输入单

46、输出系统。而且输入噪声源又可能 是相互依赖的，因此这种方法也不适合判断多输入单输出的系统5.1.2.4频谱分析法频谱分析法，由于声源噪声的形成机理不同。使每个声源的嗓声特点有差异， 能量分布的频谱范围有的位于全频带，有的位于中、低频或高频。因而，在了解各 组成声源噪声频谱特点的情况下，测量出总噪声频谱。再取得各部分声源的频谱， 并将其与总声源频谱比较，就可以分析出各组成部分噪声贡献幅度。从而找出需控 制的主要声源，以有效降低高能频率段的噪声。频谱分析法中，除应用幅值频谱外, 还常应用功率谱。但是由于采用了传统的功率谱估计，所以存在频率分辨率低，旁 瓣泄漏严重的缺点。由于现代谱估计方法为信号建立

47、了一个准确或近似的模型，从根本上摒弃了对 数据序列加窗的隐含假设，因而能够更好地估计信号功率谱.利用现代谱方法估计 功率谱时，不要求噪声源中只有一个信号为高斯分布，这与盲信号分离有很大区别 （盲信号分离要求源信号中只能有一个信号为高斯分布）。只要利用能量随频率的分 布（即功率谱）来求取噪声源的比重。现代谱估计的上述优点为噪声源识别提供了更 为精确的方法。5.2研究内容由于缺少实际的噪声源，在此采用仿真的噪声信号来模拟噪声。我们知道AR 模型的参数与其阶次有关。阶选得太低，功率谱受到的平滑太厉害，平滑后的谱已 经分辨不出真实谱中的两个峰了。阶选的太高，固然会提高谱估计的分辨率，但同 时会产生虚假

48、的谱峰或谱的细节。所以必须选择合适的AR模型的阶次。本文采用 选用Akaike(alC)准则作为确定模型阶的依据。将信号xl,x2按不同的比例混合在一 起组成混合噪声以模仿生活中常见的噪声信号。5.2.1仿真信号的选取N=512;fs = 1;a1=5;a2=3;a3=2;w=2*pi/fs;x1=a1*sin(0.5*w*(0:N-1)+a2*cos(w*(0:N-1).A2+randn(1,N);x2=a2*sin(w*(0:N-1)+a2*cos(0.8*w*(0:N-1).A2+randn(1,N);y=1:512;plot(y,x1);figure(2);plot(y,x2);save x1 x1;save x2 x2;5.2.2模型选取首先判断AR模型的阶数。这里应用的就是A1C准则。我们用仿真信号xl为 例进行了阶数测定，与模型阶的关系曲线如图所示。从图中可以发现模型的AIC存 在最小值，此时信号模型对应的阶数为12。所以选取12阶的模型。clear;load xl.mat;x=x1;N=length(x);M=20;e1=zeros(M,N);e2=zeros(M,N);e1(1,1:N)=x;e2(1,1:N)