20_经典谱估计与参数模型法

第四章功率谱估计§4.1引言研究信号与系统有两类基本方法:时域法和频域法.由于随机信号不存在傅氏变换,因此在频域主要研究它的功率谱.功率谱的两种定义:（1）定义之一——维纳-辛钦定理根据维纳-辛钦定理,信号的功率谱和它的自相关函数构成一对傅氏变换:正变换(4.1.1)逆变换(4.1.2))(nx)(jxxeP)(mRxxmjmxxjxxemReP)()(deePmRmjjxxxx)(21)(该式建立了平稳随机序列的自相关函数与功率谱密度之间的关系.分别由Wiener(1930)和Khinchin(1934)建立.（2）定义之二式(4.1.1)中的通常根据自相关函数定义计算,设为复信号,则(4.1.3)对于平稳随机信号,由于服从各态历经定理,因而集合平均可以用时间平均代替.因此上式可表示为(4.1.4)将上式代入式(4.1.1),可得的另一定义:(4.1.5)通常可能获得的观测数据为有限个,因此只能计算自相关函数和功率谱的估值,这就是功率谱估计问题.功率谱估计的两种方法:功率谱的估计方法有多种,一般可分为两大类:经典谱估计和现代谱估计.)(mRxx)(nx)]()([)(mnxnxEmRxxNNnNxxmnxnxNmR)()(121lim)()(jxxeP2)(121lim)(NNnnjNjxxenxNEeP)(nx§4.2经典谱估计经典谱估计建立在传统的傅里叶变换基础之上.通常有两种方法:BT法和周期图法.4.2.1BT法1958年由布莱克曼(R.Blackman)和杜基(J.Tukey)提出.BT法先是根据有限个观测数据估计自相关函数,然后再利用傅氏变换得到功率谱估值.通常采用有偏自相关函数估计(方差较小),公式为,(4.2.1)称为取样自相关函数,这种估计是渐近一致估计.对上式进行傅氏变换,得BT法的功率谱估计值为,(4.2.2)为了进一步减小估计方差,还可对上式进行加窗处理.可以证明,按上式估计是非一致估计,因而不是好的估计.1||0)()(1)(ˆmNnxxmnxnxNmR1||Nm)(ˆmRxxmjMMmxxjxxemReP)(ˆ)(ˆ1||NM取二人姓氏的第1个字母,故称BT法.又称自相关法.该式在无偏估计的基础上通过三角加窗得到的.4.2.2周期图法1.周期图定义周期图法是通过对随机数据直接进行傅里叶变换的一种功率谱估计方法.1965年Cooley和Tukey提出FFT以后,使DFT的运算量显著降低,从而为周期图法的应用创造了条件.根据定义,的自相关函数可表示为与的卷积,即据此,式(4.2.1)可写成,)(nx)(nx)(nxnnxxmxmxnmxnxnmxnxmR)()()]([)()()()()]()([1)(ˆmxmxNmRxx1||Nm上式的傅氏变换用符号表示,则(4.2.3)其中,是的点DFT:(4.2.4)由于是周期函数,所以用式(4.2.3)估计,称为“周期图法”,可用FFT实现.实际上,根据功率谱的第二种定义式:若忽略求统计平均的运算,并设观测数据为:,所得结果与式(4.2.3)相同.)(NI2)(1)]()([1)(jjjNeXNeXeXNI)(jeX)(nxN10)()(NnnjjenxeX)(jeX)(jxxeP2)(121lim)(NNnnjNjxxenxNEeP)(nx10Nn2.周期图与BT法的等价性BT估计算法可归结为:由随机序列,首先求得有偏自相关函数估计,然后求的傅氏变换:周期图法可归结为:由随机序列,首先求的傅氏变换再求周期图即取多个样本傅氏变换幅值模的平方,再求平均值作为功率谱估计.利用以上关系式可以证明,周期图就是BT法中有偏自相关估值所得的功率谱估计.)}1(,),1(),0({Nxxx1||0)()(1)(ˆmNnxxnmxnxNmR1||Nm)(ˆmRxx)](ˆ[DFT)(ˆmRePxxjxx)}1(,),1(),0({Nxxx)(nx)]([DFT)(nxeXj2)(1)(jNeXNI)(ˆmRxx●证明如下:由周期图表达式令,即,则由于1010)(1010210)()(1)()(1)(1)(NkNnnkjNnnjNkkjNnnjNenxkxNenxekxNenxNInkmnmk1)1(||10)()(1)(NNmmjmNnNenmxnxNI||101ˆ()()()NmxxnRmxnxnmN令m=k-n:(1)k=0m=-nn=0N-1m=-(N-1)0(2)k=N-1m=N-1-nn=0N-1m=0(N-1)综合(1),(2)可得m的数值范围是:-(N-1)≤m≤N-1有偏估计所以周期图法不需要估计自相关函数,可以利用FFT进行计算,因此比BT法简单,在出现FFT以后,应用较为广泛.1BT(1)ˆˆ()[()]()NjjmNxxxxmNIPeRme即同期图与采用有偏自相关估计的BT法是一致的.4.2.3经典谱估计的问题与改进1.主要问题经典谱估计的基础是傅氏变换.因傅里叶变换域为无穷大,而观测数据是有限的,对观测不到的数据实际上都被强制地当作0处理,这相当于无限长样本用矩形窗加以截断.由于窗口以外的数据仍有很强的相关性,因此,用有限长样本序列估计出来的功率谱必然存在很大的偏差.有限长序列的傅氏变换,等于无限长序列与矩形窗函数各自傅氏变换的卷积.矩形窗函数的频谱是sinc函数,卷积结果带来的影响是:（1）sinc函数主瓣的影响——分辨率下降sinc函数的主瓣将引起信号功率谱向邻近频域扩展(即“频谱泄漏”),使谱的主瓣展宽,谱的分辨率下降.（2）sinc函数旁瓣的影响——谱间干扰谱间干扰的后果是,强信号的功率谱旁瓣会影响弱信号功率谱的检测，造成谱的模糊甚至畸变,严重时可能掩盖弱信号,或者把旁瓣误认为是信号,以至造成假信号.由于在采用自相关函数的有偏估计时,周期图法即是BT法,因此,当序列为有限长时,BT法与周期图法及其改进方法,都具有相同的问题.2.改进方法（1）平均周期图法对随机信号进行组独立观测,一次观测得到个记录数据.共有个观测数据.其中第组的周期图为(4.2.5)的功率谱是个独立周期图的平均,即(4.2.6)改进效果:的估计方差是的.（2）窗函数法选择适当的窗函数对周期图进行平滑.设点序列周期图为,窗口的谱函数为,加窗后的功率谱是与的卷积,即(4.2.7))(nxLMLMNi210)(1)(MnnjiienxMI)(nxL11ˆ()()LjxxiiPeIL)(ˆjxxeP)(iIL1L个周期图之间近似互不相关.N)(nx)(NI)(nw)(jeW)(NI)(jeWdeWIePjNjxx)()(21)(ˆ)(其中具有低通特性.上式等效于对周期图进行频域修正,即让周期图通过一个线性非频变系统,以滤除周期图中的快变成分.（3）修正的周期图求平均法这种方法与平均周期图法类似.首先把长度为点的信号分为成段,每段含个观测数据,对每一段周期图进行加窗修正,得,(4.2.8)式中(4.2.9)然后对个近似互不相关的求平均,得到信号的功率谱为(4.2.10))(jeWN)(nxLM210)()(1)(MnnjiienwnxUILi,,2,1102)(1MnnwMUL)(iILiijxxILeP1)(1)(ˆ[评价]（1）以上三种改进方法都是以降低分辨率为代价,换取估计方差的减少,没有从根本上解决提高分辨率的问题.（2）现代谱估计法很好地解决了经典谱估计上述固有问题,特点是:●舍弃了传统的“傅里叶变换”这一处理环节;●排除了有限数据长度带来的真实谱与数据截取窗函数频谱之间的卷积效应.§4.3谱估计的参数模型法4.3.1现代谱估计基础现代谱估计以信号模型为基础.平稳随机序列的信号模型如图4.3.1所示.图中,输入白噪声的均值为0,方差为,的功率谱由下式决定:(4.3.1)由于模型传输函数与模型参数有关,因此,信号的功率谱估计问题,变成了由观测数据估计信号模型的参数问题.按照上述思路,功率谱估计可分下列三个步骤:（1）选择合适的信号模型;（2）根据的有限个观测数据(或它的有限个自相关函数),估计模型参数（3）计算模型输出功率谱.)(nw)(nx图4.3.1平稳随机序列的信号模型)()()(zAzBzH)(nx)(nw2w)(nx22|)(|)(jwjxxeHeP)(nx4.3.2AR模型谱估计根据第三章(§3.3)的讨论,已知AR()模型的差分方程:(4.3.2)当为实序列时,已导出AR(p)模型的Yule-Walker方程(即AR模型参数与信号自相关函数之间的关系)为(4.3.3)上式的矩阵形式为(4.3.4))(nw)(nxAR模型)(1)(zAzHppkkpnwknxanx1,)()()()(nxpiwxxipixxixxmiRamimRamR1210,)(1,)()(001)0()1()()1()0()1()()1()0(21wpxxxxxxxxxxxxxxxxxxaaRpRpRpRRRpRRR令(4.3.5)可进一步将式(4.3.4)表示成如下的矢量形式:(4.3.6)式中,为维全零列矢量.由上可见,只要已知(或估计出)个自相关函数值,即可由该方程式解出个模型参数.AR模型的系统函数为(4.3.7))1()1()1()0()1()()1()0()1()()1()0(ppxxxxxxxxxxxxxxxxxxpxxRpRpRpRRRpRRRRT21],,,,1[paaaa12)1(pwpxx0aR1p01p1p1p},,,,{221wpaaapkkkzazAzH1AR11)(1)(采用Levinson-durbin递推算法求解.由式(4.3.1)和式(4.3.7),可得以AR模型为基础的谱估计计算式:(4.3.8)212221|)(|)(pkkjkwjwjxxeaeAeP4.3.3ARMA模型谱估计在ARMA(p,q)模型中,其输入与输出满足如下差分方程:(4.3.9)其中.现用乘以上式两边,并取数学期望,得(4.3.10)式中(4.3.11)(4.3.12)设ARMA模型为稳定因果系统,其冲激响应函数为,则有(4.3.13)本节与下节内容参考王世一(北理工).数字信号处理.p.308311)(nw)(nxpkqkkkknwbknxanx10)()()(10b)(mnxpkqkwxkxxkxxkmRbkmRamR10)()()()]()([)(mnxnxEmRxx)]()([)(mnxnwEmRwx)(nh0)()()()()(llnwlhnwnhnx代入式(4.3.12)得(4.3.14)注:因此有(4.3.15)将上式代入式(4.3.12),得(4.3.16)020)()()]()([)()(lwlwxlmlhlmnwnwElhmRmlmlmhlmlhmRwlwwx),(,0)()()(20200,)(,0)(2mmmhmRwwx为满足系统的因果性,应再对