数字信号处理_第四章

1、第四章第四章 功功 率率 谱谱 估估 计计 4.1 引言引言 4.2 经典谱估计经典谱估计 4.3 现代谱估计中的参数建模现代谱估计中的参数建模 4.4 AR模型谱估计的性质模型谱估计的性质 4.5 AR谱估计的方法谱估计的方法 4.6 最大熵谱估计与最大似然谱估计最大熵谱估计与最大似然谱估计 4.7 特征分解法谱估计特征分解法谱估计 NO.24.1 引引 言言 对信号和系统进行分析研究、处理有两类方法：对信号和系统进行分析研究、处理有两类方法：一类是在时域进行，我们前面学习的维纳滤波、卡一类是在时域进行，我们前面学习的维纳滤波、卡尔曼滤波和自适应滤波都属于这种方法；本章则是尔曼滤波和自适应滤

2、波都属于这种方法；本章则是在频率域进行研究的一类方法。这两类方法都是信在频率域进行研究的一类方法。这两类方法都是信号处理的重要方法。对确定性信号傅里叶变换是在号处理的重要方法。对确定性信号傅里叶变换是在频率域分析研究的理论基础，但对于随机信号，其频率域分析研究的理论基础，但对于随机信号，其傅里叶变换不存在，因此转向研究它的功率谱。傅里叶变换不存在，因此转向研究它的功率谱。谱，就是信号的某些特征在频域随频率的分布。谱，就是信号的某些特征在频域随频率的分布。功率谱反映了随机信号功率的分布特性，有着很功率谱反映了随机信号功率的分布特性，有着很广泛的应用：在雷达信号处理中，回波信号的频广泛的应用：在雷

3、达信号处理中，回波信号的频谱提供了运动目标的位置、强度、速度等信息；谱提供了运动目标的位置、强度、速度等信息；在声纳系统中，为了寻找舰艇或潜艇也要对混有在声纳系统中，为了寻找舰艇或潜艇也要对混有噪声的信号进行谱分析；在语音处理中，谱分析噪声的信号进行谱分析；在语音处理中，谱分析用来探测语音语调共振峰；在电子战中，还利用用来探测语音语调共振峰；在电子战中，还利用频谱对目标进行分类。频谱对目标进行分类。mjmxxjxxmrPe )()e ( 按照按照Weiner-Khintchine定理，信号的功率谱定理，信号的功率谱和其自相关函数服从傅里叶变换关系：和其自相关函数服从傅里叶变换关系： 实际所能得

4、到的随机信号的长度总是有限的，实际所能得到的随机信号的长度总是有限的，用有限长度的信号计算得到的功率谱只是真实功率用有限长度的信号计算得到的功率谱只是真实功率谱的一个估计。功率谱估计分为谱的一个估计。功率谱估计分为经典谱估计经典谱估计和和现代现代谱估计谱估计两大类。两大类。dPmrnjjxxxxee)(21)()()()(*mnxnxEmrxx 上式称做功率谱的定义，对于平稳随机信号，服从上式称做功率谱的定义，对于平稳随机信号，服从各态历经定理，集合平均可以用时间平均代替，由上各态历经定理，集合平均可以用时间平均代替，由上式还可以推出功率谱的另一个定义：式还可以推出功率谱的另一个定义：)()(

5、121lim)(*mnxnxNmrNNnNxx将计算自相关函数中的集合平均用时间平均代替：将计算自相关函数中的集合平均用时间平均代替：mjNNnNmjxxmnxnxNPee)()(121lim)(*代入功率谱定义式，得代入功率谱定义式，得 mmnjnjNNnNmnxnxN)(*)()(121limee令令l=n+m, ,则则 ljlNNnnjNxxlxenxNPeej)()(121lim)(*mmnjnjNNnNjxxmnxnxNeP)(*)()(121lim)(ee2121jnNNnNnxNe)(lim 上式中上式中x(n)是观测数据，是观测数据，Pxx(ej)是随机变量，必是随机变量，必须

6、对须对Pxx(ej)取统计平均值，得到取统计平均值，得到 2j)(121lim)e (NNnnjNxxenxNEP该式被认为是功率谱的另一定义，周期图法谱估计该式被认为是功率谱的另一定义，周期图法谱估计 Weiner-khintchine定理表明功率谱是无限多个自定理表明功率谱是无限多个自相关函数的函数，但观测数据只有有限个相关函数的函数，但观测数据只有有限个, ,只能得到有只能得到有限个自相关函数。按照上面的定义公式求功率谱限个自相关函数。按照上面的定义公式求功率谱, ,也需也需要无限多个观测数据。因此根据有限个样本数据，计要无限多个观测数据。因此根据有限个样本数据，计算随机信号的功率谱，是

7、一个功率谱的估计问题算随机信号的功率谱，是一个功率谱的估计问题 现代谱估计是以信号模型为基础，下图表示现代谱估计是以信号模型为基础，下图表示x(n)的信号模型，输入白噪声的信号模型，输入白噪声w(n)的均值为的均值为0，方差为，方差为2w，x(n)的功率谱为：的功率谱为： 2j2j| )e (|)e (HPwxxH(z)w(n)x(n) 如果由观测数据能够估计出信号模型的参数，如果由观测数据能够估计出信号模型的参数，则信号的功率谱可以按照上式计算出来，这样则信号的功率谱可以按照上式计算出来，这样, ,估计估计功率谱的问题就变成了估计信号模型参数的问题。功率谱的问题就变成了估计信号模型参数的问题

8、。信号模型有很多种，如信号模型有很多种，如AR模型、模型、MA模型等等，针模型等等，针对不同的情况，需要选择不同的模型。现代谱估计对不同的情况，需要选择不同的模型。现代谱估计的质量比经典谱估计的质量有很大的提高。但遗憾的质量比经典谱估计的质量有很大的提高。但遗憾的是，尚无任何理论能指导选择一个合适的模型，的是，尚无任何理论能指导选择一个合适的模型，只能根据功率谱的一些先验知识，或者说一些重要只能根据功率谱的一些先验知识，或者说一些重要的谱特性，来选择模型。的谱特性，来选择模型。4.2 经经 典典 谱谱 估估 计计 BT法是先估计自相关函数法是先估计自相关函数, , 然后进行傅里叶变换然后进行傅

9、里叶变换得到功率谱。设对随机信号得到功率谱。设对随机信号x(n)，只观测到一段样本，只观测到一段样本数据，数据，n=0, 1, 2, , N-1。根据这一段样本数据估计自。根据这一段样本数据估计自相关函数相关函数, , 有两种估计方法，即有偏自相关函数估计有两种估计方法，即有偏自相关函数估计和无偏自相关函数估计。有偏自相关函数估计的误差和无偏自相关函数估计。有偏自相关函数估计的误差相对较小，这种估计是一种渐近一致估计：相对较小，这种估计是一种渐近一致估计： 1|0*)()(1)(mNnxxmnxnxNmr4.2.1 BT法法对上式进行傅里叶变换，得到对上式进行傅里叶变换，得到BT法的功率估计法

10、的功率估计: :mnjxxjmrP-BTe )()e ( 为了减少谱估计的方差，经常用窗函数为了减少谱估计的方差，经常用窗函数w(m)对自对自相关函数进行加权，此时谱估计公式为相关函数进行加权，此时谱估计公式为 1)1(-)()()(MMmnjxxjmwmrPeeBT式中式中 0)()(mwmw-(M-1)m(M-1) 其它其它 为了采用为了采用FFT计算傅里叶变换，必须将求和域计算傅里叶变换，必须将求和域(-M+1, M-1)移到移到( (0L-1)，功率谱的计算公式为：功率谱的计算公式为： 10-)()(LmmxxmSPjjBTee 上式也被称为加权协方差谱估计。它要求加窗后上式也被称为加

11、权协方差谱估计。它要求加窗后的功率谱仍是非负的，这样窗函数的功率谱仍是非负的，这样窗函数w(m)的选择必须满的选择必须满足一个原则，即它的傅里叶变换必须是非负的，例如足一个原则，即它的傅里叶变换必须是非负的，例如巴特利特窗就满足这一条件。巴特利特窗就满足这一条件。 k=0, 1, 2, , L-1 )()(mSFFTkPxxBT )()(0)()()(LmwLmrmwmrmSxxxxxx 0mM-1 MmL-M L-M+1mL-1 按照有偏自相关函数公式估计自相关函数，已经证明按照有偏自相关函数公式估计自相关函数，已经证明这是渐近一致估计，但经过傅里叶变换后得到功率谱这是渐近一致估计，但经过傅

12、里叶变换后得到功率谱的估计，功率谱估计却不一定仍是渐近一致估计，可的估计，功率谱估计却不一定仍是渐近一致估计，可以证明它是非一致估计，是一种不好的估计方法。以证明它是非一致估计，是一种不好的估计方法。BT法中用有偏自相关函数进行估计时，它和用周期图法法中用有偏自相关函数进行估计时，它和用周期图法估计功率谱是等价的，因此估计功率谱是等价的，因此BT法的估计质量和周期图法的估计质量和周期图法的估计质量是一样的。法的估计质量是一样的。 将功率谱的另一定义式重写如下将功率谱的另一定义式重写如下: : 2jje )(121lim)e (NNnnNxxnxNEP如果忽略上式中求统计平均的运算，假设观测数据

13、如果忽略上式中求统计平均的运算，假设观测数据为为: :x(n) 0nN-1，便得到周期图法的定义，便得到周期图法的定义: : 210j -je )(1)e (NnnxxnxNP4.2.2 周期图法周期图法用周期图法计算功率谱框图用周期图法计算功率谱框图 观 测 数 据 x(n)FFT取 模的 平 方1/NPxx(ej)由周期图法功率谱估计公式推导它与由周期图法功率谱估计公式推导它与BT法的等价关系法的等价关系210)(1)(NnnjjxxnxNPee令令 m = k-n, 即即k = m+nmjmNnNNmxxnmxnxNPee10*1)1(j)()(1)(1.周期图与周期图与BT法的等价关系

14、法的等价关系10*10)()(1NnnjNkkjnxkxNee1010)(*)()(1NkNnnkjnxkxNe方括号中的部分是有偏自相关函数的计算公式：方括号中的部分是有偏自相关函数的计算公式：mjxxNNmxxmrP)e()e (1)1(j)()(1)(*|10nmxnxNmrmNnxx)e ()e (BTjjxxPP利用有偏自相关函数的利用有偏自相关函数的BT法和周期图法是等价的关系法和周期图法是等价的关系已知自相关函数的估计值已知自相关函数的估计值 ，m= -(N-1),-1,0,1, , N-1，按照，按照BT法求功率谱的统计平均值：法求功率谱的统计平均值：)(mrxxmxxNNmx

15、xmrEPEj -1)1(je)()e (有偏自相关函数的统计平均值在第一章中已确定，将有偏自相关函数的统计平均值在第一章中已确定，将结果代入上式结果代入上式, ,得到得到 mxxNNmxxmrNmNPEj -1)1(j)e(|)e (2.周期图法谱估计质量分析周期图法谱估计质量分析1) )周期图的偏移周期图的偏移mmxxBxxmrmwPEj -j)e()()e (式中式中 其它0|)(NmNmNmwB两序列乘积的傅里叶变换，其频域服从卷积关系两序列乘积的傅里叶变换，其频域服从卷积关系)()e (mrFTPxxjxx)()()e (jxxjBjxxePeWPE2)2/sin()2/sin(1)

16、()e (NNmwFTWBjBWB(ej)称为三角窗的谱函数。上式表明，周期图的统称为三角窗的谱函数。上式表明，周期图的统计平均值等于它的真值与三角窗函数频谱的卷积，因计平均值等于它的真值与三角窗函数频谱的卷积，因此周期图是有偏估计，但当此周期图是有偏估计，但当N时，时，wB(m)1，三，三角窗函数的频谱趋近于角窗函数的频谱趋近于函数，周期图的统计平均值函数，周期图的统计平均值趋于它的真值，因此周期图属于渐近无偏估计。趋于它的真值，因此周期图属于渐近无偏估计。 周期图的方差的精确表达式很繁冗，为分析简单：假周期图的方差的精确表达式很繁冗，为分析简单：假设设x(n)是实的零均值的正态白噪声信号，

17、方差是是实的零均值的正态白噪声信号，方差是x2，即功率谱是常数即功率谱是常数x2 ，其周期图用，其周期图用IN()表示，表示，N表示观表示观测数据的长度。按照周期图的定义测数据的长度。按照周期图的定义: :2| )(|1)(jeXNIN推导得周期图的方差：推导得周期图的方差：2) )周期图的方差周期图的方差24)sin()sin(1)(varNNIxNnNnNkknxkxNj-jee1010)()(1 当当N趋于无限大时，周期图的方差并不趋于趋于无限大时，周期图的方差并不趋于0，而，而趋于功率谱真值的平方，即趋于功率谱真值的平方，即 4)(varxNNI 所以无论怎样选择所以无论怎样选择N，周

18、期图的方差总是和，周期图的方差总是和4x同同一个数量级。信号功率谱的真值是一个数量级。信号功率谱的真值是2x，这说明周期，这说明周期图的方差很大。用这种方法估计的功率谱在图的方差很大。用这种方法估计的功率谱在2x附近附近起伏很大，所以周期图是非一致估计，是一种很差起伏很大，所以周期图是非一致估计，是一种很差的功率谱估计方法。的功率谱估计方法。 图图4.2.2 白噪声的周期图白噪声的周期图 可以看到，随着可以看到，随着N的增大，功率谱曲线并没有趋的增大，功率谱曲线并没有趋向真值，而是起伏变的越来越剧烈。向真值，而是起伏变的越来越剧烈。 周期图法估计功率谱不是一致估计，均方误差周期图法估计功率谱不

19、是一致估计，均方误差很大，使估计出的功率谱不可靠。其频率分辨力低很大，使估计出的功率谱不可靠。其频率分辨力低是根本缺点，原因在于观测数据只有一段。由于是根本缺点，原因在于观测数据只有一段。由于BT法和周期图法具有等效的功率谱估计质量，因此法和周期图法具有等效的功率谱估计质量，因此BT法也不是一致估计，分辨率低，估计误差大。法也不是一致估计，分辨率低，估计误差大。4.2.3 经典谱估计方法的改进经典谱估计方法的改进 基本思想基本思想：对随机变量进行观测，得到：对随机变量进行观测，得到L组独立的记组独立的记录数据，用每一组数据求其周期图，然后将录数据，用每一组数据求其周期图，然后将L个周期个周期图

20、加起来求平均。这样得到的周期图，其方差将是用图加起来求平均。这样得到的周期图，其方差将是用一组数据得到的周期图的方差的一组数据得到的周期图的方差的1/L。210)(1)(MnnjiienxMI1.平均周期图法平均周期图法 假设随机信号假设随机信号x(n)的观测区间为：的观测区间为：0nM-1，共，共进行了进行了L次独立观测，得到次独立观测，得到L组记录数据，某一组记录组记录数据，某一组记录数据用数据用xi(n), i=1,2,L表示，第表示，第i组的周期图为：组的周期图为： 将得到的将得到的L个周期图进行平均，作为信号个周期图进行平均，作为信号x(n)的功的功率谱估计：率谱估计：LiixxIL

21、P1j)(1)e (为了分析偏移，对上式求统计平均为了分析偏移，对上式求统计平均: : 其中，窗函数的频谱为其中，窗函数的频谱为 2j)2/sin()2/sin(1)()e (MMmwFTWBB)()(1)(1IEIELPELiixxje)()()(jxxjBePeWIE表明平均周期图仍然是有偏估计，偏移和每一段的数表明平均周期图仍然是有偏估计，偏移和每一段的数据个数据个数M有关；由于有关；由于M N, ,平均周期图的偏移比周期平均周期图的偏移比周期图的偏移大，表现在三角谱窗主瓣的宽度比周期图主图的偏移大，表现在三角谱窗主瓣的宽度比周期图主瓣的宽度宽。所以其分辨率更加降低，因此也可以说，瓣的宽

22、度宽。所以其分辨率更加降低，因此也可以说，偏移的大小反映分辨率的低与高。偏移的大小反映分辨率的低与高。 对平均周期图求方差，由于是对平均周期图求方差，由于是L次独立观测，次独立观测，L个个周期图相互独立，因此平均周期图的方差为周期图相互独立，因此平均周期图的方差为 )(var1)e (varijxxILP即平均周期图的估计方差是周期图的方差的即平均周期图的估计方差是周期图的方差的1/L。即：。即：是以分辨率的降低为代价换取了估计方差的减少是以分辨率的降低为代价换取了估计方差的减少平均周期图法平均周期图法 Pxx(ej)/032101230321012303210123Pxx(ej)/Pxx(e

23、j)/(a)(b)(c)N 256， L 2N 256， L 4N 256， L 8 这种方法是用一适当的窗函数这种方法是用一适当的窗函数W(ej)与周期图进与周期图进行卷积，来达到使周期图平滑的目的。行卷积，来达到使周期图平滑的目的。 dee)()(21)()(jNlwxxWIP式中式中 mxxNNmNmrIj -1)1(e )()( 是有偏自相关函数是有偏自相关函数 )(mrxxde )e (21)(jjnWnw-(M-1)nM-1 2.窗函数法窗函数法mxxMMmjxxemwmrPje)()()(1)1( 周期图的窗函数法和前面提到的周期图的窗函数法和前面提到的BT法的加权协法的加权协方

24、差谱估计是类似的。窗函数法中，周期图和窗函方差谱估计是类似的。窗函数法中，周期图和窗函数的频谱卷积得到功率谱，等效于在频域对周期图数的频谱卷积得到功率谱，等效于在频域对周期图进行修正，使周期图通过一个线性系统，滤除掉周进行修正，使周期图通过一个线性系统，滤除掉周期图中的快变成分，谱窗函数需具有低通特性。期图中的快变成分，谱窗函数需具有低通特性。对上式求统计平均对上式求统计平均, , 得到得到 mxxMMmxxmwmrEPEj -1)1(je )()( )e (有偏自相关函数的期望等效于真值与三角函数的乘积有偏自相关函数的期望等效于真值与三角函数的乘积mBxxMMmxxmwmwmrePEj -1

25、)1(je )()()()(NmNmmwB|1)( 上式表明上式表明, , 周期图的窗函数法仍然是有偏估计周期图的窗函数法仍然是有偏估计, , 其其偏移和偏移和wB(m)、w(m)两个窗函数有关，如果两个窗函数有关，如果w(m)窗的窗的宽度比较窄，当宽度比较窄，当M比比N小得多时，小得多时，|M| q时，时，上式是一个线性方程上式是一个线性方程)()(0*mlhlhmqlBqmlmrlhmrplxxAxx1)()()(go上式共有上式共有p个方程。可以用该方程首先计算出个方程。可以用该方程首先计算出AR部分部分的的p个系数个系数hA(i), i = 1, 2, 3, p ; 然后代入关系式，设

26、然后代入关系式，设法求出法求出MA部分的系数。部分的系数。 )()2() 1()()2() 1 ()()2() 1()2()() 1() 1() 1()(pqrqrqrphhhqrpqrpqrpqrqrqrpqrqrqrxxxxxxAAAxxxxxxxxxxxxxxxxxx将上式以矩阵形式表示：将上式以矩阵形式表示：qmlmrlhmrplxxAxx1)()()(AR模型的系统函数为：模型的系统函数为：H(z)=1/A(z)，相当于，相当于ARMA模模型中型中B(z)=1的情况，这样在公式中的情况，这样在公式中 )()(nnhBPlwxxAPlxxAxxlrlhlmrlhmr121)()()()

27、()(m1 m=0=0 2.AR模型的系数和信号自相关函数之间的关系模型的系数和信号自相关函数之间的关系plxxAplmqlBwxxAxxlmrlhmlhlhlmrlhmr110*2)()()()()()()(m=0, 1, 2, 3, , q mq+1 h(0) = ?整个关系式用矩阵表示为整个关系式用矩阵表示为 00)() 1 (1)0() 1()() 1()0() 1 ()() 1()0(2wAAxxxxxxxxxxxxxxxxxxphhrprprprrrprrr可以用模型参数表示为：可以用模型参数表示为： 001)0() 1()() 1()0() 1 ()() 1()0(21wpxxx

28、xxxxxxxxxxxxxxxaarprprprrrprrr 为自相关矩阵，它满足为自相关矩阵，它满足 ，且沿任一对，且沿任一对角线的元素相等，它也是正定矩阵角线的元素相等，它也是正定矩阵 上面推导出的关系式或模型系数关系式确定了上面推导出的关系式或模型系数关系式确定了AR模型参数模型参数( (包括模型输入噪声方差包括模型输入噪声方差) )和信号自相关和信号自相关函数之间的关系。这是一个线性方程，如果能够由信函数之间的关系。这是一个线性方程，如果能够由信号的观测数据求出信号的自相关函数，可以按照该公号的观测数据求出信号的自相关函数，可以按照该公式，通过解线性方程得到模型参数式，通过解线性方程得

29、到模型参数. .)1( pxxRHxxxxRRMA模型的系统函数模型的系统函数H(z)=B(z)，相当于，相当于ARMA模型中模型中A(z)=1, hA(n)=(n)的情况，此时的情况，此时h(n)=hB(n)，由，由ARMA模型可得到模型可得到MA模型的系数和自相关函数的关系模型的系数和自相关函数的关系0)()()(0*2mqlBBwxxmlhlhmrm=0, 1, , q mq+1 MA模型的参数和自相关函数之间也是非线性关系模型的参数和自相关函数之间也是非线性关系3. MA模型的系数和信号自相关函数之间的关系模型的系数和信号自相关函数之间的关系plxxAplmqlBwxxAxxlmrlh

30、mlhlhlmrlhmr110*2)()()()()()()(m=0, 1, 2, 3, , q mq+1 三种信号模型的参数和信号自相关函数之间的关系。三种信号模型的参数和信号自相关函数之间的关系。这些关系式提供了一种估计功率谱的方法，即首先根这些关系式提供了一种估计功率谱的方法，即首先根据信号观测数据估计信号自相关函数，然后按照所选据信号观测数据估计信号自相关函数，然后按照所选择的信号模型，解上面相应的方程，求出模型参数，择的信号模型，解上面相应的方程，求出模型参数，最后按照下式求出信号的功率谱最后按照下式求出信号的功率谱: : 21022j2je1e| )e (|)e (Piiiqiii

31、wwxxabHPH(z)w(n)x(n)参数模型参数模型信号的功率谱信号的功率谱传统的功率谱估计方法无论采取哪一种改进方法，总传统的功率谱估计方法无论采取哪一种改进方法，总是以减少分辨率为代价来换取估计方差的减少，其最是以减少分辨率为代价来换取估计方差的减少，其最大缺点是：频谱的分辨率低大缺点是：频谱的分辨率低MA,ARMA模型：自相关函数和模型参数之间呈模型：自相关函数和模型参数之间呈非非线性的线性的关系，增加了估计功率谱的困难关系，增加了估计功率谱的困难21211)(PiiiwxxaPeej用用n时刻时刻之前的之前的p个数据为：个数据为：x(n-1)，x(n-2), ,x(n-p)，预测预

32、测n时刻的数据时刻的数据x(n)：pipiinxanx1)()( ia预测预测误差误差为：为：)( )()(nxnxnepipiinxanx1)()(4.4.1 AR模型的线性预测模型的线性预测4.4 AR模型谱估计的性质模型谱估计的性质其中，其中， 为预测值，为预测值， 为预测系数为预测系数)( nxpiae(n)表示线性一步预测误差表示线性一步预测误差api(i=1, 2, 3, ，p)表示预测器的系数，它和线性预表示预测器的系数，它和线性预测器单位脉冲响应测器单位脉冲响应h(n)差一负号差一负号pnanhpn, 3 , 2 , 1)(对差分方程进行对差分方程进行z变换变换 piipiza

33、zXzXzE1)()()(pipiinxanxne1)()()(输入是观测数据，输出的是预测误差输入是观测数据，输出的是预测误差令令He(z)=E(z)/X(z)，得到预测误差滤波器的系统函数，得到预测误差滤波器的系统函数 piipiezazH11)(称称He(z)为一步为一步线性线性预测误差滤波器，其作用是将信预测误差滤波器，其作用是将信号号x(n)转换成预测误差转换成预测误差e(n)。一般认为。一般认为e(n)具有白噪声具有白噪声的性质，因此的性质，因此He(z)也称为白化滤波器。也称为白化滤波器。 预测误差滤波器预测误差滤波器/ /白化滤波器白化滤波器 He(z)x(n)e(n)x(n)

34、e(n)z1z1z1ap1ap2ap, p1ap, pAR模型的系统函数为模型的系统函数为: : )(111)(1zAzazHPiii当当api=ai(i=1, 2, 3, ,p)时，时，He(z)和和H(z)互为逆滤波器，互为逆滤波器， He(z) = 1/H(z) = A(z)白化滤波器的系统函数白化滤波器的系统函数: :piipiezazH11)(预测误差滤波器预测误差滤波器AR模型模型 当当AR模型的阶数与线性预测误差滤波器的阶数相模型的阶数与线性预测误差滤波器的阶数相同时，二者互为逆滤波器。同时，二者互为逆滤波器。AR模型和预测误差滤波器模型和预测误差滤波器的级联相当于一个全通网络。

35、的级联相当于一个全通网络。w(n)是是AR模型输入白噪模型输入白噪声，故预测误差声，故预测误差e(n)具有白噪声的性质，所以，预测具有白噪声的性质，所以，预测误差滤波器也叫白化滤波器。并且，最小预测误差功误差滤波器也叫白化滤波器。并且，最小预测误差功率就是激励源白噪声的方差率就是激励源白噪声的方差w2。 由于由于AR模型具有这种特性，因而模型具有这种特性，因而AR模型法也称模型法也称为线性预测为线性预测AR模型法。模型法。 AR模型必须是因果稳定的，即极点均在单位圆模型必须是因果稳定的，即极点均在单位圆内，才能保证信号内，才能保证信号x(n)是平稳随机信号，由于是平稳随机信号，由于AR模型模型

36、H(z)和预测误差滤波器和预测误差滤波器He(z)互为逆滤波器，互为逆滤波器，所以所以He(z)应为最小相位系统。应为最小相位系统。4.4.2 预测误差滤波器的最小相位特性预测误差滤波器的最小相位特性 自相关函数和自相关函数和AR模型的模型参数之间的关系服从模型的模型参数之间的关系服从尤勒尤勒- -沃克沃克( (Yule-Walker) )方程：方程： plwxxAplxxAxxlrlhlmrlhmr121)()()()()(m1 m=0 上式中，对于上式中，对于m1的情况，公式本身就是一个递推方的情况，公式本身就是一个递推方程，如果已由观测数据计算出程，如果已由观测数据计算出p+1个自相关函

37、数，个自相关函数，用用 , ,m = 0,1,2, , ,p表示，对于表示，对于mp的情况的情况, , 可以用该公式外推得到：可以用该公式外推得到： )(mrxx4.4.3 AR模型隐含自相关函数延拓特性模型隐含自相关函数延拓特性PlxxAxxxxlmrlhmrmr1)( )()()(0mp m p 式中，系数式中，系数hA(l)是是用前用前p+1个自相关函数求出的参数个自相关函数求出的参数因此因此AR模型隐含着自相关函数外推的特性。在经典模型隐含着自相关函数外推的特性。在经典谱估计谱估计BT法中，自相关函数只能限于由观测数据计算法中，自相关函数只能限于由观测数据计算出的有限个自相关函数，其它

38、的认为是出的有限个自相关函数，其它的认为是0，造成了谱，造成了谱估计分辨率低、模糊。正是估计分辨率低、模糊。正是AR模型具有自相关函数模型具有自相关函数外推的特性，使外推的特性，使AR估计的谱具有高分辨率的优点。估计的谱具有高分辨率的优点。 信号频域的分辨率与时域长度的关系：信号频域的分辨率与时域长度的关系：TF14.5 AR谱估计的方法谱估计的方法 自相关法的出发点是选择自相关法的出发点是选择AR模型的参数使预测模型的参数使预测误差功率最小，预测误差功率为误差功率最小，预测误差功率为 npipininxanxNneN122| )()(|1| )(|1假设信号假设信号x(n)的数据区在的数据区

39、在0nN-1范围，有范围，有p个预测系个预测系数，数，N个数据经过冲激响应为个数据经过冲激响应为api(i=0,1, 2, , p)的滤波的滤波器，输出预测误差器，输出预测误差e(n)的长度为的长度为N+P4.5.1 自相关法自相关法列文森列文森( (Levenson) )递推法递推法由信号的观测数据，估计信号的功率谱由信号的观测数据，估计信号的功率谱102| )(|1pNnneN显然显然, ,e(n)的长度长于数据的长度，上式中数据的长度长于数据的长度，上式中数据x(n)的的两端需补充零点，这相当于无穷长的信号经过加窗处两端需补充零点，这相当于无穷长的信号经过加窗处理，得到长度为理，得到长度

40、为N的数据。用上式对系数的数据。用上式对系数api求微分的求微分的方法使预测误差功率最小，得到方法使预测误差功率最小，得到 )()2() 1 ()0()2() 1()2()0() 1 () 1() 1()0(21prrraaarprprprrrprrrxxxxxxppppxxxxxxxxxxxxxxxxxx上面的矩阵就是上面的矩阵就是Yule-Walker方程方程式中的自相关函数采用有偏自相关估计式中的自相关函数采用有偏自相关估计)()()(1)(*10mrmnxnxNmrxxmNnxxm=0, 1, 2, , p m=-p, -p+2, , -1 因此自相关法是基于解因此自相关法是基于解Yu

41、le-Walker方程的一种方法。方程的一种方法。首先由信号的观测数据估计出自相关函数，再解该方首先由信号的观测数据估计出自相关函数，再解该方程，得到模型参数，便可求出信号的功率谱。因此该程，得到模型参数，便可求出信号的功率谱。因此该方法也称为方法也称为Yule-Walker法。但是直接解该方程，需要法。但是直接解该方程，需要计算逆矩阵，很不方便。利用计算逆矩阵，很不方便。利用Yule-Walker方程中自相方程中自相关矩阵的性质，可以导出关矩阵的性质，可以导出Levenson-Durbin递推法，这递推法，这是一种高效的解方程方法。简称为列文森递推法。是一种高效的解方程方法。简称为列文森递推

42、法。2111, 1/)()(kklxxlkxxkklkrakra*, 1, 1ikkkkikkiaaaai=1, 2, 3, , k-1 2122)|1 (kkkka由由k=1开始递推，递推到开始递推，递推到k=p，依次得到，依次得到a11,21，a21,a22, 22，ap1,ap2,app,2p。AR模型的系数以模型的系数以及模型输入白噪声方差求出后，信号功率谱为：及模型输入白噪声方差求出后，信号功率谱为： 2122211| )(|)(PiiiwwxxaHPj-jjeeeu上式表明：上式表明： ，说明随着阶数增加，说明随着阶数增加， 预测误差功率将减少或者不变，为此要求预测误差功率将减少或

43、者不变，为此要求|akk|1,akk称称为反射系数。为反射系数。u递推公式提供了一种确定模型阶数的方法，如模型递推公式提供了一种确定模型阶数的方法，如模型的阶数未知，由低阶开始递推，当递推到的阶数未知，由低阶开始递推，当递推到M阶时，预阶时，预测误差满足要求，则停止递推，选测误差满足要求，则停止递推，选AR模型的阶数为模型的阶数为M。递推法效率高，当阶数变化时，无需重新计算递推法效率高，当阶数变化时，无需重新计算u如果已知如果已知N个观测数据个观测数据（x(n),0nN-1），利用列文），利用列文森递推法计算功率谱的计算流程图如下图所示森递推法计算功率谱的计算流程图如下图所示2232221p利

44、用列文森递推法计算功率谱的流程图利用列文森递推法计算功率谱的流程图 输入： x(n),n0,1,2,3,N1,2计算：)0(1),0(/ ) 1 (21 , 1211 , 1xxxxxxrarra)(mrxxp=22122*, 1, 12111, 111, 2 , 1 , 0)()(ppppkppppkppkppixxipxxppapkaaaaiprapra22ppiipiwxxaP1j2je11)(e结束p=p+1NYk=1, 2, 3, , p-1 该方法和自相关法一样，仍利用预测误差功率最小方该方法和自相关法一样，仍利用预测误差功率最小方法求模型参数，但求预测误差功率的公式不同：法求模型

45、参数，但求预测误差功率的公式不同： 11212| )()(|1| )(|1NpnNpnpkpkknxanxpNnepN对比自相关法求预测误差功率的公式，不同的是求和对比自相关法求预测误差功率的公式，不同的是求和限不同。该公式使用的观测数据是已知的，不需要在限不同。该公式使用的观测数据是已知的，不需要在数据两端补充零点，因此与自相关法相比，去掉了加数据两端补充零点，因此与自相关法相比，去掉了加窗处理的不合理假设。窗处理的不合理假设。4.5.2 协方差法与修正协方差法协方差法与修正协方差法1. .协方差法协方差法)0 ,()0 , 2()0 , 1 (),()2 ,() 1 ,(), 2()2 ,

46、 2() 1 , 2(), 1 ()2 , 1 () 1 , 1 (21pcccaaappcpcpcpcccpcccxxxxxxppppxxxxxxxxxxxxxxxxxx1*)()(1),(NpnxxknxjnxpNkjc白噪声的方差为白噪声的方差为 1min2), 0()0 , 0(Npkxxpkxxwkcac仍然使用梯度最小的方法求模型参数：仍然使用梯度最小的方法求模型参数： 协方差函数协方差函数 观测数据观测数据x(n)(n=0, 1, 2, , N-1)，利用上面公式，利用上面公式可以求出模型的参数：可以求出模型的参数：api(i=1, 2, 3, , p); 2w。式。式中的协方差

47、函数中的协方差函数cxx(j, k)，有两个变量，因此也适合，有两个变量，因此也适合于非平稳随机信号。式中的协方差矩阵是埃尔米特于非平稳随机信号。式中的协方差矩阵是埃尔米特（Hermitian）矩阵，是半正定的。这种方法近似于）矩阵，是半正定的。这种方法近似于自相关法。一些实验结果说明它的分辨率优于自相自相关法。一些实验结果说明它的分辨率优于自相关法，另外对于纯正弦信号数据，可以有效地估计关法，另外对于纯正弦信号数据，可以有效地估计正弦信号的频率。正弦信号的频率。 修正协方差法使用前向和后向预测误差平均值修正协方差法使用前向和后向预测误差平均值最小的方法，估计最小的方法，估计AR模型的参数，从

48、而得到信号模型的参数，从而得到信号的功率谱。信号的前向和后向预测分别如下：的功率谱。信号的前向和后向预测分别如下： pkpkpkpkknxanxknxanx11)()( )()( 2.修正协方差法修正协方差法式中式中apk是是AR模型的参数模型的参数最小预测误差平均功率是模型输入白噪声的方差，即最小预测误差平均功率是模型输入白噪声的方差，即p=2w，前、后向预测误差平均功率为，前、后向预测误差平均功率为 )(5 . 0pbpep前向和后向预测误差功率前向和后向预测误差功率pe、pb分别用下式表示分别用下式表示2101211)()(1)()(1pNnpkpkpbNpnpkpkpeknxanxpN

49、knxanxpN和协方差法一样，上式仅对用到的观测数据的预测误和协方差法一样，上式仅对用到的观测数据的预测误差求和。为了使预测误差平均功率最小，求差求和。为了使预测误差平均功率最小，求p对对apk(k=1, 2, 3, , p)的微分的微分, ,或者用复梯度法求，得到或者用复梯度法求，得到 pllnxknxanxlnxknxanxpNapkpkpNnNpnpkpkplp, 3 , 2 , 10)()()()()()(111011化简并写成矩阵形式为：化简并写成矩阵形式为： )0 ,()0 , 2()0 , 1 (),()2 ,() 1 ,(), 2()2 , 2() 1 , 2(), 1 ()

50、2 , 1 () 1 , 1 (21pcccaaappcpcpcpcccpcccxxxxxxppppxxxxxxxxxxxxxxxxxx白噪声的方差估计值为白噪声的方差估计值为 pkxxpkxxpNnpkpkNpnpkpkpwkcacnxknxanxnxknxanxpN1101*1*1min,2), 0()0 , 0()()()()()()()(21例例: :已知信号的四个观察数据为已知信号的四个观察数据为x(n)=x(0), x(1), x(2), x(3)=2, 4, 1, 3，分别用自相关法和协方差法估计，分别用自相关法和协方差法估计AR(1)模型参数。模型参数。 011a210| )(

51、|1neNPNn1) )自相关法：自相关法： 解解402| )(|41nne2110| )()(|1pipiPNninxanxN09)3()41 (4)24(211111111aaaa) 1()()(11nxanxne0) 1()(21)()(21401140nxneanenenn5 . 011a2) )协方差法协方差法: : ) 1()()(| )(|3111312nxanxnenen11212| )()(|1| )(|1NpnNpnpkpkknxanxpNnepN311131110) 1()(32)()(32nnnxneanenea0)3()41 (4)24(2111111aaax(n)=

52、2, 4, 1, 3714. 011a4.5.3 伯格伯格( (Burg) )递推法递推法设信号设信号x(n)观测数据区间为：观测数据区间为：0nN-1，前向、后向，前向、后向预测误差功率分别用预测误差功率分别用p,e和和p,b表示，预测误差平均功表示，预测误差平均功率用率用p表示，公式分别为表示，公式分别为 12,| )(|1NpnfpepnepN12,| )(|1NpnbpbpnepN)(21,bpepp前向、后向预测误差递推公式如下：前向、后向预测误差递推公式如下： )() 1()() 1()()(1*111nekneneneknenefppbpbpbppfpfp将上式带入平均误差公式中

53、，得到将上式带入平均误差公式中，得到 121*1211| )() 1(| ) 1()(|)(21NpnfppfpbppfppnekneneknepN求预测误差平均功率求预测误差平均功率p最小时的反射系数最小时的反射系数kp， 0ppk)| ) 1(| )(|) 1()(2121211*11NpnbpfpNpnbpfppnenenenek上式就是利用伯格递推法求第上式就是利用伯格递推法求第p个反射系数的公式个反射系数的公式将伯格递推法求将伯格递推法求AR模型参数的递推公式总结如下：模型参数的递推公式总结如下： )| ) 1(| )(|) 1()(21, 1 , 0)()(1, 1 , 0)()(

54、)0(| )(|1)0(211211*11000102NpnbpfpNpnbpfppbfxxNnxxnenenenekNnnxneNnnxnernxNr2, 1,)() 1()(1, 2, 1)()()(1, 3 , 2 , 1)|1 (1*1111,*, 1, 112NppnnekneneNppnneknenekapiakaakfppbpbpbppfpfppppipppippippp利用伯格递推法求利用伯格递推法求AR模型参数的流程图如下图所示模型参数的流程图如下图所示 伯格递推法流程图伯格递推法流程图书书159页页 输 入 ： x(n),n 0,1,2,3,N 1阶 数 IP1)(1)()

55、()()(210000pnxNnxnenxneNnbfp = IP结 束1 pp12121*111)()()1()(2NpnbpfpbpNpnfppnenenenekpwpppk2121pppipppippikaakaa*, 1, 11, 2, 1),()1()(1, 2, 1),1()()(1*111NppnnekneneNppnneknenefppbpbpbpppfp2,;, 3 ,2, 1:wippia输出21j2je1)(epiipiwxxaPNYn = p N-1 p = 0 计算计算AR模型参数应用较多的为以下三种方法：模型参数应用较多的为以下三种方法：1、自相关法、自相关法 -（

56、Levinson法）法）2、协方差法、协方差法 / / 改进协方差法改进协方差法3、Burg法法自相关法（自相关法（Levinson法）法）计算简单，但由于假设数据计算简单，但由于假设数据外为零使得分辨率相对较差外为零使得分辨率相对较差协方差法协方差法与自相关法的主要区别是：求和范围不同。与自相关法的主要区别是：求和范围不同。数据段两端不需要添加任何零取样值，没有假设数据数据段两端不需要添加任何零取样值，没有假设数据段段N以外的数据等于零。以外的数据等于零。谱估计性能较好，但潜在不稳定因素谱估计性能较好，但潜在不稳定因素改进协方差法改进协方差法piNpnbpfpanenepN|min)()(1

57、122谱估计性能最好，但计算过于复杂谱估计性能最好，但计算过于复杂Burg法法不直接估计不直接估计AR参数，而是先估计反射系数，参数，而是先估计反射系数，然后利用然后利用Levinson算法由反射系数计算算法由反射系数计算AR参数参数计算反射系数时，使用前向、后向预测误差功率的平计算反射系数时，使用前向、后向预测误差功率的平均值最小准则。均值最小准则。ppNpnbpfpanene|min)()(122计算不太复杂，且给出了较好的谱估计质量计算不太复杂，且给出了较好的谱估计质量. .是较为通用的方法是较为通用的方法在在AR模型谱估计中，模型阶数的选择是一个关键问题。模型谱估计中，模型阶数的选择是

58、一个关键问题。一般模型的最好选择是先验未知的，实际中需预先选一般模型的最好选择是先验未知的，实际中需预先选定模型阶次。如果是纯定模型阶次。如果是纯p阶阶AR信号，选择模型阶次信号，选择模型阶次kp时，将产生对谱的平滑作用，降低谱的分辨率，如时，将产生对谱的平滑作用，降低谱的分辨率，如下图所示。图中，下图所示。图中，AR信号信号p=4，选择模型的阶次，选择模型的阶次k=2，产生的平滑作用使两个峰变成一个峰，分辨率明显降产生的平滑作用使两个峰变成一个峰，分辨率明显降低。如果选择低。如果选择kp，且假定观测的数据没有误差（没，且假定观测的数据没有误差（没有干扰），估计的参数应是：有干扰），估计的参数

59、应是： kppipiaaipki, 2, 10, 3 , 2 , 1,4.5.4 关于关于AR模型阶次的选择模型阶次的选择AR模型阶次太小时的平滑作用模型阶次太小时的平滑作用50.0030.0010.0010.0030.000.0000.400.50真实AR(4)的PSDAR(2)的PSD模型PSD/dB/因此因此, ,对于纯对于纯p阶阶AR信号，应选择阶次信号，应选择阶次kP。如果是白。如果是白噪声中的噪声中的AR信号信号( (观测数据有误差或者信号中含有白观测数据有误差或者信号中含有白噪声噪声) )，此时选择，此时选择ARMA模型合适，如选择了模型合适，如选择了AR模

60、型，模型，其阶次应加大，较低的阶次会使谱估计产生偏移，降其阶次应加大，较低的阶次会使谱估计产生偏移，降低分辨率。当然低分辨率。当然, , 这也和信噪比有关，信噪比愈低，这也和信噪比有关，信噪比愈低，平滑作用愈严重，愈需要高的阶次，因此信噪比低应平滑作用愈严重，愈需要高的阶次，因此信噪比低应选高的阶次。一般来说，阶次愈高，分辨率愈高；但选高的阶次。一般来说，阶次愈高，分辨率愈高；但阶次太高，会因参数过多而使估计误差加大，谱峰分阶次太高，会因参数过多而使估计误差加大，谱峰分裂，因此，对于白噪中的裂，因此，对于白噪中的AR信号，其阶次的选择应折信号，其阶次的选择应折衷考虑衷考虑4.6 最大熵谱估计与