Chapter+4中国科学院大学现代数字信号处理课程课件

1、4.1 引言 4.2 经典谱估计 4.3 现代谱估计中的参数建模 4.4 AR模型谱估计方法 4.5 最大熵谱估计方法4.6 特征分析法谱估计功率谱定义功率谱估计的方法功率谱估计的应用信号的功率谱和其自相关函数服从一对傅里叶变换关系mjmxxjxxmrPe )()e (dPmrnjjxxxxe )e (21)(经典谱估计方法间接方法：BT法直接方法：周期图法现代谱估计方法参数法：ARMA模型法（AR模型、MA模型、ARMA模型）非参数法：Pisarenko谐波分解法、MUSIC方法在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数)；常常利用功率谱估计来得到线性系统的参数估计；从

2、宽带噪声中检测窄带信号。 22)()(HPyyBT法周期图法改进的周期图法BT法是先估计自相关函数, 然后进行傅里叶变换得到功率谱。有偏自相关函数估计的误差相对较小，是一种渐近一致估计：1|0*)()(1)(mNnxxmnxnxNmr-BT(e)( )ejjmxxmPrmn 周期图法的定义如下: 212j-j011(e )( )eNnjxxnPx nX eNN观 测 数 据 x(n)FFT取 模的 平 方1/NPxx(ej )1. 1. 周期图与周期图与BTBT法的等价关系法的等价关系 2101111*()00001(e)( )e11( )e( )e( )( )eNjj nxxnNNNNj k

3、j njk nknknPx nNx kx nx k x nNN令 m=k-n, 即k=m+n，则 1*01j(1)1(1)BT(e )ee(e1( ) () ( )Nj mxxmNmnxNNj mmxNjPPx n x mnNrm 利用有偏自相关函数的BT法和周期图法是等价的。 2. 周期图法谱估计质量分析周期图法谱估计质量分析 1) 周期图的偏移 1j-j(1)1-j(1)-j ( )|(e)ee( )( )e)NmxxmNNmmNmxxBxmxxxE rmNmrE Pwm rmmN其它0|)(NmNmNmwB式中 上式在频域表示为： j-1(e)(e)ed2jjxxBxxE PWP式中 )

4、()e (mrFTPxxjxx2)2/sin()2/sin(1)()e (NNmwFTWBjB 周期图的统计平均值等于它的真值卷积三角谱窗函数，因此周期图是有偏估计，但当N时，wB(m)1，三角谱窗函数趋近于函数，周期图的统计平均值趋于它的真值，因此周期图属于渐近无偏估计。 2) 周期图的方差 用这种方法估计的功率谱在2x附近起伏很大，故周期图是非一致估计。 白噪声的周期图 Bartlett平均周期图法窗口处理法平均周期图Welch法（修正的周期图求平均法） 主要思想：主要思想：对序列x(n)进行L次独立观测或将其分成L段，计算每组观测数据的周期图，再将L个周期图加和后求平均。 LiixxIL

5、P1j)(1)e ( 假设随机信号x(n)的观测数据区间为：0nN-1，共进行了L次独立观测，得到L组记录数据，每一组记录数据用xi(n), i=1, 2, 3, ,L表示； 或对长为N的数据x(n)分成L段，每段有M个数据，N=LM，第i段数据表示为xi(n)= x(n+iM-M)。 第i组的周期图用下式表示： 210)(1)(MnnjiienxMIn 估计方法：估计方法： 将得到的L个周期图进行平均，作为信号x(n)的功率谱估计， 公式如下： LiixxILP1j)(1)e (n 估计效果分析：估计效果分析： 平均周期图的估计方差是周期图的方差的1/L，L越大方差越小，功率谱越平滑；相应的

6、，M越小，偏移越大，分辨率越低；估计的均方误差也减少； 以分辨率的降低换取了估计方差的减少，估计量的方差和分辨率是一对矛盾。平均周期图法 Pxx(ej)/032101230321012303210123Pxx(ej)/Pxx(ej)/(a)(b)(c)N256， L2N256， L4N256， L8主要思想：主要思想：用一适当的功率谱窗函数W(ej)与周期图进行卷积，来达到使周期图平滑的目的。 (e)( )(e)jjxxNPIW()0 jW e其中，()1(e)( )(e)( )(e)d2jjjxxNNPIWIW 式中 mxxNNmNmrIj -1)1(e )()(de)e(21)(jjnWn

7、w-(M-1)nM-1 n 估计方法：估计方法：那么 mxxMMmjxxemwmrPj1)1()()()e (mxxMMmxxmwmrEPEj -1)1(je )()( )e (又 1j-j(1)()( )( ) ( )eMmxxxxBmME Perm wm w m 偏移分析： n 估计效果分析：估计效果分析：| ( )( )( )( )xxxxxxBNmE rmrmrmwmN可得 周期图的窗函数法仍然是有偏估计, 其偏移和wB(m)、w(m)两个窗函数有关。 如果w(m)窗的宽度比较窄，M比N小得多，这样|m|p AR谱估计方法可归结为求解AR模型系数或线性预测器系数的问题。AR模型参数估计

8、方法：信号预测误差最小原则（或预测误差功率最小）自相关法（Levison递推法）Burg法协方差法修正协方差法（前后向线性预测最小二乘法）最大熵原则最大熵谱估计方法1、 自相关法自相关法列文森（列文森（Levinson）递推）递推 估计估计方法：方法：自相关法的出发点是选择AR模型的参数使预测误差功率最小；采用Levison-Durbin递推方法求解Yule-Walker方程得到AR模型参数。 21j -22j2je11| )e (|)e (PiiiwwxxaHP 预测误差功率为 22111| ( )|( )()ppinnie nx na x niNN假设信号x(n)的数据区间在0nN-1范围

9、，有P个预测系数，N个数据经过冲激响应为api(i=0,1, 2, , P)的滤波器， 输出预测误差e(n)的长度为N+P， 因此应用下式计算: 2110210| )()(|1| )(|1pipiPNnPNninxanxNneN0( )( )( )()pepiie nx nh na x ni预测误差功率最小，得到 )()2() 1 ()0()2() 1()2()0() 1 () 1() 1()0(21prrraaarprprprrrprrrxxxxxxppppxxxxxxxxxxxxxxxxxx0pkpkaa 采用Levinson-Durbin递推法求解Yule-Walker方程： 由k=1开

10、始递推，递推到k=p，依次得到a11,21，a21,a22, 22，ap1,ap2,app,2p。 AR模型的各个系数以及模型输入白噪声方差求出后， 信号功率谱用下式计算： 21j -22j2je11| )e (|)e (PiiiwwxxaHP利用列文森递推法计算功率谱的流程图 输入： x(n),n0,1,2,3,N1,2计算：)0(1),0(/ ) 1 (21 , 1211 , 1xxxxxxrarra)(mrxxp=22122*, 1, 12111, 111, 2 , 1 , 0)()(ppppkppppkppkppixxipxxppapkaaaaiprapra22ppiipiwxxaP1

11、j2je11)(e结束p=p+1NY性能分析：性能分析：该方法需要基于有限的观测数据估计自相关序列，当数据长度较短时，估计误差会比较大，AR参数的计算就会引入很大的误差。从而导致功率谱估计出现谱线分裂与谱峰频率偏移等现象。2、 伯格（伯格（Burg）递推法）递推法 n 估计方法：估计方法：直接由时间序列计算AR模型参数的方法，求前、后向预测误差平均功率最小时的反射系数kp，进而求AR模型参数ak和2w。20,ppkwpkak 设信号x(n)观测数据区间为：0nN-1，前向、后向预测误差功率分别用p,e和p,b表示，预测误差平均功率用p为 1122,11|( )|( )|1()2NNfbp ep

12、p bpn pn ppp ep benenNpNppkpkbppkpkfpkpnxapnxneknxanxne1*1)()()()()()(其中，前向、后向预测误差公式分别为 求预测误差平均功率p最小时的反射系数kp，令 0ppk1*11122112( )(1)(|( )|(1)| )Nfbppn ppNfbppn pen enkenen 基于反射系数kp，由Levinson-Durbin递推关系求AR模型参数ak和2w，进而求得功率谱Pxx 伯格递推法流程图 输入： x(n),n0,1,2,3,N1阶数IP1)(1)()()()(210000pnxNnxnenxneNnbfp = IP结束1

13、 pp12121*111)()() 1()(2NpnbpfpbpNpnfppnenenenekpwpppk2121pppipppippikaakaa*, 1, 11, 2, 1),() 1()(1, 2, 1),1()()(1*111NppnnekneneNppnneknenefppbpbpbpppfp2,;, 3 , 2 , 1:wippia输出21j2je1)(epiipiwxxaPNY性能分析：性能分析：该方法避免了采用有限数据估计自相关函数的计算，适合短序列参数估计，克服了L-D递推中的某些缺点，计算量小。但对正弦信号的谱估计，仍存在某些谱线分裂与频率偏移现象。3 3、 协方差法与修正

14、协方差法协方差法与修正协方差法（1 1）. . 协方差法协方差法n 估计方法：估计方法：利用使预测误差功率最小的方法求模型参数11212| )()(|1| )(|1NpnNpnpkpkknxanxpNnepN 该公式中使用的观测数据均已得到，不需要在数据两端补充零点， 因此比较自相关法去掉了加窗处理的不合理假设。0pkpkaan 性能分析：性能分析： 一些实验结果说明它的分辨率优于自相关法，另外对于纯正弦信号数据，可以有效地估计正弦信号的频率。 （2 2）. . 修正协方差法（前后向线性预测最小二乘法）：修正协方差法（前后向线性预测最小二乘法）：n 估计方法：估计方法：修正协方差法使用前向和后

15、向预测误差平均值最小的方法， 估计AR模型的参数，进而估计信号的功率谱。,0pp kp kaa,2p fp bp 前向和后向预测误差功率p,f、p,b分别用下式表示： 21,121*,011( )()1( )()pNp fpkn pkNppp bpknkx na x nkNpx na x nkNp 预测误差平均功率最小1*,1,1*,011( )()()( )()()01,2,3,pNpp kn pkp lNppp knkx nax nkx nlaNpx nax nkx nllp n 性能分析：性能分析： 该方法去掉了Burg法所用的Levinson的约束条件，估计得到的谱在谱线分裂和频率偏移

16、时较Burg法有较大改善。 自相关法可以用Levinson递推算法，运算量小，但分辨率受窗长度的限制；协方差法，去除了自相关法加窗处理的不合理假设，分辨率高，运算量较大；修正协方差法，分辨率高，在谱线分裂和偏移上较Burg法有较大改善，运算量大；Burg算法，可用改进的Levinson递推算法，分辨率高，但对正弦信号存在谱线分裂和偏移现象。 例例4.5.1 已知信号的四个观察数据为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)

17、()(21) 1()()(| )(|4111111111114011401111402aaaaanxneaneneanxanxnenennn (2) 协方差法: 714. 00)3()41 (4)24(20) 1()(32)()(32) 1()()(| )(|31| )(|111111111311131111131212aaaanxneaneneanxanxnenenepNnnnNpn4 4、 关于关于ARAR模型阶次的选择模型阶次的选择n对于白噪中的AR信号，其阶次的选择应折衷考虑。如选择AR模型，其阶次应加大，较低的阶次会使谱估计产生偏移， 降低分辨率。信噪比愈低，平滑作用愈严重，愈需要高

18、的阶次， 因此信噪比低应选高的阶次。阶次愈高，分辨率愈高；但阶次太高，会使估计误差加大，谱峰分裂。AR模型阶次太小时的平滑作用50.0030.0010.0010.0030.000.0000.400.50真实AR(4)的PSDAR(2)的PSD模型PSD/dB/n 最终预测误差（最终预测误差（FPEFPE）准则：）准则： 用下式计算后再取最小，估计模型阶次： kmin( )kNkFPE kNk式中, 表示k阶AR模型的白噪声方差（预测误差功率）的估计值，N是观测数据的长度，k是模型阶次。上式中, 随k增加而减少，但(N+k)/(N-k)却随k增加而加大。kkn 阿凯克信息论

19、准则：阿凯克信息论准则： kminAIC( )ln2kkNk 选择使上式最小的k值作为模型的阶次。但这种方法有过高估计模型阶次的趋势。为改进, 提出用k lnN代替公式中的2k的方法。 AIC和FPE估计阶次的性能是相近的， 但对于短记录数据，建议使用AIC准则，对于长记录（N）, 两种方法得到相同的模型阶次。 n 自回归传递函数准则（自回归传递函数准则（CATCAT）：）： k1111minCAT( )kiikkN式中 iiiNN同样, 模型的阶次取使上式最小的k值。最大熵谱估计与AR模型谱估计是一致的，采用最大熵原则外推自相关函数，得到自相关函数与模型参数服从Yule-Walker方程；通

20、过解Yule-Walker方程，解出模型参数，进而计算信号功率谱。1 1 利用最大熵的原则外推自相关函数利用最大熵的原则外推自相关函数 按照Shannon对熵的定义， 当随机变量X取离散值时，熵的定义为 iiippHln式中pi是出现状态i的概率。当X取连续值时，熵的定义为 )(lnd)(ln)(xpExxpxpH式中, p(x)是X的概率密度函数， 假设x(n)是零均值正态分布的平稳随机序列，它的N维高斯概率密度函数为 XNRXNRxxxpxxHxxNN12/12/21)(21exp)(det)2(),(式中 H21,NxxxX)0() 1()() 1()0() 1 ()() 1()0()(

21、xxxxxxxxxxxxxxxxxxxxrNrNrNrrrNrrrNR先讨论一维高斯分布的信号的熵，然后推广到N维。 222222222( )ln( )( )ln22( )ln2( )21ln22ln2lnln2xHp xp xdxp xdxx p x dxp x dxee 同理可求得N维高斯分布信号的熵为/21/2ln(2e)(det()NxxHRN式中det(Rxx(N)表示矩阵Rxx(N)的行列式，由上式表明为使熵最大，要求det(Rxx(N)最大。 maxmaxdet()xxHRN由观测序列x(n)估计自相关函数 ；用 代替Levison-Durbin算法里的rxx(m)，求解AR模型

22、参数 和 ；计算功率谱 ；(e)jxxP xxrmka2w21(e)1ejwxxNj kkkPa xxrm该方法特别适用于对多正弦加白噪声序列进行谱分析，可以得到比AR模型法更高的分辨率和更准确的频率估计，尤其在信噪比低时更有效。估计思想：将白噪声加正弦波作为一特殊的ARMA模型，用特征方程求它的参数，计算出正弦波的频率和功率，及噪声功率。1 1 正弦波用退化正弦波用退化ARAR模型模型表示表示令随机相位实正弦信号x(n)=sin(n+), 由下面三角恒等式：) 1(sincos2)2(sin)sin(nnn- )2() 1(cos2)(nxnxnx可得：x(n)=-a1x(n-1)-a2x(

23、n-2)122cos,1aa 其中，正弦波信号用一个特殊的AR（2）模型表示： P个正弦波组合的模型用2P阶差分方程描述 21( )()Pkkx na x nk Z变换 2111PkkkH za z 211PkkkA za z 对于多项式该方程有2P个根，这些根在z平面的单位圆上，它们是 iizjei=1, 2, 3, , 2P 式中的i即是正弦波的频率。2 2 白噪声白噪声中正弦波组合用一特殊的中正弦波组合用一特殊的ARMAARMA模型模型表示表示)()sin()()()(1nwnqnwnxnyPiiii 21( )( )( )()Pkkx ny nw nx na x nk 将代入21221

24、1( )- ( )-()()( )()( )()PiiPPiiiiy n w nay niw niy na y niw na w ni即上述差分方程是一特殊的ARMA（2p,2p）模型，其AR和MA参数相同，这种对称性是白噪声中正弦信号的一个特点。白噪声中正弦波组合的信号为 受到加性白噪声干扰的AR(p)模型等价于一个ARMA(p,p)模型。3 3 特征特征分解法谱估计分解法谱估计 YTA=WTA 式中 TT221T)2(,),1(),(, 1 )2(,),1(),(pnwnwnwWaaaApnynynyYpPiiPiiinwanwinyany2121)()()()(上式的矩阵形式 用向量Y左

25、乘上式, 并取数学期望, 得到 E EYY YY TA A=EYWYW TA 式中 yyyyyyyyyyyyyyyyyyyyRrprprprrrprrrYYE)0() 12()2() 12()0() 1 ()2() 1()0(TIWWEWWXEYWEw2TTT)( 特征向量A的求解：假定W是零均值白噪声信号，X是确定性信号。AARwyy2式中,w2是自相关矩阵Ryy的特征值; A是对应w2的特征向量。 对于白噪声中含有p个正弦信号的情况，当自相关矩阵Ryy的尺寸等于或大于(2p+1) (2p+1)的时候，方差w2是自相关矩阵Ryy的最小特征值。20yywRI A 求解各正弦波的频率：特征方程为

26、 0122221120PPPkkkzazazaza该方程有2P个根，这些根在z平面的单位圆上，它们是 iizjei=1, 2, 3, , 2P 式中的i即是正弦波的频率。4 Pisarenko4 Pisarenko谐波分解谐波分解mPmriPiiyycos)(1m0 式中， 22iiqP ，qi是频率为i的正弦波的幅度，Pi是其功率。对于白噪声中P个实正弦组合信号, 它的自相关函数为 PiiwyyPr12)0( 求各正弦波的幅值和功率： 12112212coscoscos12cos2cos2cos2coscoscosPyyyyPyyPPrPrPFPrPPPPPFP=r 由此，可以得到其矩阵表征形式如下： 求解噪声功率: piiyywPr12)0( 以上就是皮萨论科谱分解法的全过程。其中关键的一步是求自相关函数R Ryy的最小特征值及对应的特征向量。 Pisarenko方法的计算过程， 具体计算步骤归纳如下： (1) 由y(n)的观测数据估计出2p+1个相关函数ryy(m), m=0, 1,2, ,2p值。 (2) 求Ryy的最小特