【学习课件】第四章功率谱估计

1、 对于确定性信号，傅里叶变换是在频率分析研究的理论基础，但对于随机信号，其傅里叶变换不存在，因此研究它的功率谱。第四章功率谱估计4.1 引言 在实际应用中，通常只能采集或观测到平稳随机过程的一个抽样序列的一段(有限个)数据，如果根据这有限个已知数据来估计随机过程的功率谱的问题，简称谱估计(谱分析)问题。对信号和系统进行分析研究、处理有两类方法：一类是时域进行，维纳-卡尔曼滤波和自适应滤波都属于时域处理方法；另一类是频域方法。1精选课件按照Weiner-Khintchine定理，信号的功率谱和其自相关函数服从一对傅里叶变换关系，称为功率谱的定义。对于平稳随机信号，服从各态历经性，集合平均可以用时

2、间平均代替，可推出功率谱的另一个定义。(4.1.1)(4.1.2)(4.1.3)将(4.1.3)式中的集合平均用时间平均代替，得(4.1.4)2精选课件将(4.1.4)式代入(4.1.1)式，得到令 ，则 是随机变量，必须对 取统计平均值，得到将(4.1.4)式代入(4.1.1)式，得到上式被认为是功率谱的另一定义。(4.1.5)(4.1.6)3精选课件谱估计方法经典谱估计，也称为线性谱估计现代谱估计，也称为非线性谱估计经典谱估计BT法：1958年，R.Blackmant和J.Tukey提出， 先估计自相关函数，再计算功率谱。周期图法：1898年，Schuster利用傅里叶级数去拟合待分析的信

3、号，提出周期图的术语，但直到FFT出现，周期图法才受到人们的重视。这种方法直接对观测数据进行FFT，取模平方，除以N得到功率谱。经典谱估计致命的缺点是频率分辨率低，原因是傅里叶变换域是无限大，而观测数据是有限长，观测不到的数据被认为是0。这相当于将信号在时域加了矩形窗，在频域使真正的功率谱卷积一个sinc函数。4精选课件现代谱估计如果由观测数据能够估计出信号模型的参数，信号的功率谱可以计算出来。谱估计问题变成了由观测数据估计信号模型参数的问题。模型种类很多，如AR模型、MA模型等。合适地选择模型，功率谱估计质量比经典谱估计的估计质量有很大提高。参数模型法以信号模型为基础。参数模型法：AR、MA

4、、ARMA模型非参数模型法：Pisarenko、MUSIC、ESPRIT法5精选课件图a BT法图b 最大熵谱估计法图c Pisarenko谐波分解法(自回归PSD法)6精选课件1. 无偏自相关函数的估计4.2 经典谱估计4.2.1BT法BT法是先估计自相关函数，再进行傅里叶变换得到功率谱。下面介绍自相关函数的估计。利用观测到的实随机序列 ，估计自相关函数的两种方法是：无偏自相关函数估计和有偏自相关函数估计。估计器为7精选课件也可写成一个表达式估计性能分析：估计量的偏差：，这是一种无偏估计。估计量的方差：8精选课件可以证明，一般观测数据量N很大， 此式表明，只有当 时，估计量的方差才趋于0，此

5、估计是渐近一致估计。但是当 时，方差将很大，因此，这种方法在一般情况下不是一种好的估计方法。虽然是无偏的，但不能算是一致的。9精选课件2. 有偏自相关函数的估计有偏自相关函数用 表示，估计器为估计性能分析：估计量的偏差：因为 是无偏估计，两边取均值，得有偏估计10精选课件由 可见，只有当 时， 才是无偏的，其它 都是有偏的，但当 时， ，因此 是渐近无偏的。偏移量为估计量的方差：由两种估计的关系式知估计量 的方差为11精选课件当 时， ，并且 虽然是有偏估计，但是渐近一致估计，估计量的方差小于 的方差。实际应用中多用这种有偏自相关估计。也用符号 表示。12精选课件BT法功率谱估计采用有偏自相关

6、函数估计法，对上式进行傅里叶变换，得到BT法的功率谱估计值为为了减少谱估计的方差，经常用窗函数 对自相关函数进行加权，此时谱估计公式为式中加权协方差谱估计要求加窗后的功率谱仍是非负的，这样窗函数的选择必须满足一个原则，即它的傅里叶变换必须是非负的，如巴特利特窗13精选课件为了采用FFT计算 ，设FFT的变换域为 ，必须将求和域 移到 ，功率谱的计算公式如下：这里自相关函数是渐近一致估计，但经过傅里叶变换得到功率谱的估计却不一定是渐近一致估计，可以证明它是非一致估计，是一种不好的估计方法。14精选课件4.2.2周期图法功率谱的另一种定义式重写如下：忽略上式中求统计平均的运算，观测数据为 ，便得到

7、周期图法的定义：上式的绝对值符号内的部分可以用FFT计算，这样周期图法的计算框图如下图所示。15精选课件下面证明周期图和BT法的等价关系，再分析周期图的估计质量。1.周期图和BT法的等价关系令 ，即 ，则上式中的方括号部分正是有偏自相关函数的计算公式，则等价16精选课件已知自相关函数的估计值 ， ，求功率谱的统计平均值，2.周期图法谱估计质量分析1)周期图的偏移17精选课件 称为三角谱窗函数。 式表明，周期图的统计平均值等于它的真值卷积三角谱窗函数，因此周期图是有偏估计，但当 时， ，三角谱窗函数趋近于 函数，周期图的统计平均值趋于它的真值，因此周期图属于渐近无偏估计。式 中乘积的傅里叶变换，

8、在频域服从卷积关系，得到式中式中18精选课件2)周期图的方差 由于周期图方差的精确表示式很复杂，为分析简单，通常假设 是实的零均值的正态白噪声信号，方差是 ，即功率谱是常数 ，其周期图用 表示， 表示观测数据的长度，则下面先求周期图的均值，再求其均方值：19精选课件式中上式说明周期图是无偏估计，这是由于假设信号是实白噪声。求均方值时，先求两个频率 和 处的均方值，最后令20精选课件利用正态白噪声、多元正态随机变量多阶矩公式，有将上式代入周期图的均方值公式中，得21精选课件将 代入上式，得显然，当N趋于无限大时，周期图的方差并不趋于0，而是趋于功率谱真值的平方，即无论N取多长，周期图的方差总是和

9、 同一数量级。信号功率谱真值是 ，说明周期图的方差很大，周期图的均方误差也非常大。用这种方法估计的功率谱在 附近起伏很大，故周期图是非一致估计，是一种很差的功率谱估计方法。显然，当N趋于无限大时，周期图的方差并不趋于0，而是趋于功率谱真值的平方，即22精选课件图 白噪声的周期图23精选课件则第 组的周期图为： 假设随机信号 的观测数据区间为： ，共进行L次独立观测，得到L组记录数据，表示为 ，4.2.2经典谱估计方法改进一般有三种改进方法：平均周期图法、窗函数法和修正的周期图平均法。1.平均周期图法平均周期图的思想：对一个随机变量进行观测，得到L组独立记录数据，用每一组数据求其均值，然后将L个

10、均值加起来求平均。这样得到的均值，其方差将是用一组数据得到的均值的方差的1/L。24精选课件可见，平均周期图仍然是有偏估计，偏移和每一段数据的个数M有关；由于 ，平均周期图的偏移比周期图的偏移大。由于平均周期图三角谱窗主瓣的宽度变宽，分辨率更加降低，因此，偏移的大小反映分辨率的低与高。偏移越大分辨率越低。将得到的L个周期图进行平均，作为信号 的功率谱估计，有为分析偏移，对上式求统计平均，得前边已求出周期图的统计平均值为25精选课件但实际中数据是连续的，数据段之间不能认为是完全不相关的，估计方差的减少小于 ，或认为近似为原来的 。平均周期图的方差是周期图的方差的 。显然是以分辨率的降低换取了估计

11、方差的减少，当然，估计的均方误差也减少。 求平均周期图的方差，由于是L次独立观测，L个周期图相互独立，因此平均周期图的方差为实际中，很难得到独立的数据组，一般常用的方法是将长度为N的数据分成L段，每段有M个数据，N=LM，第i段数据表示为：如果各段数据互不相关，利用平均周期图法估计功率谱的公式、偏移和方差的计算公式和上边结果一样，方差是周期图法的26精选课件按照 式，段数L愈多，方差愈小，即功率谱愈平滑，如果 ，则 。但如果数据量N一定，L加大，每一段的数据量M将减少，因此估计量方差减少了，使偏移加大，分辨率降低。估计量的方差和分辨率是一对矛盾。可根据实际情况适当选择L和M，如果对分辨率要求不

12、高，可以取L大些。最好将N取大些，分辨率和估计误差都能适当满足要求。下图中，信号是均值为0，方差为1的白噪声，观察数据长度为N=256，L=2,4,8段，按平均周期图法估计功率谱曲线如图所示。27精选课件表明，随着分段数的增加，功率谱估计值在1附近摆动的幅度越来越小，显示出分段平均对周期图方差减少有明显效果。图 平均周期图法28精选课件2.窗函数法 这种方法是用一适当的功率谱窗函数 与周期图进行卷积，来达到使周期图平滑的目的。即式中 是有偏自相关函数。那么此式和前面讲的BT法的加权协方差谱估计方法是相同的。29精选课件 两个窗函数有关。将 和 代入上式，得在 式中，周期图和谱窗函数卷积得到功率

13、谱，等效于在频域对周期图进行修正，使周期图通过一个线性非频变系统，滤除掉周期图中的快变成分，谱窗函数需具有低通特性。求 的统计平均，得上式表明，周期图的窗函数法仍然是有偏估计，其偏移和 、 30精选课件比较上式和周期图法的 式，等于将 式中的窗函数 改变成 。 ，则如果 窗的宽度比较窄，M比N小得多， ，则如果 比 的主瓣宽度宽，利用窗函数法可进一步平滑周期图，减少估计方差，但是偏移加大，分辨率降低。由于功率谱是频率的非负函数，要求 是非负函数，三角窗函数是适用的，哈明窗、汉宁窗并不适用。31精选课件3.修正的周期图求平均法 这种方法和平均周期图法一样，首先把数据长度为N的信号 分成L段，每一

14、段数据长度为M，N=LM，然后把窗函数 加到每一个数据段上，求出每一段的周期图，形成修正的周期图，再对每一个修正的周期图进行平均。第i段的修正周期图为式中U称为归一化因子。将每一段的修正周期图之间近似看成互不相关，最后功率谱估计为32精选课件对上式求统计平均，得式中 这种在计算周期图之前，先对各数据段加窗函数的方法，使平均周期图方法的估计方差减少，分辨率降低。但这种方法对窗函数没有限制；分段时，相邻的两段可以重叠，进一步使方差减少，可以重叠50%。为使功率谱估计渐近无偏，归一化因子U是不可缺少的。33精选课件由图可见，加哈明窗，大大压低了旁瓣，使低电平信号清晰可见，但由于主瓣加宽，功率谱波峰变

15、宽了，降低了信号的分辨率。一般来说，两个等幅的正弦信号的频率相隔很近，可以不加数据窗，频率间隔应该大于 才能分辨。图 利用数据窗减少窄带过程周期图的旁瓣(a)没有数据窗 (b)加哈明窗34精选课件 (1)选择合适的信号模型： (2)根据 有限的观测数据，或者它的有限个自相关函数估计值，估计模型的参数； (3)计算模型输出功率谱。 如果能确定信号 的信号模型，根据信号观测数据求出模型参数，系统函数用 表示，模型输入白噪声方差为 信号的功率谱可用下式求出4.3 现代谱估计中的参数建模功率谱估计可分成三个步骤： 这种方法隐藏着数据的外推，不是假设观测区以外的数据为0，它的估计分辨率都比经典谱估计高，

16、故现代谱估计也统称为高分辨率谱估计。35精选课件 主要考虑原则：模型能够表示谱峰、谱谷能力。对于具有尖峰的谱，应该选用具有极点的模型，如AR和ARMA模型；对于具有平坦的谱峰和深谷的信号，可以选用MA模型；既有极点也有零点的谱应该选用ARMA模型。对于没有模型准确反映其谱特性的信号，可以选用高阶的AR模型近似表示。 在选择模型合适基础上，应尽量减少模型的参数。4.3.1模型选择下图表示的是用MA模型估计二阶AR信号功率谱的例子。36精选课件图 对AR(2)信号的模型选择37精选课件假设模型的差分方程和系统函数分别用下式表示：4.3.2模型参数和自相关函数之间的关系式中，模型输入白噪声 均值为0

17、，方差为 ，B(z)称为MA部分，是一q阶多项式，单位脉冲响应为 ， ， A(z)称为AR部分，是一个p阶多项式， ， ，38精选课件将上式两边同乘以A(z),得1.ARMA模型的系数和信号自相关函数之间的关系将上式进行Z反变换，左边的反变换为公式右边的Z反变换为式中39精选课件此式即是ARMA模型输出自相关函数与模型参数之间的关系式。如果能由信号的观测数据估计出信号的自相关函数，即可求出ARMA模型的参数，再由 求出信号功率谱。由于模型H(z)是因果的， ， ，可以得到代入上式得这样得到40精选课件 ， ，然后再设法求出MA部分的系数。当 时，上式是一个线性方程，用矩阵表示如下上式共有p个方

18、程。可以用该方程首先计算出AR部分的p个系数41精选课件AR模型的系统函数 ，相当于ARMA模型中 的情况，这样公式中2.AR模型的系数和信号自相关函数之间的关系因为 是因果的，因此 时 ，得42精选课件 式的矩阵形式如下也可以将上式中 的情况写成矩阵形式：43精选课件 称为自相关矩阵，满足 ，是一个Hermitian矩阵，也是一具Toeplitz矩阵，是正定矩阵。令或者用模型参数表示44精选课件 的情况，此时， ，这样可得MA模型系数和信号相关函数的关系为3.MA模型的系数和信号自相关函数之间的关系MA模型的系统函数 ，相当于ARMA模型中 ，45精选课件4.4 AR模型谱估计的性质4.4.

19、1AR的线性预测 第二章，已推导出维纳线性一步预测器系数和信号自相关函数之间的关系式，这就是著名的Yule-Walker方程。式中， 表示线性一步预测误差，公式为 表示预测器的系数，它和线性预测器的 关系为46精选课件令 ，由上式，得对 式进行Z变换，得称 为线性一步预测误差滤波器，其作用是将 转换成预测误差 ，如图所示。一般认为 具有白噪声的性质，因此 也称为白化滤波器。图 预测误差滤波器47精选课件对比AR模型和预测误差滤波器的系统函数当 时， 和 互为逆滤波器， 。基于以上分析，可得AR模型定义为式中， 是基于信号前p个样本的最佳一步线性预测，为 是模型输入白噪声。48精选课件对比当AR

20、模型的阶数与线性预测器的阶数相同时，得到信号自相关函数和它的AR模型参数之间的关系服从Yule-Walker方程。因此AR模型法也称为线性预测AR模型法。49精选课件如果 具有白噪声性质，可以利用AR模型与预测滤波器之间的关系对上式进行解卷积，得到 信号。方法是：先由 的观测数据估计AR模型的参数，得到AR模型的系统函数 ，再让 通过 的逆滤波器 ，便得到信号 。AR模型与线性预测之间的关系，可以被用来解卷积，假设信号 通过一个AR系统，系统单位脉冲响应 ，响应 ，50精选课件 AR模型必须因果稳定，即极点均在单位圆内，才能保证信号 是平稳随机信号，于是AR模型 和预测误差滤波器 互为逆滤波，

21、那么 应为最小相位系统。 但是由解Yule-Walker方程得到AR模型的参数，其极点不一定在单位圆内。 下面证明当最佳p阶线性预测系数与AR模型参数相同时，由此得到的极点保证在单位圆内，AR滤波器稳定，预测误差滤波器 或者 是最小相位系统。 4.4.2预测误差滤波器的最小相位系统51精选课件设 的第i个零点 在单位圆外部，即 ，用 代替 ，这时 的幅度函数不变，按照上式计算出的观测误差滤波器输出功率 应不变，仍是最小。下面将用证明预测误差滤波器输出功率 不是最小的，这一矛盾的结论只能说明 的零点不可能在单位圆外部。式中 ， 是最佳预测系数。 解Yule-Walker方程得到的是最佳线性预测滤

22、波器的系数，此时预测误差滤波器输出功率 达到最小，用 表示，即 下面用反证法证明 的全部零点不在单位圆外部。52精选课件说明将零点 换成 时，预测误差滤波器的输出功率将更小,与假设矛盾。 零点必须在单位圆内或单位圆上。又因为自相关矩阵是正定矩阵， 的零点必须在单位圆内部。因为 ，则式中 ，代入 ，得式中53精选课件对于 情况，公式本身就是一个递推方程，如果已由观测数据计算出p+1个自相关函数，用 ， 表示，对于 的情况，可以用该公式外推得到，AR模型的自相关函数和模型系数服从Yule-Walker方程4.4.3 AR模型隐含自相关函数延拓特性AR模型的自相关函数隐含着自相关函数外推的特性，它具

23、有很高分辨率。BT法中自相关函数只由观测数据计算出有限个值，其它认为是0。54精选课件4.5 AR谱估计的方法 已知信号的观测数据，估计功率谱的方法：自相关法、伯格法、协方差和修正协方差法。 它们共同的特点是由信号观测数据求模型系数时，均使用信号预测误差最小的原则。 对于长记录数据，这些方法的估计质量相似；但对于短记录数据，不同的方法之间存在差别。4.5.1 自相关法-列文森(Levinson)递推法自相关法的出发点是选择AR模型的参数使预测误差功率最小。55精选课件假设信号 的数据区在 范围，有 个预测系数， 个数据经过冲激响应为 的滤波器，输出预测误差 的长度为 ，因此有预测误差功率为 的

24、长度长于数据的长度，上式中数据 在 以外补充零点，相当于对无穷长的信号加窗处理，会引入误差。56精选课件上式对系数 的实部和虚部求微分使预测误差功率最小，得式中自相关函数采用有偏自相关估计，即Yule-Walker方程57精选课件Levinson-Durbin算法的推导：最直观的方法是根据AR(1)、 AR(2)、 AR(3)各阶模型的Yule-Walker方程的求解结果归纳出一般的迭代计算公式。第二种方法是利用自相关矩阵对称性和Toeplitz性质，通过构造增广矩阵导出。 用线性方程组的常用解法(如高斯消元法)求解Yule-Walker方程，运算量数量级为 。利用自相关矩阵的对称性和Toep

25、litz性质，可构成一些高效算法(Levinson-Durbin算法)，运算量数量级为 。Levinson-Durbin算法：是一种按阶次递推的算法。它以AR(0)和AR(1)模型参数作为初始条件，计算AR(2)模型参数；再用AR(2)模型参数计算AR(3)模型参数，k阶模型参数由k-1阶模型参数计算得到。一直计算出AR(p)模型参数为止。58精选课件一阶AR模型(p=1)的Yule-Walker方程为由该方程解出实平稳信号59精选课件由上面方程解出增加一阶，即令p=2，得60精选课件 称为反射系数， 。 ，随着阶数增加，预测误差功率将减少或不变。然后令p=3,4,.,以此类推，可以得到一般递

26、推公式如下：61精选课件由k=1开始递推，递推到k=p，依次得到各阶模型的参数， 递推公式提供了一种确定模型阶数的实验方法，如模型的阶数不知道，由低阶开始递推，当递推到M阶时，预测误差满足允许的值，停止递推，选AR模型的阶数为M。 这种方法递推效率高，当阶数变化时，无需从头计算。但需要预先估计出信号自相关函数，当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂；如数据很长，估计自相关函数较准确。计算流程图如下AR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算：62精选课件图 利用Levinson-Durbin递推法计算功率谱的流程图63精选课件1. 协方差法4.5

27、.2协方差法与修正协方差法 与自相关法一样，仍然利用使预测误差功率最小的方法求模型参数。但由观测数据求预测误差功率的公式不同，为 与自相关法相比，求预测误差功率的求和限不同。该公式中使用的观测数据均已知，不需要在数据两端补充零点，因此比较自相关法去掉了加窗处理的不合理假设。64精选课件为求模型参数仍然应用复梯度法使上式达到最小，公式如下式中白噪声的方差为称为协方差函数65精选课件 ，但不具有Topelitz性质( )上式由 组成的 维矩阵是对称的，即所以协方差方程的解不能借助于自相关法，可用乔里斯基(Choleskey)分解法求解上述方程组。用协方差法估计出的极点不能保证在单位圆内。66精选课

28、件 和前向预测误差 用下式表示： 修正协方差方法使用前向和后向观测误差平均值最小的方法，估计AR模型参数，进而估计信号的功率谱。 该方法最初由Nuttall在1976年提出，称为前向后向法；同年，Ulrych和Ctayton也独立提出，称之为最小二乘法，因此该法也称为前后向线性预测误差功率最小的最小二乘法。2. 修正协方差法下面对前向和后向线性预测误差滤波器进行分析。前向线性预测误差滤波器 为分析简单，假设信号是实平稳随机信号。前向线性预测误差滤波器由 预测 。其估计值67精选课件由于假设信号是实的，式中预测误差和系数均是实的。对上式进行Z变换，得到令则 称为前向预测误差滤波器的系统函数68精

29、选课件将 代入，得用均方误差最小准则求前向预测误差滤波器的最佳系数图 前向预测误差滤波器前向预测误差与用于预测的数据正交，正交原理。69精选课件式中 。已知自相关函数，可解出 和 。前向预测误差滤波器的最佳系数 和信号的自相关函数之间的关系为Yule-Walker方程，如下用矩阵表示为70精选课件 进行预测，为此，将上式改为后向线性预测误差滤波器利用 数据预测 ，则称为后向线性预测。其估计值为 一般前向、后向预测用同一数据进行，即利用这样，前向预测是由 预测 ，后向预测是由 预测预测 ，这两种预测数据之间的关系如图。71精选课件后向预测误差 表示的是信号在 时刻的预测误差，图 前后向预测数据之

30、间的关系同样，利用最小均方准则，可得关于后向预测时的正交原理及Yule-Walker方程72精选课件式中， 是后向预测误差的最小误差功率。利用Toeplitz性质，可得以下重要关系：上式表明，前、后向预测的最小误差功率相等，系数也相等(如果是复数，则是共轭关系)。73精选课件式中，当 时， ，上式可写成由 ， ，得图 后向预测误差滤波器后向预测误差滤波器的系数排序是前向预测误差滤波器系数排序的逆转排列。74精选课件对 进行Z变换，得后向预测误差滤波器的系统函数为前、后向预测误差滤波器的系统函数之间关系为75精选课件前面已推导出信号的前向和后向预测，分别为式中 是AR模型的参数。在这两种情况下，

31、最小预测误差功率都是白噪声方差 。修正协方差使用估计的前向和后向预测误差功率的平均值为最小的方法来估计AR参数。即76精选课件式中前向和后向预测误差功率 、 分别用下式表示和协方差法一样，仅对那些用到观测数据样本的预测误差求和。77精选课件经过简化，得到令将上式写成矩阵形式78精选课件白噪声的方差估计值为协方差法和修正协方差法相同，除了 的定义不同。修正协方差法还可以得到高分辨率的统计稳定的谱估计值。对于AR谱估计的谱峰由于附加观测噪声而偏离真实频率位置的现象很少。这种方法不能保证 ，因而不能保证滤波器的稳定性。79精选课件例：已知信号的四个观测数据为 ，分别用自相关法和协方差法估计AR(1)

32、模型参数。解：(1)自相关法80精选课件(2)协方差法81精选课件Levinson-Durbin递推法需要先由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型的参数。4.5.3伯格(Burg)递推法Levinson-Durbin递推法重写如下82精选课件 伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差和后向预测误差，并按照使它们最小的原则求出 ，从而实现不用估计自相关函数，直接用观测数据得出结果。首先由前向、后向预测误差滤波器推导预测误差格型滤波器前向预测误差公式重写如下将系数的递推公式 代入，得83精选课件 ，方括号部分可以写成 式右边最后一项中，因为 时，将上式与 对比，方程右边前两项是p-1阶前向预测误差，即84精选课件将上式与 式对比，该部分就是n-1时刻p-1阶的后向预测误差，即这样，得到前向预测误差的递推公式为类似地，得到后向预测误差的递推公式为格型滤波器第p节以及整个预测误差格型滤波器的结构如图。85精选课件图 全零点格型滤波器对于p=0的情况， 。86精选课件 设信号 观测数据区间为： ，前向、后