数字信号处理-第4章

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

2、一类方法。 这两类方法都是信号处理的重要方法。对确定性信号傅里叶变换是在频率域分析研究的理论基础，但对于随机信号，其傅里叶变换并不存在，因此转向研究它的功率谱。按照Weiner-Khintchine定理，信号的功率谱和其自相关函数服从一对傅里叶变换关系， 公式如下 mjmxxjxxmrPe )()e (4.1.1) 第四章第四章 功功 率率 谱谱 估估 计计 dPmrnjjxxxxe )e (21)()()()(*mnxnxEmrxx(4.1.2) (4.1.3) （4.1.1）式被称做功率谱的定义，对于平稳随机信号，服从各态历经定理，集合平均可以用时间平均代替，由（4.1.1）式还可以推出功

3、率谱的另一个定义，推导如下： 将（4.1.3）式中的集合平均用时间平均代替， 得到 )()(121lim)(*mnxnxNmrNNnNxx(4.1.4) 第四章第四章 功功 率率 谱谱 估估 计计 将（4.1.4）式代入(4.1.1)式, 得到 mmnjnjNNnNnjmNmjxxmnxnxNmnxnxNP)(*e )(e )(121lime )()(121lim)e (令l=n+m, 则 2nj*je )(121lime )()(121lim)e (nxNlxenxNPNNnNljlNNnnjNxx(4.1.5) 第四章第四章 功功 率率 谱谱 估估 计计 上式中x(n)是观测数据，Pxx(

4、ej)是随机变量，必须对Pxx(ej)取统计平均值， 得到 2j)(121lim)e (NNnnjNxxenxNEP(4.1.6) 上式被认为是功率谱的另一定义。 （4.1.1）式表明功率谱是无限多个自相关函数的函数，但观测数据只有有限个,只能得到有限个自相关函数。按照（4.1.6）式求功率谱,也需要无限个观测数据。因此根据有限个样本数据，分析计算随机序列的真正功率谱，是求功率谱的中心问题，毫无疑问，这是一个功率谱的估计问题。在第一章已介绍了统计估计的一般估计准则，主要有偏移、估计量方差和估计量的均方误差（有效性），这里不再重复，下面直接用它们分析估计质量。 第四章第四章 功功 率率 谱谱 估

5、估 计计 现代谱估计以信号模型为基础，图4.1.1表示的是x(n)的信号模型，输入白噪声w(n)均值为0，方差为2w，x(n)的功率谱由下式计算： 2j2j| )e (|)e (HPwxx（4.1.7） 如果由观测数据能够估计出信号模型的参数，信号的功率谱可以按照（4.1.7）式计算出来，这样,估计功率谱的问题变成了由观测数据估计信号模型参数的问题。模型有很多种类，例如AR模型、 MA模型等等，针对不同的情况，合适地选择模型，功率谱估计质量比较经典谱估计的估计质量有很大的提高。遗憾的是，尚无任何理论能指导我们选择一个合适的模型，我们只能根据功率谱的一些先验知识，或者说一些重要的谱特性，选择模型

6、。第四章第四章 功功 率率 谱谱 估估 计计 H(z)w(n)x(n)图4.1.1 平稳随机序列的信号模型第四章第四章 功功 率率 谱谱 估估 计计 4.2 经经 典典 谱谱 估估 计计 4.2.1 BT法法 BT法是先估计自相关函数, 然后按照(4.1.1)式进行傅里叶变换得到功率谱。设对随机信号x(n)，只观测到一段样本数据， n=0, 1, 2, , N-1。 关于如何根据这一段样本数据估计自相关函数, 第一章已经作了详细介绍，结果是共有两种估计方法， 即有偏自相关函数估计和无偏自相关函数估计。有偏自相关函数估计的误差相对较小，这种估计是一种渐近一致估计， 将该估计公式重写如下： 1|0

7、*)()(1)(mNnxxmnxnxNmr（4.2.1） 第四章第四章 功功 率率 谱谱 估估 计计 对上式进行傅里叶变换，得到BT法的功率估计值为 mjxxjmrP-BTe )()e (4.2.2) 为了减少谱估计的方差，经常用窗函数w(m)对自相关函数进行加权， 此时谱估计公式为 jxxMmjmrP-)1(BTe )()e (（4.2.3） 式中 0)()(mwmw-(M-1)m(M-1) 其它 , MN （4.2.4） 第四章第四章 功功 率率 谱谱 估估 计计 有时称(4.2.3)式为加权协方差谱估计。它要求加窗后的功率谱仍是非负的，这样窗函数w(m)的选择必须满足一个原则，即它的傅里

8、叶变换必须是非负的， 例如巴特利特窗就满足这一条件。 为了采用FFT计算（4.2.3）式，设FFT的变换域为(0L-1)，必须将求和域(-M+1, M-1)移到（0L-1），功率谱的计算公式如下： )()(e )()e (BT10j -jBTmSFFTkPmSPxxLmmxx k=0, 1, 2, , L-1 第四章第四章 功功 率率 谱谱 估估 计计 (4.2.7) )()(0)()()(LmwLmrmwmrmSxxxxxx0mM-1 MmL-M-1 L-MmL-1 按照（4.2.1）式估计自相关函数，我们已经证明这是渐近一致估计，但经过傅里叶变换得到功率谱的估计，功率谱估计却不一定仍是渐近

9、一致估计，可以证明它是非一致估计，是一种不好的估计方法。下面我们将证明：BT法中用有偏自相关函数进行估计时，它和用周期图法估计功率谱是等价的，因此BT法估计质量和周期图法的估计质量是一样的。 第四章第四章 功功 率率 谱谱 估估 计计 4.2.2 4.2.2 周期图法周期图法 将功率谱的另一定义（4.1.6）式重写如下: 2jje )(121lim)e (NNnnNxxnxNEP如果忽略上式中求统计平均的运算，观测数据为:x(n) 0nN-1， 便得到周期图法的定义: 210j -je )(1)e (NnnxxnxNP(4.2.8) 第四章第四章 功功 率率 谱谱 估估 计计 图 4.2.1

10、用周期图法计算功率谱框图 观 测 数 据 x(n)FFT取 模的 平 方1/NPxx(ej)第四章第四章 功功 率率 谱谱 估估 计计 1. 1. 周期图与周期图与BTBT法的等价关系法的等价关系 （不讲）（不讲） 周期图法的功率谱估计公式用（4.2.8）式表示，下面由该公式出发推导它们的等价关系。 1010)(*10*10210e )()(1e )(e )(1e )(1)e (NkNnnkjNnnjNknjNnnjjxxnxkxNnxkxNnxNP令 m=k-n, 即k=m+n，则 )(10*je)()(1)e (nkjmNnxxnmxnxNP第四章第四章 功功 率率 谱谱 估估 计计 上式

11、中的方括号部分正是有偏自相关函数的计算公式, 因此得到 mjxxNNmxxmrP)e()e (1)1(j)()(1)(*|10nmxnxNmrmNnxx)e ()e (BTjjxxPP因此证明了利用有偏自相关函数的BT法和周期图法的等价关系。 第四章第四章 功功 率率 谱谱 估估 计计 2. 周期图法谱估计质量分析周期图法谱估计质量分析 1) 周期图的偏移 已知自相关函数的估计值 ，m=-(N-1), -N, -N+1, , 0, 1, 2, ,N-1，按照（4.2.2）式求功率谱的统计平均值， 得到 )(mrxxnxxNNmxxmrEPEj -1)1(je)()e (有偏自相关函数统计平均值

12、已由第一章(1.3.30)式确定，将该式代入上式,得到 nxxNNmxxmrNmNPEj -1)1(j)e(|)e (第四章第四章 功功 率率 谱谱 估估 计计 mnxxBxxmrmwPEj -j)e()()e (（4.2.9） 式中 其它0|)(NmNmNmwB（4.2.10） (4.2.9)式中， 两序列乘积的傅里叶变换， 在频域服从卷积关系， 得到 1)e (jxxPEWPEjBjxxd)e (21)e ()(（4.2.11） 式中 )()e (mrFTPxxjxx第四章第四章 功功 率率 谱谱 估估 计计 将得到的L个周期图进行平均，作为信号x(n)的功率谱估计， 公式如下： Liix

13、xILP1j)(1)e (4.2.19) 为了分析偏移，对上式求统计平均，得到 2j)-j(j-j)2/sin()2/sin(1)()e ()e ()e (21)e (MMmwFTWdPWPEBBxxBxx(4.2.20) (4.2.21) 周期图的统计平均值已经求出， 如（4.2.11）、 (4.2.12)式所示， 重写如下： 第四章第四章 功功 率率 谱谱 估估 计计 2)2/sin()2/sin(1)()e (NNmwFTWBjB(4.2.12) WB(ej)称为三角谱窗函数。（4.2.11）式表明，周期图的统计平均值等于它的真值卷积三角谱窗函数，因此周期图是有偏估计，但当N时，wB(m

14、)1，三角谱窗函数趋近于函数，周期图的统计平均值趋于它的真值，因此周期图属于渐近无偏估计。 第四章第四章 功功 率率 谱谱 估估 计计 2) 周期图的方差 由于周期图的方差的精确表示式很繁冗，为分析简单起见， 通常假设x(n)是实的零均值的正态白噪声信号，方差是x2，即功率谱是常数x2 ，其周期图用IN()表示，N表示观测数据的长度。 按照周期图的定义，周期图表示为 )()()(varee )()(1| )e (|1)(22j -1010j2jNNNnNnNkkNIEIEInxkxNXNI下面先求周期图的均值，再求其均方值： )(j -e )()()()(1)(n-kNkNnNkRnRkxnx

15、ENIE第四章第四章 功功 率率 谱谱 估估 计计 式中 knknknknknrkxnxExxx01)()()()()(2222)(1)(xnxNNnRNIE（4.2.13） 上式说明周期图是无偏估计，但前面已推导出周期图是有偏估计（一般情况），这里由于对信号作了实白噪声的假设，才有无偏估计的结果。在求均方值时，先求两个频率1和2处的均方值， 最后令=1=2。 第四章第四章 功功 率率 谱谱 估估 计计 )(j -)(j -22j2j2212121ee)()()()()()()()(1)e ()e (1)()(qpknNNNqNpknNNNNqxpxkxnxEqRpRkRnRNXXNEIIE利

16、用正态白噪声、多元正态随机变量的多阶矩公式，有 )()()()()()()()()()()()()()()()(pxkxEqxnxEqxkxEqxpxEqxpxEkxnxEqxpxkxnxE其它0,;,)()()(4kpnqkqnpqpnkqxkxnxEx第四章第四章 功功 率率 谱谱 估估 计计 将上式代入周期图的均方值公式中， 得到 22121221214101010)(10)(224212/ )sin(2/)sin(2/ )sin(2/)sin(1ee)()(2121NNNNNNIIExNnNnNkn-kjNkn-kjxNN(4.2.14) 将=1=2代入上式，得到 2442)sin()

17、sin(1)(var)sin()sin(2)(NNINNIExNxN(4.2.15) 第四章第四章 功功 率率 谱谱 估估 计计 显然，当N趋于无限大时，周期图的方差并不趋于0，而趋于功率谱真值的平方，即 4)(varxNNI (4.2.16) 这里无论怎样选择N，周期图的方差总是和4x同一个数量级。我们知道,信号的功率谱真值是2x,说明周期图的方差很大，周期图的均方误差也是非常大。用这种方法估计的功率谱在2x附近起伏很大，故周期图是非一致估计，是一种很差的功率谱估计方法。 第四章第四章 功功 率率 谱谱 估估 计计 为了进一步说明数据长度N对功率谱估计的影响，下面求两个频率处的协方差函数。

18、)()()()()(),(cov212121NNNNNNIIEIIIII将（4.2.13）式和(4.2.14)式代入上式, 得到 22121221214212/ )sin(2/ )(sin2/ )sin(2/ )(sin)()(covNNNNIIxNN令: 1=2k/N, 2=2l/N, 式中k, l均是整数，得到 22421/ )(sin)(sin/ )(sin)(sin)()(covNlkNlkNlkNlkIIxNN(4.2.17) 第四章第四章 功功 率率 谱谱 估估 计计 图 4.2.2 白噪声的周期图 第四章第四章 功功 率率 谱谱 估估 计计 4.2.3 经典谱估计方法改进经典谱估

19、计方法改进 1. 平均周期图法平均周期图法 平均周期图法是基于这样的思想：对一个随机变量进行观测，得到L组独立记录数据，用每一组数据求其均值，然后将L个均值加起来求平均。这样得到的均值，其方差将是用一组数据得到的均值的方差的1/L。 假设随机信号x(n)的观测数据区间为：0nM-1，共进行了L次独立观测，得到L组记录数据，每一组记录数据用xi(n), i=1, 2, 3, ,L表示，第i组的周期图用下式表示： 210)(1)(MnnjiienxMI(4.2.18) 第四章第四章 功功 率率 谱谱 估估 计计 上式表明,平均周期图仍然是有偏估计， 偏移和每一段的数据个数M有关；由于MN,平均周期

20、图的偏移比周期图的偏移大， 表现在三角谱窗主瓣的宽度比周期图主瓣的宽度宽。由于三角谱窗主瓣的宽度变宽，分辨率更加降低，因此也可以说， 偏移的大小反映分辨率的低与高。 按照（4.2.19）式求方差，由于是L次独立观测，L个周期图相互独立，因此平均周期图的方差为 )(var1)e (varijxxILP（4.2.22） 即平均周期图的估计方差是周期图的方差的1/L。显然，是以分辨率的降低换取了估计方差的减少，当然, 估计的均方误差也减少。 第四章第四章 功功 率率 谱谱 估估 计计 图 4.2.3 平均周期图法 Pxx(ej)/032101230321012303210123Pxx(ej)/Pxx

21、(ej)/(a)(b)(c)N256， L2N256， L4N256， L8第四章第四章 功功 率率 谱谱 估估 计计 2. 2. 窗函数法窗函数法（不讲）（不讲） 这种方法是用一适当的功率谱窗函数W(ej)与周期图进行卷积，来达到使周期图平滑的目的的。 d)e ()(21)e ()(jNlwxxWIP（4.2.23） 式中 mxxNNmNmrIj -1)1(e )()( 是有偏自相关函数 )(mrxxde )e (21)(jjnWnw-(M-1)nM-1 （4.2.24） 第四章第四章 功功 率率 谱谱 估估 计计 那么 mxxMMmjxxemwmrPj1)1()()()e (4.2.25)

22、 将(4.2.25)式和(4.2.3)式进行对比，它们是一样的，说明周期图的窗函数法就是前面提到的BT法的加权协方差谱估计。在（4.2.23）式中， 周期图和谱窗函数卷积得到功率谱， 等效于在频域对周期图进行修正，使周期图通过一个线性非频变系统， 滤除掉周期图中的快变成分，谱窗函数需具有低通特性。 第四章第四章 功功 率率 谱谱 估估 计计 对（4.2.25）式求统计平均, 得到 mxxMMmxxmwmrEPEj -1)1(je )()( )e (将(1.3.30)式和(1.3.32)式代入上式，得到 mBxxMMmxxmwmwmrePEj -1)1(je )()()()(4.2.26) 式中

23、 NmNmmwB|1)(第四章第四章 功功 率率 谱谱 估估 计计 上式表明, 周期图的窗函数法仍然是有偏估计, 其偏移和wB(m)、w(m)两个窗函数有关，如果w(m)窗的宽度比较窄， M比N小得多，这样|m|q时，(4.3.4)式却是一个线性方程， 用矩阵方程表示如下： )()(0*mlhlhmqlB第四章第四章 功功 率率 谱谱 估估 计计 )()2() 1()()2() 1 ()()2() 1()2()() 1() 1() 1()(pqrqrqrpqhhhqrpqrpqrpqrqrqrpqrqrqrxxxxxxAAAxxxxxxxxxxxxxxxxxx上式共有p个方程。可以用该方程首先

24、计算出AR部分的p个系数hA(i), i=1, 2, 3, p; 然后代入(4.3.4)式，设法求出MA部分的系数。 第四章第四章 功功 率率 谱谱 估估 计计 2.AR2.AR模型的系数和信号自相关函数之间的关系模型的系数和信号自相关函数之间的关系 AR模型的系统函数为：H(z)=1/A(z)，相当于ARMA模型中B(z)=1的情况，这样在公式中 mlmlmlmlhnnhBB01)()()()(因为h(n)是因果性的，因此m0时 h*(-m)=0, 将上面公式代入到（4.3.4）式中, 得到 PlwxxAPlxxAxxlrlhmrlhmr121)()() 1()()(m1 m=0 (4.3.

25、6) 第四章第四章 功功 率率 谱谱 估估 计计 也可以将上式中m1的情况写成矩阵形式： )()2() 1 ()()2() 1 ()0()2() 1()2()0() 1 () 1() 1()0(prrrphhhrprprprrrprrrxxxxxxAAAxxxxxxxxxxxxxxxxxx(4.3.7) （4.3.6）式的矩阵形式如下式： 001)0() 1()()2()0() 1 () 1() 1()0(21wpxxxxxxxxxxxxxxxxxxharprprprrrprrr(4.3.8) 第四章第四章 功功 率率 谱谱 估估 计计 或者用模型参数表示： 001)0() 1()() 1()

26、0() 1 ()() 1()0(21wpxxxxxxxxxxxxxxxxxxharprprprrrprrr(4.3.9) 令 )0() 1()() 1()0() 1 ()() 1()0()1(xxxxxxxxxxxxxxxxxxpxxrprprprrrprrrR(4.3.10) 第四章第四章 功功 率率 谱谱 估估 计计 称为自相关矩阵，它满足（H表示共轭转置）， 是一个埃尔米特（Hermitian）矩阵，且沿任一对角线的元素相等，是一个托布列斯（Toeplitz）矩阵。 也是正定矩阵。 上面推导出的（4.3.6）式或(4.3.9)式确定了AR模型参数（包括模型输入噪声方差）和信号自相关函数之

27、间的关系。 我们注意到这是一个线性方程，如果能够由信号的观测数据求出其自相关函数，可以按照（4.3.7）式,通过解一组线性方程得到模型参数,相对ARMA模型，这是AR模型的优点。 )1( pxxRHxxxxRR)1( pxxR第四章第四章 功功 率率 谱谱 估估 计计 3.MA3.MA模型的系数和信号自相关函数之间的关系模型的系数和信号自相关函数之间的关系 MA模型的系统函数H(z)=B(z)，相当于ARMA模型中A(z)=1, hA(n)=(n)的情况，此时h(n)=hB(n)，由（4.3.4）式得到MA模型系数和信号自相关函数的关系为 0)()()(0*2mqlBBwxxmlhlhmrm=

28、0, 1, , q mq+1 (4.3.11) 上式表明，MA模型的参数和信号自相关函数之间也是非线性关系。 第四章第四章 功功 率率 谱谱 估估 计计 上面我们分别推导了三种信号模型的参数和信号自相关函数之间的关系。这些关系式为我们提供了一种估计功率谱的方法，即首先根据信号观测数据估计信号自相关函数，再按照所选择信号模型，解上面相应的方程，求出模型参数，最后按照下式求出信号的功率谱: 21022j2je1e| )e (|)e (PiiiqiiiwwxxabHP(4.3.12) 第四章第四章 功功 率率 谱谱 估估 计计 4.4 AR模型谱估计的性质模型谱估计的性质4.4.1 AR4.4.1

29、AR模型的线性预测模型的线性预测 在第二章中，我们已推导出维纳线性一步预测器系数和信号自相关函数之间的关系式(称为Yule-Walker方程)，重写如下： 00)(1)0() 1()() 1()0() 1 ()() 1()0(min21neEaarprprprrrprrrpppxxxxxxxxxxxxxxxxxx(4.4.1) 第四章第四章 功功 率率 谱谱 估估 计计 式中, e(n)表示线性一步预测误差，其公式为 pipiinxanxnxnxne1)()()( )()(4.4.2) Ee2(n)min表示e(n)的均方差最小值；api(i=1, 2, 3, ，p)表示预测器的系数，它和线性

30、预测器单位脉冲响应h(n)之间差一符号，即pnanhpn, 3 , 2 , 1)( 对（4.4.2）式进行Z变换, 得到 piipizazXzXzE1)()()(4.4.3) 第四章第四章 功功 率率 谱谱 估估 计计 令He(z)=E(z)/X(z)，由上式，得到 piipiezazH11)(4.4.4) 称He(z)为线性一步预测误差滤波器，其作用是将信号x(n)转换成预测误差e(n)，如图4.4.1所示。一般认为e(n)具有白噪声的性质，因此He(z)也称为白化滤波器。 第四章第四章 功功 率率 谱谱 估估 计计 图 4.4.1 预测误差滤波器 He(z)x(n)e(n)x(n)e(n)

31、z1z1z1ap1ap2ap, p1ap, p第四章第四章 功功 率率 谱谱 估估 计计 我们知道AR模型的系统函数为 PiiizazH111)(4.4.5) 对比（4.4.4）、(4.4.5)两式，当api=ai(i=1, 2, 3, ,p)时， He(z)和H(z)互为逆滤波器，He(z)=1/H(z), 基于以上分析，也可以将AR模型定义为 )()( )(nwnxnx(4.4.6) 式中, 是基于信号前p个样本的最佳一步线性预测，公式为 )( nxpiipizanx1)( 第四章第四章 功功 率率 谱谱 估估 计计 w(n)是模型输入白噪声。将AR模型参数和信号自相关函数之间的关系式（4

32、.3.9）和(4.4.1)式进行对比, 得到: api=ai, w(n)=e(n),Ee2(n)min=2w，信号自相关函数和它的AR模型参数之间的关系服从Yule-Walker方程。注意:只有当AR模型的阶数与线性预测器的阶数相同时，以上结论才是正确的。由于AR模型具有这种特性，因而AR模型法也称为线性预测AR模型法。 第四章第四章 功功 率率 谱谱 估估 计计 AR模型与线性预测之间的关系，可以被用来解卷积，假设信号s(n)通过一个AR系统，系统单位脉冲响应为h(n)，响应是x(n)，即 x(n)=s(n) * h(n) 如果s(n)具有白噪声性质，可以利用AR模型与预测滤波器之间的关系对

33、上式进行解卷积，得到s(n)信号。方法是：先由x(n)的观测数据估计AR模型的参数，得到AR模型的系统函数H(z)，再让x(n)通过H(z)的逆滤波器H-1(z)，便得到信号s(n)。 第四章第四章 功功 率率 谱谱 估估 计计 4.4.2 预测误差滤波器的最小相位特性预测误差滤波器的最小相位特性 我们知道, AR模型必须因果稳定，即极点均在单位圆内， 才能保证信号x(n)是平稳随机信号，于是AR模型H(z)和预测误差滤波器He(z)互为逆滤波，那么He(z)应为最小相位系统。但是由解Yule-Walker方程得到AR模型的参数，其极点不一定在单位圆内。下面将证明当最佳P阶线性预测系数与AR模

34、型参数相同时，由此得到的极点保证在单位圆内，AR滤波器稳定，预测误差滤波器He(z)或者A(z)是最小相位系统。这里自相关矩阵 是正定的。 )1( pxxR第四章第四章 功功 率率 谱谱 估估 计计 解解Yule-Walker方程得到的是最佳线性预测滤波器的系数， 此时预测误差滤波器输出功率Pe达到最小，用Pe min表示，即 )d(e| )(e|j-2jminxxePAP(4.4.7) 式中, ,ai是最佳预测系数。下面先用反证法证明A(z)的全部零点不在单位圆外部（即全部零点在单位圆上或者单位圆内部）。设A(z)的第i个零点zi在单位圆外部， 即|zi|1，用1/z*i代替zi，这时A(z

35、)的幅度函数不变，按照（4.4.7）式计算出的预测误差滤波器输出功率Pe应不变，仍是最小的。 但是下面将用公式证明预测误差滤波器输出功率Pe不是最小的， 这一矛盾的结论只能说明A(z)的零点不可能在单位圆外部。 piiizazA11)(第四章第四章 功功 率率 谱谱 估估 计计 )( )1 ()1 ()1 ()1 (1)(1111111zAzzzzzzzzzazApijjjjipjjpiii(4.4.8) 式中 pijjjzzzA11)1 ()( 将(4.4.8)式代入（4.4.7）式, 得到 d)e (| )e ( |e1j2j2j -minxxiePAzP(4.4.9) 第四章第四章 功功

36、 率率 谱谱 估估 计计 式中 2*22*j22j22j11|1|e1|e1 |jiiiiiiiezzzezzzz因为 1iz因此 2j -*2j -e11|e1 |iizz(4.4.10) 第四章第四章 功功 率率 谱谱 估估 计计 4.4.3 AR4.4.3 AR模型隐含自相关函数延拓特性模型隐含自相关函数延拓特性 AR模型的自相关函数和模型系数之间的关系服从Yule-Walker方程，重写如下： plwxxAplxxAxxlrlhlmrlhmr121)()()()()(m1 m=0 上式中，对于m1的情况，公式本身就是一个递推方程，如果已由观测数据计算出p+1个自相关函数，用 ,m=0,

37、1,2, , ,p表示，对于mp的情况, 可以用该公式外推得到，公式如下： )(mrxx第四章第四章 功功 率率 谱谱 估估 计计 PlxxAxxxxlmrlhmrmr1)( )()()(0mp 0p (4.4.11) 上式中，系数hA(l)需用（4.3.6）式求出。因此AR模型隐含着自相关函数外推的特性。我们知道,经典谱估计BT法中，自相关函数只能限于由观测数据计算出的有限个自相关函数，其它的认为是0，造成了谱估计分辨率低、模糊。也正是AR模型具有自相关函数外推特性， 使它具有高分辨率的优点。 第四章第四章 功功 率率 谱谱 估估 计计 4.5 AR谱估计的方法谱估计的方法4.5.1 4.5

38、.1 自相关法自相关法列文森（列文森（LevensonLevenson）递推法）递推法 自相关法的出发点是选择AR模型的参数使预测误差功率最小，预测误差功率为 npipininxanxNneN122)()(|1| )(|1(4.5.1a) 假设信号x(n)的数据区在0nN-1范围，有p个预测系数，N个数据经过冲激响应为api(i=0, 1, 2, , p)的滤波器， 输出预测误差e(n)的长度为N+P， 因此应用下式计算: 第四章第四章 功功 率率 谱谱 估估 计计 2110210| )()(|1| )(|1pipiPNnPNninxanxNneN(4.5.1b) 显然,e(n)的长度长于数据

39、的长度，上式中数据x（n）的两端需补充零点，这相当于无穷长的信号经过加窗处理，得到长度为N的数据。 用（4.5.1b）式对系数api的实部和虚部求微分的方法使预测误差功率最小，得到 )()2() 1 ()0()2() 1()2()0() 1 () 1() 1()0(21prrraaarprprprrrprrrxxxxxxppppxxxxxxxxxxxxxxxxxx(4.5.2) 第四章第四章 功功 率率 谱谱 估估 计计 式中自相关函数采用有偏自相关估计， 即 )()()(1)(*10mrmnxnxNmrxxmNnxxm=0, 1, 2, , p m=-p+1, -p+2, , -1 (4.5

40、.3) 对比（4.3.7）式，（4.5.2）式就是已推导出的Yule-Walker方程，因此自相关法也是基于解Yule-Walker方程的一种方法。 首先由信号的观测数据估计出其自相关函数，再解该方程，得到模型参数，便可求出信号的功率谱。 因此该方法也称为Yule-Walker法。但是直接解该方程，需要计算逆在矩阵，不方便。在第三章自适应滤波器中，曾介绍了基于Yule-Walker方程中自相关矩阵的性质，导出Levenson-Durbin递推法，这是一种高效的解方程方法。下面把已推出的LevensonCD*2Durbin递推法简化重写如下： 第四章第四章 功功 率率 谱谱 估估 计计 2111

41、, 1/)()(kklxxlkxxkklkrakra（4.5.4） *, 1, 1ikkkkikkiaaaai=1, 2, 3, , k-1 （4.5.5） 2122)|1 (kkkka（4.5.6） 由k=1开始递推，递推到k=p，依次得到a11,21，a21,a22, 22，ap1,ap2,app,2p。AR模型的各个系数以及模型输入白噪声方差求出后， 信号功率谱用下式计算： 21j -22j2je11| )e (|)e (PiiiwwxxaHP（4.5.7） 第四章第四章 功功 率率 谱谱 估估 计计 （4.5.6）式表明： ，说明随着阶数增加， 预测误差功率将减少或者不变，为此要求|a

42、kk|1,akk称为反射系数。另外, 递推公式提供了一种确定模型阶数的实验方法， 如模型的阶数不知道，由低阶开始递推，当递推到M阶时，预测误差满足允许的值，停止递推，选AR模型的阶数为M。 这种递推法效率高， 且当阶数变化时， 无需从头计算。 如果知道信号的 N个观测数据（x(n),0nN-1 ），利用列文森递推法计算功率谱的计算流程图如图4.5.1所示。图中 采用有偏自相关估计（4.5.3）式计算，2是要求的方差。 2232221p)(mrxx第四章第四章 功功 率率 谱谱 估估 计计 图 4.5.1 利用列文森递推法计算功率谱的流程图 输入： x(n),n0,1,2,3,N1,2计算：)0

43、(1),0(/ ) 1 (21 , 1211 , 1xxxxxxrarra)(mrxxp=22122*, 1, 12111, 111, 2 , 1 , 0)()(ppppkppppkppkppixxipxxppapkaaaaiprapra22ppiipiwxxaP1j2je11)(e结束p=p+1NY第四章第四章 功功 率率 谱谱 估估 计计 4.5.2 4.5.2 协方差法与修正协方差法协方差法与修正协方差法 1. 1. 协方差法协方差法 这种方法和自相关法一样，仍然利用使预测误差功率最小的方法求模型参数 ，但由观测数据求预测误差功率的公式如下式： 11212| )()(|1| )(|1Np

44、nNpnpkpkknxanxpNnepN（4.5.8） 将该式对比自相关法中求预测误差功率的公式（4.5.1b），不同的是求和限不同。 该公式中使用的观测数据均已得到，不需要在数据两端补充零点， 因此比较自相关法去掉了加窗处理的不合理假设。为求得模型参数仍然应用复梯度法使（4.5.8）式达到最小， 公式如下： 第四章第四章 功功 率率 谱谱 估估 计计 )0 ,()0 , 2()0 , 1 (),()2 ,() 1 ,(), 2()2 , 2() 1 , 2(), 1 ()2 , 1 () 1 , 1 (21pcccaaappcpcpcpcccpcccxxxxxxppppxxxxxxxxxxx

45、xxxxxxx（4.5.9） 式中 1*)()(1),(NpnxxknxjnxpNkjc（4.5.10） 白噪声的方差为 1min2), 0()0 , 0(Npkxxpkxxwkcac（4.5.11） 第四章第四章 功功 率率 谱谱 估估 计计 由观测数据x(n)(n=0, 1, 2, , N-1)，利用上面三个公式可以求出模型的参数：api(i=1, 2, 3, , p); 2w。按照定义， （4.5.9）式中的cxx(j, k)可以称作协方差函数， 它有两个变量， 因此也适合于非平稳随机信号。（4.5.9）式中的协方差矩阵是埃尔米特（Hermitian）矩阵，cxx(k,j)=c*xx(j

46、,k)，是半正定的。 这种方法近似于自相关法。一些实验结果说明它的分辨率优于自相关法5，另外对于纯正弦信号数据，可以有效地估计正弦信号的频率。 进一步的内容可参考文献14、15 第四章第四章 功功 率率 谱谱 估估 计计 2. 2. 修正协方差法修正协方差法 修正协方差法使用前向和后向预测误差平均值最小的方法， 估计AR模型的参数，进而估计信号的功率谱。该方法最初由Nuttall在1976年提出，称为前向后向法；同年,Ulrych和Ctayton也独立提出，称之为最小二乘法，因此该法也称为前后向线性预测误差功率最小的最小二乘法。 前面我们已推导出信号的前向和后向预测， 分别重写如下： pkpk

47、pkpkknxanxknxanx1*1)()( )()( (4.5.12) (4.5.13) 第四章第四章 功功 率率 谱谱 估估 计计 式中apk是AR模型的参数。最小预测误差平均功率是模型输入白噪声的方差，即p=2w，前后向预测误差平均功率为 )(5 . 0pbpep(4.5.14) 式中前向和后向预测误差功率pe、pb分别用下式表示： 2101*211)()(1)()(1pNnpkpkpbNpnpkpkpeknxanxpNknxanxpN(4.5.15) (4.5.16) 第四章第四章 功功 率率 谱谱 估估 计计 和协方差法一样，上式中仅对那些用到的观测数据预测误差求和。为了使预测误差

48、平均功率最小，求p对apk(k=1, 2, 3, , p)的实部和虚部的微分,或者用复梯度法求，得到 pllnxknxanxlnxknxanxpNapkpkpNnNpnpkpkplp, 3 , 2 , 10)()()()()()(11*101*1(4.5.17) 第四章第四章 功功 率率 谱谱 估估 计计 经过简化，得到 1*10*11*10*)()()()()()()()(NpnpNnpkNpnpNnpklnxnxlnxnxlnxknxlnxknxa令 1*10*)()()()()(21),(NpnpNnxxknxjnxknxjnxpNkjc(4.5.18) 第四章第四章 功功 率率 谱谱

49、估估 计计 将上式写成矩阵形式： )0 ,()0 , 2()0 , 1 (),()2 ,() 1 ,(), 2()2 , 2() 1 , 2(), 1 ()2 , 1 () 1 , 1 (21pcccaaappcpcpcpcccpcccxxxxxxppppxxxxxxxxxxxxxxxxxx(4.5.19) 白噪声的方差估计值为 pkxxpkxxpNnpkpkNpnpkpkpwkcacnxknxanxnxknxanxpN1101*1*1min,2), 0()0 , 0()()()()()()()(21(4.5.20) 第四章第四章 功功 率率 谱谱 估估 计计 例例4.5.1 已知信号的四个观

50、察数据为x(n)=x(0), x(1), x(2), x(3)=2, 4, 1, 3， 分别用自相关法和协方差法估计AR（1）模型参数。 解解（1） 自相关法： 5 . 009)3()41 (4)24(20) 1()(21)()(21) 1()()(| )(|4111111111114011401111402aaaaanxneaneneanxanxnenennn第四章第四章 功功 率率 谱谱 估估 计计 (2) 协方差法: 714. 00)3()41 (4)24(20) 1()(32)()(32) 1()()(| )(|31| )(|111111111311131111131212aaaanx

51、neaneneanxanxnenenepNnnnNpn第四章第四章 功功 率率 谱谱 估估 计计 4.5.3 伯格（Burg）递推法 设信号x(n)观测数据区间为：0nN-1，前向、后向预测误差功率分别用p,e和p,b表示，预测误差平均功率用p表示， 公式分别为 )(21| )(|1| )(|1,12,12,bpeppNpnbpbpNpnfpepnepNnepN(4.5.21) (4.5.22) (4.5.23) 第四章第四章 功功 率率 谱谱 估估 计计 下面解释（4.5.21）、(4.5.22)式中的求和限。我们已经知道， 前向、后向预测误差公式分别为 pkpkbppkpkfpkpnxap

52、nxneknxanxne1*1)()()()()()(4.5.24) (4.5.25) 在上面两个公式中，信号项的自变量最大的是n，最小的是n-p。为了保证计算范围不超出给定的数据范围，在（4.5.21）和(4.5.22)式中， 选择求和范围为：pnN-1。 第四章第四章 功功 率率 谱谱 估估 计计 为求预测误差平均功率p最小时的反射系数kp，令 0ppk(4.5.26) 将第三章中已推导出的前向、后向预测误差递推公式重写如下： )() 1()() 1()()(1*111nekneneneknenefppbpbpbppfpfp(4.5.27) (4.5.28) 将上面两个公式带入（4.5.2

53、1）、 (4.5.22)、（4.5.23）式中，得到 121*1211| )()(| ) 1()(|)(21NpnfppfpbppfppnekneneknepN第四章第四章 功功 率率 谱谱 估估 计计 利用复梯度法，求p对kp的实部和虚部的微分，并令其结果等于0，得到 0)()() 1()1() 1()()(11*1*11*111nenekneneneknepNkfpfppbpNpnbpbppfppp解出kp，得到 )| ) 1(| )(|) 1()(2121211*11NpnbpfpNpnbpfppnenenenek(4.5.29) 第四章第四章 功功 率率 谱谱 估估 计计 上式就是利用

54、伯格递推法求第p个反射系数的公式。将伯格递推法求AR模型参数的递推公式总结如下： )| ) 1(| )(|) 1()(21, 3 , 2 , 1 , 0)()(1, 3 , 2 , 1 , 0)()()0(| )(|1)0(211211*11000102NpnbpfpNpnbpfppbfxxNnxxnenenenekNnnxneNnnxnernxNr(4.5.30) (4.5.31) (4.5.32) (4.5.33) (4.5.34) 第四章第四章 功功 率率 谱谱 估估 计计 2, 1,)() 1()(1, 2, 1)()()(1, 3 , 2 , 1)|1 (1*1111,*, 1, 1

55、12NppnnekneneNppnneknenekapiakaakfppbpbpbppfpfppppipppippippp(4.5.36) (4.5.35) (4.5.37) (4.5.38) (4.5.39) 对于p阶AR模型的输入白噪声方差2w=p。利用伯格递推法求AR模型参数的流程图如图4.5.2所示。 第四章第四章 功功 率率 谱谱 估估 计计 图 4.5.2 伯格递推法流程图 输入： x(n),n0,1,2,3,N1阶数IP1)(1)()()()(210000pnxNnxnenxneNnbfp = IP结束1 pp12121*111)()() 1()(2NpnbpfpbpNpnfpp

56、nenenenekpwpppk2121pppipppippikaakaa*, 1, 11, 2, 1),() 1()(1, 2, 1),1()()(1*111NppnnekneneNppnneknenefppbpbpbpppfp2,;, 3 , 2 , 1:wippia输出21j2je1)(epiipiwxxaPNY第四章第四章 功功 率率 谱谱 估估 计计 4.5.4 4.5.4 关于关于ARAR模型阶次的选择模型阶次的选择 在AR模型谱估计中，模型阶次的选择也是一个关键问题。 一般模型的最好选择是先验未知的，实际中需预先选定模型阶次。如果是纯P阶AR信号，选择模型阶次kP时，将产生对谱的平

57、滑作用，降低谱的分辨率，如图4.5.3所示。图中，AR信号P=4，选择模型的阶次k=2，产生的平滑作用使两个峰变成一个峰，分辨率明显降低。如果选择kP，且假定观测的数据没有误差（没有干扰），估计的参数应是： kppipiaaipki, 2, 10, 3 , 2 , 1,第四章第四章 功功 率率 谱谱 估估 计计 图 4.5.3 AR模型阶次太小时的平滑作用50.0030.0010.0010.0030.000.0000.400.50真实AR(4)的PSDAR(2)的PSD模型PSD/dB/第四章第四章 功功 率率 谱谱 估估 计计 因此,对于纯P阶AR信号，应选择阶次kP。

58、一般对于纯AR信号，只要记录数据的长度不十分短， 模型阶次的估计结果比较满意21。但是如果是白噪声中的AR信号（观测数据有误差或者信号中含有白噪声），此时选择ARMA模型合适，如选择了AR模型，其阶次应加大，较低的阶次会使谱估计产生偏移， 降低分辨率。当然, 这也和信噪比有关， 信噪比愈低，平滑作用愈严重，愈需要高的阶次， 因此信噪比低应选高的阶次。一般来说，阶次愈高，分辨率愈高；但阶次太高，会使估计误差加大，谱峰分裂，因此， 对于白噪中的AR信号， 其阶次的选择应折衷考虑5。 下面介绍三种模型阶次的估计方法。 第四章第四章 功功 率率 谱谱 估估 计计 为了确定模型的阶次，所提出的各种方法都

59、是基于估计出的预测误差功率，即选择模型的阶次，使预测误差功率极小化， 这是因为所讨论的AR模型参数估计方法，预测误差功率都是随模型的阶次的增加而减少，或者不变的。但是不能仅仅考虑预测误差功率的减少，也要考虑因模型阶次增高引起模型参数加多，引起估计误差的加大。根据这一原理，阿凯克提出两种估计方法，第一种方法称为最终预测误差（FPE）准则，用下式计算后再取最小，估计模型阶次： kkNkNkFPE)(（4.5.40） 第四章第四章 功功 率率 谱谱 估估 计计 式中, 表示k阶AR模型的白噪声方差（预测误差功率）的估计值;N是观测数据的长度; k是模型阶次。上式中, 随k增加而减少，但(N+k)/(

60、N-k)却随k增加而加大。第二种方法应用更广泛， 称为阿凯克信息论准则： kp kp knkk2ln)(AIC（4.5.41） 选择使上式最小的k值作为模型的阶次。但这种方法有过高估计模型阶次的趋势，例如真实的阶次是2，而AIC(k)最小值对应的是4，为改进, 提出用k lnN代替公式中的2k的方法22、 23。 AIC和FPE估计阶次的性能是相近的， 但对于短记录数据，建议使用AIC准则24，对于长记录（N）, 两种方法得到相同的模型阶次。 第四章第四章 功功 率率 谱谱 估估 计计 第三种方法由Parzen在1976年提出， 称为自回归传递函数准则（CAT）25，定义是： kkiiNk11