生医信号处理课件7现代功率谱

1、信自学院生医系信自学院生医系生物医学信号处理第七章 功率谱估计的现代方法7.1 概述经典谱估计：以傅立叶变换为基础，优点：计算效率高；缺点：现代谱估计Yule于1927年提出用线性回归方程来模拟时间序列参数模型法的基础1938年首次建立自回归模型参数与自相关函数关系的Yule-Walker方程Bartlett于1948年首次提出用自回归模型系数来计算功率谱1967年，Burg提出的最大熵谱估计参数模型谱估计 AR模型、MA模型、ARMA模型、PRONY指数模型等非参数模型谱估计 最小方差方法、MUSIC方法1965年，Cooley和Tukey的FFT算法问世，促进了现代谱估计的迅速发展7.2

2、谱估计的参数模型法7.2.1 三个步骤1、为被估计的随机过程确定一个合理的模型，一般用“有理传输函数模型”2、根据已知的观测数据估计模型的参数3、用模型参数计算功率谱描述线性时不变系统的差分方程为：两边进行Z变换得输出x(n)的功率谱为7.2.2 三种参数模型一、MA(q)模型（又称q阶滑动平均模型）二、AR(p)模型（又称p阶自回归模型）三、ARMA(p,q)模型应用较多的是AR模型，因为建立这种模型的计算比较容易，作为数学逼近，三种模型可以相互转换。7.3 AR模型7.3.1 AR模型的Yule-Walker方程N个样值x(0),x(1)x(N)自相关函数R(0),R(1).R(N)AR模

3、型参数和激励源方差a(1),a(2)a(p)功率谱密度YW方程 Yule-Walker方程的推导对 进行求逆z变换 把模型的差分方程代入x(n)的自相关函数 )()()(mnxnxEmRx)()()(12zAzAzSx)()()(1nuknxanxpkk)(nh0, 00,)(22mmmh如果 是因果的，则结论：结论：AR模型的Yule-Walker方程Yule-Walker方程的矩阵形式0.0.1)0(.)1()(.)1(.)0()1()(.)1()0(21paaRpRpRpRRRpRRR已知（或估计出）p+1个自相关函数值R(0),R(1),R(p)可求出 p+1个模型参数 221,.,p

4、aaaYule-Walker方程的求解采用高斯消元法，解方程组，运算的量级为p的三次方。用Levinson-Durbin算法，运算量的量级为p的二次方。7.3.2 Levinson-Durbin算法 以AR(0),AR(1)模型参数为初始条件，求出AR(2)的模型参数，再以这个为基础，求出AR(3)的模型参数，最后求出AR(p)的模型参数) 1 (001)0()2() 1()() 1() 1 ()0() 1 ()()2() 1 ()0(1paarprprprprrrrprrrr核心是求解如下方程：核心是求解如下方程：这个方程实际上是联合这个方程实际上是联合ARAR模型法和预测滤波法得出的。模型

5、法和预测滤波法得出的。我们发现，方程我们发现，方程(1)(1)有如下特点：有如下特点： 系数矩阵是一个系数矩阵是一个ToplitzToplitz矩阵，利用矩阵，利用ToplitzToplitz矩阵的性矩阵的性质可简化方程求解。质可简化方程求解。 实际问题中，一般只知道信号的某些观测值，而不知实际问题中，一般只知道信号的某些观测值，而不知道其道其ARAR模型阶数，该阶数也需要在方程求解过程中找到。模型阶数，该阶数也需要在方程求解过程中找到。 假设已得到假设已得到k k阶线性预测系数阶线性预测系数( (预测滤波参数预测滤波参数) )，我，我们来考虑求们来考虑求k k1 1阶滤波参数。阶滤波参数。k

6、 k阶滤波参数的矩阵方程为阶滤波参数的矩阵方程为)2(001)0()2() 1()() 1() 1 ()0() 1 ()()2() 1 ()0()()()(1kkkkaarkrkrkrkrrrrkrrrr由于系数矩阵的由于系数矩阵的ToplitzToplitz性质，上式又有如下形式性质，上式又有如下形式)2(001)0()2() 1()() 1() 1 ()0() 1 ()()2() 1 ()0()()(1)(kkkkaarkrkrkrkrrrrkrrrr现考虑模型阶数增加现考虑模型阶数增加1, 1, 即从即从k k变为变为k+1k+1的情况。的情况。对于对于k+1k+1模型模型, , 有有)

7、3(0001)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0()1()1(1)1()1(1)1(kkkkkkkaaarrkrkrrrkrkrkrkrrrkrkrrrR由于系数矩阵的由于系数矩阵的ToplitzToplitz性质，性质，k+1k+1阶系数矩阵阶系数矩阵 可有两种分块形式。可有两种分块形式。 )1( kR 利用这个性质，可设利用这个性质，可设)4(001001)()()(1)()()(1)1(kkkkkkkkkaaaaR式中式中)5()1() 1(1) 1 ()() 1(1)()()(1kjkjTkkkjkrakraar

8、krkr 比较比较(3)(3)和和(4)(4)，可知，当，可知，当(3)(3)与与(4)(4)等效；且有下列两个递推关系式：等效；且有下列两个递推关系式：)6(0)(k即当下式成立时即当下式成立时)7(/ )1() 1(/)(1)()(kkjkjkjkrakr)8(10011)(1)()()(1)1(1)1()1(1kkkkkkkkkkkaaaaaaa和)9()()1(kk由由(8)(8)末式还可得末式还可得: :)10()1(1kka(5)(5)(10)(10)构成构成LevinsonLevinson算法基础算法基础。 现用现用i i表示递推过程的阶数，令表示递推过程的阶数，令i=k+i=k

9、+1 1, , 并设信号模并设信号模型的最大阶数为型的最大阶数为N N，则有如下则有如下LevinsonLevinson算法算法: : 1) 1) 由由(3)(3)式，令式，令i=k+i=k+1 10 0，得，得 2)2) 置置i=k+i=k+1 11 1; ; 3) 3) 由由(8)(8)、(10)(10)式计算式计算 4)4) 由由(7)(7)、(10)(10)式计算式计算 5)5)由由(6)(6)、(7)(7)和和(10)(10)式计算式计算 6)6) 置置i i = =i i+1+1; ; 7)7) 判别：若判别：若 转转3)3)；否则，结束程序。；否则，结束程序。)0()0(r)11

10、(/ )()()1(11)1()(iijijiijiraira)12() 1,.,1()1()()1()(ijaaaaijiiiijij)13()(1 2)()1()(iiiiaNi v 讨论讨论LevinsonLevinson算法第算法第4 4步利用了一个重要递推关系步利用了一个重要递推关系(12)(12)， 通常称为通常称为LevinsonLevinson关系式关系式递推过程产生一个滤波参数序列递推过程产生一个滤波参数序列 通常称为偏相关系数通常称为偏相关系数递推过程产生的递推过程产生的 可用来监视可用来监视i i阶信号模型的阶信号模型的均方误差估值。均方误差估值。递推结果的最终解为递推结

11、果的最终解为 和和),.,1()(Niaii),.,1()(NjaNj)(N)(i递推过程及结果递推过程及结果 讨论讨论 优点：计算简单优点：计算简单 缺点：需根据有限观测数据估计自相关序列缺点：需根据有限观测数据估计自相关序列r(n)r(n) 短数据序列时，自相关估计值误差很大，引起预测短数据序列时，自相关估计值误差很大，引起预测 滤波参数误差，导致滤波参数误差，导致“谱峰飘移谱峰飘移”和和“谱线分裂谱线分裂”( (即出即出 现虚假谱线现虚假谱线) ) 长数据序列时，自相关估计值虽精确，但计算量大。长数据序列时，自相关估计值虽精确，但计算量大。优缺点优缺点确定AR模型的阶次piakppijk

12、ipa,.2, 1.2, 1.0,估计一个AR(p）过程，选取AR(k），要求k=p , k不能太大，如果精确的话，kp,确定AR模型的阶的方法不断增加阶数，观察预测误差功率，下降到一定的程度不断增加阶数，观察各阶模型预测误差序列的周期图，最接近平坦时发现后面的很多AR模型的参数接近于0FPE(最终预测误差）N为观测数据长度，为观测数据长度， 为方差，为方差，使使FPE最小最小 2kBergBerg算法算法 BergBerg算法原理算法原理根据前面的基本概念，可知根据前面的基本概念，可知m m阶前向预测误差为阶前向预测误差为)14(.,2 , 1, )()(0apminxanfmimim类似地

13、，类似地， m m阶后向预测误差为阶后向预测误差为)14(.,2 , 1, )()(0*bpminxangmimimm再利用再利用LevinsonLevinson关系式：关系式：)15() 1,.,1( ,)()1()1()(mjaaaammmmjmmmjmj有有)16()() 1()()16() 1()()(1*111bnfngngangnfnfmmmmmmmm其中其中pmnxngnf,.,2 , 1; )()()(00 定义定义m m阶前、后向预测误差的功率为阶前、后向预测误差的功率为)17()()(21P22pmnmmmngnf将将(16)(16)代入代入(17)(17)，并令，并令Pm

14、对对 的偏导数为零的偏导数为零, ,得最佳得最佳mm)18(,.,2 , 1,) 1()() 1()(12*121*1pmngnfngnfpmnmmpmnmmm 设已知有限数据序列设已知有限数据序列x(n),n=0,1,Nx(n),n=0,1,N，则可按下列步骤计算预，则可按下列步骤计算预测滤波器系数，并在此基础上计算功率谱。测滤波器系数，并在此基础上计算功率谱。1. 1. 置置m=0, m=0, 计算初值计算初值)()()()(00200nxngnfnxPNn2 2. m=m+1m=m+1,并按并按(19)(19)计算反射系数计算反射系数m3.3.计算滤波器系数：计算滤波器系数：) 1,.,1( ,)1()1()(mjaaamjmmmjmj4. 4. 计算预测误差功