1、11.1 概述11.2 自相关函数的估计11.3 经典谱估计的基本方法11.4 经典谱估计的质量11.5 经典谱估计的改进11.6 经典谱估计算法比较11.7 短时傅里叶变换请抓住并搞清楚如下四个问题： 功率谱为什么要估计？ 如何估计？ 如何评价估计质量？ 如不理想，如何改进？11.1 11.1 概述概述221()lim( )21()lim21Mjj nXMnMjMP eEx n eMX eEM随机信号的单个样本求均值运算求极限运算()( )jj mXxmP er m e集总平均两者等效平稳信号( )X n单一样本( , )( )x n ix n可将 看作能量信号，因此，可对它作傅立叶变换，并
2、得到功率谱：21()( )21Mjj nMMnMPexn eM( )Mxn问题 ： 的功率谱 和单个样本的功率谱 有何关系？和整个随机信号的功率谱 有何关系？( )Mxn()jxP e()jMPe()jXP e()jxP e()jXP e( )Mxn截短()jMPe1. 求极限：21()lim()lim( )21Mjjj nxMMMMnMP ePexn eM()()jjXxP eE P e2. 求均值：21()lim( )21Mjj nXMMnMP eExn eM单一样本的功率谱不能收敛到所有样本的功率谱，因此必须有求均值运算，此即如下定义的来历：各态遍历信号也是如此。2()()1()lim(
3、 )211lim( )( )211lim()21Mjj nXMMnMMMjm nMMMnM mMMMjm nxMnM mMP eExn eMExm xn eMr mn eM 22()(21) ( )MMMnM mMkMg mnMk g k 双求和变成单求和：22()lim(1) ( )21( )Mjj kXxMkMj kxkkPer k eMr k e证明了两个公式等效。所以自相关函数是集总自相关。证明：功率谱的两个定义都要求：样本无穷多，时间无限长，即需要集总平均。实际工作中，我们往往能得到的是： 1. 单一的样本； 2. 单一样本的有限长数据；问题：如何用这单一样本的有限长数据去估 计原随
4、机信号真实的自相关函数和功 率谱目的：自身估计的需要； 功率谱估计的需要*( )( )()xr mE Xn X nm集总自相关*1( )lim( ) ()21NxNnNr mx n x nmN时间自相关定义：101 ( )( ) ()Nmxnr mx n x nmN 实际求出的自相关函数近似质量如何Estimation Estimate Estimator（估计子）估计方法：从估计方法上看，实际上是把随机信号“视为”单样本有限长的确定性信号。问题是：偏差自相关函数估计的质量：101 ( )( ) ()Nmxnr mx n x nmN 估计方法单个样本1. 偏差bir ( ) ( )( )r m
5、E r mr m来自定义所有样本10101 ( )( ) ()1( ) ()( )NmnNmnE r mEx n x nmNE x n x nmNNmr mN (1),bia ( )0mNr m 固定，bia ( )( )mr mr mN所以：含义(2), ( )( )NmN r mr m给定，接近(3) ( )( ) ( )E r mw m r m渐近无偏估计对固定的N，此结论给出了m的选取原则在数据上加矩形窗，长度为 N ，该矩形窗函数的自相关函数正是三角窗！注意矩形窗加在数据上，三角窗加在相关函数上，体现在估计的自相关函数的均值上。( )(1) (1)Nmw mNNN,三角窗:(3) (
6、 )( ) ( )E r mw m r m那儿来的三角窗？方差222var ( ) ( ) ( ) ( ) ( )r mEr mE r mE rmE r m2.方差来自定义包含两项22)()( mrNmNmrE前面结果11220021 ( )( ) ()( ) ()1( ) ( ) () ()NmNmnknkE rmEx n x nmx k x kmNE x n x k x nm x kmN 四阶统计量！22( ) ( ) () ()()( )() ()E x n x k x nm x kmrnkrmr nkm r knm由：1(1)21var( )1( )() ()NmiNmmir mNNr
7、ir im r im 最后导出,var ( )0Nr m 有：,bia ( )0,var ( )0Nr mNr m 渐近一致估计零均值高斯分布3.自相关函数的计算) 1(,),1 (),0(Nxxx已知单个样本的 N 点数据估计 ( )xr m两个方法两个方法: :（1 1） 直接按定义：直接按定义：( )() ( ),(1) (1)xxxr mrmr m mNN 利用 101 ( )( ) ()Nmxnr mx n x nmN 最大长度（2） 利用FFT：Step1: 将 补 个零得 ;( )NxnN2( )NxnStep2: 对 做FFT,得 ；2( )Nxn2( )NXk22( )NXk
8、Step3: 对 求幅平方，得 ；2( )NXkStep4: 由 得 ， 对其作IFFT,得得 。221( )NXkN22( )NXk)(0mr思考： 和 有何关系)(0mr ( )xr mmNnNNmnxnxmNmr10)()(1)( 自相关函数的另一个估计方法（估计子）： ( )xr m( )xr m很容易证明： 是 的无偏估计，但方差性能不好。在一些谱估计的方法中，有时用到该公式。要求：很好掌握自相关函数的估计方法及估计性质。11.3 经典谱估计问题的提出：对随机信号 ，我们往往只能得到它的： 1. 单一的样本 ；并且仅是 2. 单一样本的有限长数据； 如何用这 N 数据去估计原随机信号
9、真实的功率谱( )X n( , )( )x n ix n) 1(,),1 (),0(Nxxx()jxP e1.周期图（Periodogram）法：1201( )( ),( )( )Nj nNPERNnXx n ePXN经典谱估计中有两个基本的方法：1201( )( ),( )( )NnkNNPERNnXkx n WPkXkN思路：对 做DTFT(DFT),得到频谱；对该频谱求幅平方，再除以N，即得到“周期图”功率谱，以此作为对真谱的估计。( )Nxn2.自相关（Blackman-Tukey BT法）法:101 ( )( ) ()Nmxnr mx n x nmN Step1( )( ),1Mj
10、mBTxmMPr m eMNStep2因为先要估计自相关函数，所以又称间接法。与此相对应，周期图法又称直接法。3.直接法和间接法的关系：需要考虑两种情况： （一）1MN （二）1MN ( ):,xr mmMM ( ):0,1,1x nnN数据的范围自相关函数的范围 （一）1MN比较用两种方法的估计出的离散谱：2(1)22(1)221200( )( ) ( ),0,1,21NjmkNNBTxmNNjmkNmPkrm er m ekN2N 点的谱，把所能估计出的自相关函数都使用上了，而估计自相关函数时，把 N 点数据也全都使用上了。2221( )( )NPERNPkXkN21( )( )PERNP
11、kXkN对 补N 个零，做DFT,得到( )NxnIFFT )(0mr221220( ),0,1,21NjmkNNBTkPk emN结论： 在 时，直接法和间接法估计的结果是一样的。1MN使用间接法时，往往取 ，这时二者是不一样的。因此，直接法可看作是间接法的特例。1MN 不补零，思考：( )Nxn即：21( )( )PERNPkXkNN点离散谱如何和 相等？2( )NBTPk2(1)22(1)( )( )0,1,21NjmkNNBTxmNPkrm ekNN点离散谱2100( )2Re( )(0)0,1,1NjmkNNBTmPkr m erkN （二）1MN( )( ) ( ),( ):Mrm
12、r m v mv mMM( )( )BTPERPP所以：( ):v mMM加在自相关函数上。目的是将其截短。第二次加窗。相当于只用了部分自相关函数1(1)( )( ) ( ) ( )( )* ( )Nj mBTmNMj mMPERmMPr m v m erm ePV直接法和间接法之间的关系11.4经典谱估计的质量也分两种情况讨论1MN1MN主要考察的是均值方差无偏估计一致估计( )( )( )BTPERPPP（一）1MN周期图和自相关法是等效的，统一考虑bia ( ) ( )( )PE PP1. 偏差估计值的均值1(1)( ) ( )Nj mmNr m w m e自相关函数估计的性质1(1)1
13、(1)( )( )( ) ( )PERBTNj mmNNj mmNE PE PEr m eE r m e201( )( )*( )( )*( )PERE PPWPDN于是有：( ):( )PX n的真实功率谱；00( ):( )Dd n0( ),0,1d nnN的频谱谱；( ):( )Ww n的频谱谱；0( )( )* ()w nd ndn三角窗；注意： 三角窗频谱恒为正20bia ( )( )*( )( )1( )*|( )|( )PPWPPDPN最后有：201( )( )*( )PERE PPDNbia ( ) ( )( )PE PP由于如何理解这一结果lim( )lim( )()BTPE
14、RNNE PE PPlim bia ( )0NP所以：201bia ( )( )*|( )|( )PPDPN周期图和自相关法都是渐近无偏估计1MN020,( )( )1( )|( )|( )NDWDN 因为：2Var ( ) ( ) ( ) PEPE P2. 方差又遇到四阶矩问题，直接求解困难。（1）假定 是高斯零均值的随机过程；思路：( )X n（2）求 在 处的协方差 ：12, P( )定义：121122covPPEPE PPE P(), ()()()()()有关方差公式的推导不作要求。主要是掌握结论，并用来说明问题。（3）令 ，则1211cov (), ()var ( )PPP 1211
15、221212covPPEPE PPE PE PPE PE P(), ()()()()()() ()()()求解的关键2200var( )1( )()()( )2PPDDdE PN推导的结果：方差（1） 时2001( )()()02PDDdNN 22lim Var ( ) ( ) ( )0NPE PP经典功率谱估计不是一致估计0( )D2B2B2001( )()()02PDDdN解释：22BB00()()0DD0()D 0()D122010220102cov(), ()1( )()()21( )()()2PPPDDdNPDDdN推导的结果：协方差 假定 在主瓣外为零； 那么，在频率范围 内：12
16、cov (),()0PP12|BB0( )D0( )D有（2）若 的主瓣宽度为 ；0( )D2B2B0201()()0DD01()D12 12B02()D12B01()D1202()D0102()()0DD在 12|B处，12cov (), ()0PP说明：随机变量 在 处不相关;12(), ()PP12, 原因：功率谱的定义中即要求极限，又要求均 值；而实际的估计方法，仅靠单次实现 的有限长，无极限、又无均值运算，因 此产生上述问题。设想：增大数据长度，效果如何后果：使估计出的谱曲线起伏加剧； 增大， 的主瓣（ ）将变窄，因此，引起不相关的区域进一步增多，从而引起谱曲线的更加起伏，实际上是方
17、差变大。N0( )D分辨率和方差（体现在曲线起伏上），是经典谱估计中的一对矛盾。BNN通常，增加 ，会提高谱的分辨率，对经典谱估计来说，增加 固然会有利于提高分辨率，但谱曲线的起伏令使用者难以接受，这是经典谱估计的一个致命缺点。对白噪声在不同长度情况下估计出的谱曲线：00.250.5-40-30-20-1001000.250.5-20-1001000.250.5-40-30-20-1001000.250.5-40-30-20-10010N16 N 32N64 N 128经典谱估计质量的讨论：(二)1 NM)(*)()(VPPPERBT( )v m：加在估计的自相关函数上 ，1,|NMMm (
18、)xr m周期图谱估计和自相关法的谱估计 不再一样！( )( )* ( )( )*( )* ( )BTPERE PE PVPWV1. 偏差bia( )( )( )BTBTPE PP谁的主瓣比较宽)(),(VW( )( )NW ( )( )* ( )BTE PPVbia( )( )( )0limBTNPPP假定1： 是慢变谱，在 的主瓣内近似为一个常数( )P( )V( )P假定2 1( )2(0)1Vdv窗函数的一般要求( )( )* ( )1( ) ()2limBTNE PPVPVd1( )( )2PVd也是渐近无偏估计！2.方差考虑特殊情况， 为白噪序列，其功率谱应为常数，即),()(nu
19、nx24224sin ()var( )1sin ( )var( )limPERPERNNPNP 时对白噪声功率谱估计的方差2( )P1MN42var( ) ( )2BTPVdN 时对白噪声功率谱估计的方差1MN ：方差改进之比 22var( )1 ( )2var( )1( )1BTrPERMmMPKVdNPvmN两种情况下估计的方差之比：rKNMKr8) 12(3/2,3/8/4,3/16rrMNKMNK若则若则取哈明窗： 1. 在 加上 后，估计的谱 的偏差劣于 M=N1 时估计的谱，而方差优于 M=N1 时估计的谱；)( mr)(mv( )BTP （2）在在 的范围上， 因为B变大， 不相
20、关的点变少。12|B12(),(),BTBTPP2. 上加窗 以后，估计谱 方差的改进体现在两个方面：( )BTP)( mr)(mv)(*)()(VPPPERBT（1）估计的谱曲线变得平滑些原主瓣宽，取决于2( ),WkN现主瓣宽，取决于2( ),VkM3. 方差的减小是以牺牲分辨率为代价的！若分辨率能满足要求，则这样做是有意义的，即既保证了分辨率，又使估计出的谱较为平滑。11.5 直接法估计的改进任务：改进 对 估计的性能；21( )|( )|PERNPXN( )P目标：主要是改进方差的性能方法：平滑与平均；)(*)()(VPPPERBT用 对 的加窗来实现)(mv)( mr1. 平滑（Sm
21、oothing)平滑1222:;:/LiXXXXXXL 理论依据: L个独立同分布随即变量和的分布，方差减小 倍，即：L将一个较长的信号分成若干段，对每一段求功率谱，每一段的功率谱都是随机变量，然后平均之。类似相干平均，用以弥补经典谱估计中缺少的求均值运算。注意：信号应是平稳的，且每一段的统计特性基本一样。2. 平均（Average）（1） Bartlett平均将 分成 段，每段 点，即)(nxLMNL M)(nx1( )x n( )ix n2( )xn3( )x n( )Lxn1( )d n?ML思考：如何确定或者2101( )( ),1,MiijnPERnPx n eiLM每一段谱平均后谱
22、211011( )( )( )LMiij nPERPERinPPx n eLML平均后估计出的功率谱的性能如何？在数据上加了数据窗 宽度是),(1ndM结果，在自相关函数上引入了窗函数)(1mw?ML思考：如何确定或者： 的自相关； 类似 引入的)(mw1( )( )*( )PERE PPW)(0nd)(1mw)(1nd统计性能分析： （1）偏差增大，分辨率进一步下降； （2）方差减小，但到不了 倍L（2） Welch平均特点：交叠分段01NM.)(nx1( )x n( )ix n2( )xn)(2nd若重叠一半，段数 /2/2NMLM变大1221012201( )|( )( )|1( )LM
23、ijnPERinMnPx n dn eMULUdnM ：不一定是矩形窗，如Hamming窗)(2nd归一化因子，保证无偏估计Welch 平均是常用的经典谱估计方法，MATLAB中有相应的命令( )( )PERE PPWelch平均法的方差比Barttlett方法有明显的减小，而偏差几乎没有减小3. Nottall 法：平滑与平均相结合假定1： 是慢变谱，在 的主瓣内近似为一个常数）( )P2( )W( )P假定221( )12Wd22( )( )*( )( )( )2PERE PPWPWd三种改进方法:如何比较每一个估计方法性能的好坏？想办法产生一个已知功率谱的标准数据:H(z)(2nu01.
24、 02uH(z)(1nu01. 02u)(1nv)(2nv11.6 总结与比较请掌握如下的方法：白噪声1白噪声2两个输出都是随机信号 由自己指定( )H z122421( )( )( )()21jj kyukky nv njv nP eh e令：则：构成一复信号得到 的功率谱；( )y n在 的基础上再加上四个复正弦，归一化频率分别是：( )y n12340.15,0.16,0.252,0.16ffff 412)()(knfjkkeAnynx即调整 ，可以得到不同的信噪比，本例取这样， 的真实功率谱可得到，并可画出。我们可以此作为比较各种算法的依据。kA123464dB,54dB2dB,30d
25、Bffff处为处为处为处为)(nx实际工作中，对信号 总取有限长，如 ，由这128点去“求”功率谱，得到的当然是估计值。127, 1 ,0，n)(nx421( )( )2()kxykkPPA -0.5-0.2500.250.5-60-40-200(a)-0.5-0.2500.250.5-60-40-200(b)-0.5-0.2500.250.5-60-40-200(c)-0.5-0.2500.250.5-60-40-200(d)（a）真实谱；（）真实谱；（b）周期图；（）周期图；（c）Welch平均，四平均，四段，无迭合，段，无迭合，Hamming窗；（窗；（d）同）同c, 但迭合但迭合16点
26、点-0.5-0.2500.250.5-60-40-200(e)-0.5-0.2500.250.5-60-40-200(f)（e）BT法，法，M32；（；（f）BT法，法，M16经典功率谱估计的特点经典功率谱估计的特点：1. 物理概念明确，可用FFT快速算法。所以 是大众化的谱估计方法；2. 对周期图，分辨率受到 的限制； 对自相关法，分辨率受到 的限制；2 / N2 /M3. 方差性能不好，不是一致估计， N 增 大时谱曲线反而起伏加剧；4. 改进方法是“平滑”与“平均”，改进的目 的是减小方差，但牺牲了分辨率； 5. 注意窗函数的作用与影响：加在数据上的窗函数：012,d d d产生加在自相
27、关函数上的延迟窗：12( ),( ),( )w m w m w m各个窗函数的作用及影响是什么？11.7 短时傅里叶变换平稳信号：均值、方差及均方都不随时间变化，自相关函数仅和两个观察时间的差有关，和观察的具体位置无关 ；非平稳信号：均值、方差都随时间变化，自相关函数也和观察的时间位置有关，信号的频率也随时间而变化，如语音、脑电及其他含有较多突变分量的信号 。其一阶、二阶统计量和功率谱的估计显然不能简单地使用平稳信号的估计方法，必须考虑其时变因素。方法：分段，每一小段可看作是平稳的。*,*STFT ( , )( )( )( )()( ), ()xtjjtxgdxgt edxgt e )()(2
28、RLtx1|)(|g其STFT定义为:并且窗函数应取对称函数。 式中 ,( )()jtggt e)()(3tgxx()0)()(1tgx)()(2tgx1t2t3tFTFTFT01t2t3tt22STFT ( ,)( ) ()( ,)jxxtxgt edS t 谱图是恒正的，且是实的。 “谱图（spectrogram）”1|)(|g由于xxEdtdtS),(所以谱图是信号能量的分布。考虑 是随机信号的一个样本，谱图可实现信号功率谱的估计。注意，它们是 的函数( )x t( ,)t 将信号 变换为一个二维函数 的方法称为信号的联合时频分析：( )x t( ,)t *,STFT ( , )( )(
29、 )xttxgd STFT2( ,) |STFT ( ,)|xxS tt 谱图*( ,)()()22jxW tx tx ted Wigner分布*()221( ,)()() ( , )2jtuxC tx ux ugedud d Cohen类分布例1信号x(n)由三个不同频率的正弦首尾相接所组成，即 1121232sin(),01( )sin(),1sin(),1nnNx nnNnNnNnN-101Real partSignal in time0797515951Linear scaleEnergy spectral density5010015020025030035000.10.20.30.
30、4|STFT|2, Lh=48, Nf=192, lin. scale, contour, Thld=5%Time sFrequency Hz21321NNN2( )exp()exp()x nj njn n例2线性频率调制信号（chirp）其频率与时间成正比非平稳 例2两个chirp信号，一个频率随时间增长，一个频率随时间减小，求它们和的谱图。与本章内容有关的MATLAB文件 本文件用Welch平均法估计一个信号的功率谱，其基本调用格式是： Px, F = pwelch(x, Nfft, Fs, window, Noverlap)式中 x 是随机信号，Fs是抽样频率，Nfft是对x作FFT时的长度，window是选用的窗函数, Noverlap是估计x的功率谱时每一段叠合的长度。缺省时，Nfft = 256 , noverlap = 0, window = Hanning(Nfft), Fs = 2。输出的Px是估计出的功率谱，按上述调用格式给出的是幅平方值，F是频率轴坐标。 1. pwelch.m 2. spectrum.m功能和pwelch.m类似，可用Welch平均法来估计一个信号的自功率谱，还可用于估计两个信号的互功率谱。3specgram.m 估计信号的谱图，但实际上估计的是其短时傅里叶变换。该文件主要针对非平稳信号，当然也可用于平稳信号，甚至确定性信号。
