第13章_經(jīng)典功率譜估計(jì)_

1、13.1 概述13.2 自相關(guān)函數(shù)的估計(jì)13.3 經(jīng)典譜估計(jì)的基本方法13.4 經(jīng)典譜估計(jì)的質(zhì)量13.5 經(jīng)典譜估計(jì)的改進(jìn)13.6 經(jīng)典譜估計(jì)算法比較13.7 短時(shí)傅里葉變換請(qǐng)抓住并搞清楚如下四個(gè)問題： 功率譜為什么要估計(jì)？ 如何估計(jì)？ 如何評(píng)價(jià)估計(jì)質(zhì)量？ 如不理想，如何改進(jìn)？13.1 13.1 概述概述221()lim( )21()lim21Mjj nXMnMjMP eEx n eMX eEM隨機(jī)信號(hào)的單個(gè)樣本求均值運(yùn)算求極限運(yùn)算()( )jj mXxmPer m e集總平均兩者等效平穩(wěn)信號(hào)( )X n單一樣本( , )( )x n ix n可將 看作能量信號(hào)，因此，可對(duì)它作傅立葉變換，并得

2、到功率譜：21()( )21Mjj nMMnMPexn eM( )Mxn問題 ： 的功率譜 和單個(gè)樣本的功率譜 有何關(guān)系？和整個(gè)隨機(jī)信號(hào)的功率譜 有何關(guān)系？( )Mxn()jxP e()jMPe()jXPe()jxP e()jXPe( )Mxn截短()jMPe1. 求極限：21()lim()lim( )21Mjjj nxMMMMnMP ePexn eM()()jjXxPeE P e2. 求均值：21()lim( )21Mjj nXMMnMPeExn eM單一樣本的功率譜不能收斂到所有樣本的功率譜，因此必須有求均值運(yùn)算，此即如下定義的來歷：各態(tài)遍歷信號(hào)也是如此。2()()1()lim( )211

3、lim( )( )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證明了兩個(gè)公式等效。所以自相關(guān)函數(shù)是集總自相關(guān)。證明：功率譜的兩個(gè)定義都要求：樣本無窮多，時(shí)間無限長(zhǎng)，即需要集總平均。實(shí)際工作中，我們往往能得到的是： 1. 單一的樣本； 2. 單一樣本的有限長(zhǎng)數(shù)據(jù)；問題：如何用這單一樣本的有限長(zhǎng)數(shù)據(jù)去估 計(jì)原隨機(jī)信號(hào)真實(shí)

4、的自相關(guān)函數(shù)和功 率譜目的：自身估計(jì)的需要； 功率譜估計(jì)的需要*( )( )()xr mE Xn X nm集總自相關(guān)*1( )lim( ) ()21NxNnNr mx n x nmN時(shí)間自相關(guān)定義：101 ( )( ) ()Nmxnr mx n x nmN 實(shí)際求出的自相關(guān)函數(shù)近似質(zhì)量如何Estimation Estimate Estimator（估計(jì)子）估計(jì)方法：從估計(jì)方法上看，實(shí)際上是把隨機(jī)信號(hào)“視為”單樣本有限長(zhǎng)的確定性信號(hào)。問題是：偏差自相關(guān)函數(shù)估計(jì)的質(zhì)量：101 ( )( ) ()Nmxnr mx n x nmN 估計(jì)方法單個(gè)樣本1. 偏差bir ( ) ( )( )r mE r m

5、r 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漸近無偏估計(jì)對(duì)固定的N，此結(jié)論給出了m的選取原則在數(shù)據(jù)上加矩形窗，長(zhǎng)度為 N ，該矩形窗函數(shù)的自相關(guān)函數(shù)正是三角窗！注意矩形窗加在數(shù)據(jù)上，三角窗加在相關(guān)函數(shù)上，體現(xiàn)在估計(jì)的自相關(guān)函數(shù)的均值上。( )(1) (1)Nmw mNNN,三角窗:(3) ( )( )

6、 ( )E r mw m r m那兒來的三角窗？方差222var ( ) ( ) ( ) ( ) ( )r mEr mE r mE rmE r m2.方差來自定義包含兩項(xiàng)22)()( mrNmNmrE前面結(jié)果11220021 ( )( ) ()( ) ()1( ) ( ) () ()NmNmnknkE rmEx n x nmx k x kmNE x n x k x nm x kmN 四階統(tǒng)計(jì)量！22( ) ( ) () ()()( )() ()E x n x k x nm x kmrnkrmr nkm r knm由：1(1)21var( )1( )() ()NmiNmmir mNNrir im

7、 r im 最后導(dǎo)出,var ( )0Nr m有：,bia ( )0,var ( )0Nr mNr m 漸近一致估計(jì)零均值高斯分布3.自相關(guān)函數(shù)的計(jì)算) 1(,),1 (),0(Nxxx已知單個(gè)樣本的 N 點(diǎn)數(shù)據(jù)估計(jì) ( )xr m兩個(gè)方法兩個(gè)方法: :（1 1） 直接按定義：直接按定義：( )() ( ),(1) (1)xxxr mrmr m mNN 利用 101 ( )( ) ()Nmxnr mx n x nmN 最大長(zhǎng)度（2） 利用FFT：Step1: 將 補(bǔ) 個(gè)零得 ;( )NxnN2( )NxnStep2: 對(duì) 做FFT,得 ；2( )Nxn2( )NXk22( )NXkStep3:

8、 對(duì) 求幅平方，得 ；2( )NXkStep4: 由 得 ， 對(duì)其作IFFT,得得 。221( )NXkN22( )NXk)(0mr思考： 和 有何關(guān)系)(0mr ( )xr mmNnNNmnxnxmNmr10)()(1)( 自相關(guān)函數(shù)的另一個(gè)估計(jì)方法（估計(jì)子）： ( )xr m( )xr m很容易證明： 是 的無偏估計(jì)，但方差性能不好。在一些譜估計(jì)的方法中，有時(shí)用到該公式。要求：很好掌握自相關(guān)函數(shù)的估計(jì)方法及估計(jì)性質(zhì)。13.3 經(jīng)典譜估計(jì)問題的提出：對(duì)隨機(jī)信號(hào) ，我們往往只能得到它的： 1. 單一的樣本 ；并且僅是 2. 單一樣本的有限長(zhǎng)數(shù)據(jù)； 如何用這 N 數(shù)據(jù)去估計(jì)原隨機(jī)信號(hào)真實(shí)的功率譜

9、( )X n( , )( )x n ix n) 1(,),1 (),0(Nxxx()jxP e1.周期圖（Periodogram）法：1201( )( ),( )( )Nj nNPERNnXx n ePXN經(jīng)典譜估計(jì)中有兩個(gè)基本的方法：1201( )( ),( )( )NnkNNPERNnXkx n WPkXkN思路：對(duì) 做DTFT(DFT),得到頻譜；對(duì)該頻譜求幅平方，再除以N，即得到“周期圖”功率譜，以此作為對(duì)真譜的估計(jì)。( )Nxn2.自相關(guān)（Blackman-Tukey BT法）法:101 ( )( ) ()Nmxnr mx n x nmN Step1( )( ),1Mj mBTxmM

10、Pr m eMNStep2因?yàn)橄纫烙?jì)自相關(guān)函數(shù)，所以又稱間接法。與此相對(duì)應(yīng)，周期圖法又稱直接法。3.直接法和間接法的關(guān)系：需要考慮兩種情況： （一）1MN （二）1MN ( ):,xr mmMM ( ):0,1,1x nnN數(shù)據(jù)的范圍自相關(guān)函數(shù)的范圍 （一）1MN比較用兩種方法的估計(jì)出的離散譜：2(1)22(1)221200( )( ) ( ),0,1,21NjmkNNBTxmNNjmkNmPkrm er m ekN2N 點(diǎn)的譜，把所能估計(jì)出的自相關(guān)函數(shù)都使用上了，而估計(jì)自相關(guān)函數(shù)時(shí)，把 N 點(diǎn)數(shù)據(jù)也全都使用上了。2221( )( )NPERNPkXkN21( )( )PERNPkXkN對(duì)

11、補(bǔ)N 個(gè)零，做DFT,得到( )NxnIFFT )(0mr221220( ),0,1,21NjmkNNBTkPk emN結(jié)論： 在 時(shí)，直接法和間接法估計(jì)的結(jié)果是一樣的。1MN使用間接法時(shí)，往往取 ，這時(shí)二者是不一樣的。因此，直接法可看作是間接法的特例。1MN 不補(bǔ)零，思考：( )Nxn即：21( )( )PERNPkXkNN點(diǎn)離散譜如何和 相等？2( )NBTPk2(1)22(1)( )( )0,1,21NjmkNNBTxmNPkrm ekNN點(diǎn)離散譜2100( )2Re( )(0)0,1,1NjmkNNBTmPkr m erkN （二）1MN( )( ) ( ),( ):Mrmr m v

12、mv mMM( )( )BTPERPP所以：( ):v mMM加在自相關(guān)函數(shù)上。目的是將其截短。第二次加窗。相當(dāng)于只用了部分自相關(guān)函數(shù)1(1)( )( ) ( ) ( )( )* ( )Nj mBTmNMj mMPERmMPr m v m erm ePV直接法和間接法之間的關(guān)系13.4經(jīng)典功率譜估計(jì)的質(zhì)量也分兩種情況討論1MN1MN主要考察的是均值方差無偏估計(jì)一致估計(jì)( )( )( )BTPERPPP（一）1MN周期圖和自相關(guān)法是等效的，統(tǒng)一考慮bia ( ) ( )( )PE PP1. 偏差估計(jì)值的均值1(1)( ) ( )Nj mmNr m w m e自相關(guān)函數(shù)估計(jì)的性質(zhì)1(1)1(1)(

13、 )( )( ) ( )PERBTNj mmNNj mmNE PE PEr m eE r m e201( )( )*( )( )*( )PERE PPWPDN于是有：( ):( )PX n的真實(shí)功率譜；00( ):( )Ddn0( ),0,1dnnN的頻譜譜；( ):( )Ww n的頻譜譜；0( )( )* ()w nd ndn三角窗；注意： 三角窗頻譜恒為正20bia ( )( )*( )( )1( )*|( )|( )PPWPPDPN最后有：201( )( )*( )PERE PPDNbia ( ) ( )( )PE PP由于如何理解這一結(jié)果lim( )lim( )()BTPERNNE P

14、E PPlim bia ( )0NP所以：201bia ( )( )*|( )|( )PPDPN周期圖和自相關(guān)法都是漸近無偏估計(jì)1MN020,( )( )1( )|( )|( )NDWDN 因?yàn)椋?Var ( ) ( ) ( ) PEPE P2. 方差又遇到四階矩問題，直接求解困難。（1）假定 是高斯零均值的隨機(jī)過程；思路：( )X n（2）求 在 處的協(xié)方差 ：12, P( )定義：121122covPPEPE PPE P(), ()()()()()有關(guān)方差公式的推導(dǎo)不作要求。主要是掌握結(jié)論，并用來說明問題。（3）令 ，則1211cov (),()var ( )PPP 1211221212c

15、ovPPEPE PPE PE PPE PE P(), ()()()()()() ()()()求解的關(guān)鍵2200var( )1( )()()( )2PPDDdE PN推導(dǎo)的結(jié)果：方差（1） 時(shí)2001( )()()02PDDdNN 22lim Var ( ) ( ) ( )0NPE PP經(jīng)典功率譜估計(jì)不是一致估計(jì)0( )D2B2B2001( )()()02PDDdN解釋：22BB00()()0DD0()D 0()D122010220102cov(), ()1( )()()21( )()()2PPPDDdNPDDdN推導(dǎo)的結(jié)果：協(xié)方差 假定 在主瓣外為零； 那么，在頻率范圍 內(nèi)：12cov (),

16、 ()0PP12|BB0( )D0( )D有（2）若 的主瓣寬度為 ；0( )D2B2B0201()()0DD01()D12 12B02()D12B01()D1202()D0102()()0DD在 12|B處，12cov (), ()0PP說明：隨機(jī)變量 在 處不相關(guān);12(), ()PP12, 原因：功率譜的定義中即要求極限，又要求均 值；而實(shí)際的估計(jì)方法，僅靠單次實(shí)現(xiàn) 的有限長(zhǎng)，無極限、又無均值運(yùn)算，因 此產(chǎn)生上述問題。設(shè)想：增大數(shù)據(jù)長(zhǎng)度，效果如何后果：使估計(jì)出的譜曲線起伏加??； 增大， 的主瓣（ ）將變窄，因此，引起不相關(guān)的區(qū)域進(jìn)一步增多，從而引起譜曲線的更加起伏，實(shí)際上是方差變大。N0

17、( )D分辨率和方差（體現(xiàn)在曲線起伏上），是經(jīng)典譜估計(jì)中的一對(duì)矛盾。BNN通常，增加 ，會(huì)提高譜的分辨率，對(duì)經(jīng)典譜估計(jì)來說，增加 固然會(huì)有利于提高分辨率，但譜曲線的起伏令使用者難以接受，這是經(jīng)典譜估計(jì)的一個(gè)致命缺點(diǎn)。對(duì)白噪聲在不同長(zhǎng)度情況下估計(jì)出的譜曲線：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經(jīng)典譜估計(jì)質(zhì)量的討論：(二)1 NM)(*)()(VPPPERBT( )v m：加在估計(jì)的自相關(guān)函數(shù)上 ，1,|NMMm ( )xr m周

18、期圖譜估計(jì)和自相關(guān)法的譜估計(jì) 不再一樣！( )( )* ( )( )*( )* ( )BTPERE PE PVPWV1. 偏差bia( )( )( )BTBTPE PP誰的主瓣比較寬)(),(VW( )( )NW ( )( )* ( )BTE PPVbia( )( )( )0limBTNPPP假定1： 是慢變譜，在 的主瓣內(nèi)近似為一個(gè)常數(shù)( )P( )V( )P假定2 1( )2(0)1Vdv窗函數(shù)的一般要求( )( )* ( )1( ) ()2limBTNE PPVPVd1( )( )2PVd也是漸近無偏估計(jì)！2.方差：考慮特殊情況， 為白噪序列，其功率譜應(yīng)為常數(shù)，即),()(nunx242

19、24sin ()var( )1sin ( )var( )limPERPERNNPNP 時(shí)對(duì)白噪聲功率譜估計(jì)的方差2( )P1MN42var( ) ( )2BTPVdN 時(shí)對(duì)白噪聲功率譜估計(jì)的方差1MN ：方差改進(jìn)之比 22var( )1 ( )2var( )1( )1BTrPERMmMPKVdNPvmN兩種情況下估計(jì)的方差之比：rKNMKr8) 12(3/2,3/8/4,3/16rrMNKMNK若則若則取哈明窗： 1. 在 加上 后，估計(jì)的譜 的偏差劣于 M=N1 時(shí)估計(jì)的譜，而方差優(yōu)于 M=N1 時(shí)估計(jì)的譜；)( mr)(mv( )BTP （2）在在 的范圍上， 因?yàn)锽變大， 不相關(guān)的點(diǎn)變少

20、。12|B12(),(),BTBTPP2. 上加窗 以后，估計(jì)譜 方差的改進(jìn)體現(xiàn)在兩個(gè)方面：( )BTP)( mr)(mv)(*)()(VPPPERBT（1）估計(jì)的譜曲線變得平滑些原主瓣寬，取決于2( ),WkN現(xiàn)主瓣寬，取決于2( ),VkM3. 方差的減小是以犧牲分辨率為代價(jià)的！若分辨率能滿足要求，這樣做是有意義的。即：既保證了分辨率，又使估計(jì)出的譜較為平滑。13.5 直接法估計(jì)的改進(jìn)任務(wù)：改進(jìn) 對(duì) 估計(jì)的性能；21( )|( )|PERNPXN( )P目標(biāo)：主要是改進(jìn)方差的性能方法：平滑與平均；)(*)()(VPPPERBT用 對(duì) 的加窗來實(shí)現(xiàn))(mv)( mr1. 平滑（Smoothi

21、ng)平滑1222:;:/LiXXXXXXL 理論依據(jù): L個(gè)獨(dú)立同分布隨即變量和的分布，方差減小 倍，即：L將一個(gè)較長(zhǎng)的信號(hào)分成若干段，對(duì)每一段求功率譜，每一段的功率譜都是隨機(jī)變量，然后平均之。類似相干平均，用以彌補(bǔ)經(jīng)典譜估計(jì)中缺少的求均值運(yùn)算。注意：信號(hào)應(yīng)是平穩(wěn)的，且每一段的統(tǒng)計(jì)特性基本一樣。2. 平均（Average）（1） Bartlett平均將 分成 段，每段 點(diǎn)，即)(nxLMNL M)(nx1( )x n( )ix n2( )xn3( )x n( )Lxn1( )d n?ML思考：如何確定或者2101( )( ),1,MiijnPERnPx n eiLM每一段譜平均后譜21101

22、1( )( )( )LMiij nPERPERinPPx n eLML平均后估計(jì)出的功率譜的性能如何？在數(shù)據(jù)上加了數(shù)據(jù)窗 寬度是),(1ndM結(jié)果，在自相關(guān)函數(shù)上引入了窗函數(shù))(1mw?ML思考：如何確定或者： 的自相關(guān)； 類似 引入的)(mw1( )( )*( )PERE PPW)(0nd)(1mw)(1nd統(tǒng)計(jì)性能分析： （1）偏差增大，分辨率進(jìn)一步下降； （2）方差減小，但到不了 倍L（2） Welch平均特點(diǎn)：交疊分段01NM.)(nx1( )x n( )ix n2( )xn)(2nd若重疊一半，段數(shù) /2/2NMLM變大1221012201( )|( )( )|1( )LMijnPE

23、RinMnPx n d n eMULUdnM ：不一定是矩形窗，如Hamming窗)(2nd歸一化因子，保證無偏估計(jì)Welch 平均是常用的經(jīng)典譜估計(jì)方法，MATLAB中有相應(yīng)的命令( )( )PERE PPWelch平均法的方差比Barttlett方法有明顯的減小，而偏差幾乎沒有減小3. Nottall 法：平滑與平均相結(jié)合假定1： 是慢變譜，在 的主瓣內(nèi)近似為一個(gè)常數(shù)）( )P2( )W( )P假定221( )12Wd22( )( )*( )( )( )2PERE PPWPWd三種改進(jìn)方法:如何比較每一個(gè)估計(jì)方法性能的好壞？想辦法產(chǎn)生一個(gè)已知功率譜的標(biāo)準(zhǔn)數(shù)據(jù):H(z)(2nu01. 02u

24、H(z)(1nu01. 02u)(1nv)(2nv13.6 總結(jié)與比較請(qǐng)掌握如下的方法：白噪聲1白噪聲2兩個(gè)輸出都是隨機(jī)信號(hào) 由自己指定( )H z122421( )( )( )()21jj kyukky nv njv nP eh e令：則：構(gòu)成一復(fù)信號(hào)得到 的功率譜；( )y n在 的基礎(chǔ)上再加上四個(gè)復(fù)正弦，歸一化頻率分別是：( )y n12340.15,0.16,0.252,0.16ffff 412)()(knfjkkeAnynx即調(diào)整 ，可以得到不同的信噪比，本例取這樣， 的真實(shí)功率譜可得到，并可畫出。我們可以此作為比較各種算法的依據(jù)。kA123464dB,54dB2dB,30dBfff

25、f處為處為處為處為)(nx實(shí)際工作中，對(duì)信號(hào) 總?cè)∮邢揲L(zhǎng)，如 ，由這128點(diǎn)去“求”功率譜，得到的當(dāng)然是估計(jì)值。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）真實(shí)譜；（）真實(shí)譜；（b）周期圖；（）周期圖；（c）Welch平均，四平均，四段，無迭合，段，無迭合，Hamming窗；（窗；（d）同）同c, 但迭合但迭合16點(diǎn)點(diǎn)-0.

26、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經(jīng)典功率譜估計(jì)的特點(diǎn)經(jīng)典功率譜估計(jì)的特點(diǎn)：1. 物理概念明確，可用FFT快速算法。所以 是大眾化的譜估計(jì)方法；2. 對(duì)周期圖，分辨率受到 的限制； 對(duì)自相關(guān)法，分辨率受到 的限制；2 / N2 /M3. 方差性能不好，不是一致估計(jì)， N 增 大時(shí)譜曲線反而起伏加??；4. 改進(jìn)方法是“平滑”與“平均”，改進(jìn)的目 的是減小方差，但犧牲了分辨率； 5. 注意窗函數(shù)的作用與影響：加在數(shù)據(jù)上的窗函數(shù)：012,d d d產(chǎn)生加在自相關(guān)函數(shù)上

27、的延遲窗：12( ),( ),( )w m w m w m各個(gè)窗函數(shù)的作用及影響是什么？13.7 短時(shí)傅里葉變換平穩(wěn)信號(hào)：均值、方差及均方都不隨時(shí)間變化，自相關(guān)函數(shù)僅和兩個(gè)觀察時(shí)間的差有關(guān)，和觀察的具體位置無關(guān) ；非平穩(wěn)信號(hào)：均值、方差都隨時(shí)間變化，自相關(guān)函數(shù)也和觀察的時(shí)間位置有關(guān)，信號(hào)的頻率也隨時(shí)間而變化，如語音、腦電及其他含有較多突變分量的信號(hào) 。其一階、二階統(tǒng)計(jì)量和功率譜的估計(jì)顯然不能簡(jiǎn)單地使用平穩(wěn)信號(hào)的估計(jì)方法，必須考慮其時(shí)變因素。方法：分段，每一小段可看作是平穩(wěn)的。*,*STFT ( , )( )( )( )()( ), ()xtjjtxgdxgt edxgt e )()(2RLtx1|)(|g其STFT定義為:并且窗函數(shù)應(yīng)取對(duì)稱函數(shù)。 式中 ,( )()jtggt e)()(3tgxx()0)()(1tgx)()(2tgx1t2t3tFTFTFT01t2t3tt22STFT ( ,)( ) ()( ,)jxxtxgt edS t 譜圖是恒正的，且是實(shí)的。 “譜圖（spectrogram）”1|)(|g由于xxEdtdtS),(所以譜圖是信號(hào)能量的分布。考慮 是隨機(jī)信號(hào)的一個(gè)樣本，譜圖可實(shí)現(xiàn)信號(hào)功率譜的估計(jì)。注意，它們是 的函數(shù)( )x t( ,)t 將信號(hào) 變換為一個(gè)二維函數(shù) 的方法稱為信號(hào)