周期圖譜估計——頻譜泄露及改進

1、周期圖法功率譜估計周期圖法功率譜估計頻譜泄漏及改進頻譜泄漏及改進2簡介n功率譜估計一般分成兩大類：功率譜估計一般分成兩大類：n經(jīng)典譜估計，也稱為非參數(shù)譜估計。經(jīng)典譜估計，也稱為非參數(shù)譜估計。n現(xiàn)代譜估計，也稱為參數(shù)譜估計?，F(xiàn)代譜估計，也稱為參數(shù)譜估計。n經(jīng)典譜估計是建立在傳統(tǒng)的傅立葉變換基礎之經(jīng)典譜估計是建立在傳統(tǒng)的傅立葉變換基礎之上的。經(jīng)典譜估計又可以分為兩種方法：上的。經(jīng)典譜估計又可以分為兩種方法：n相關圖法。相關圖法。n1958年，年，Blackman和和Tukey首先提出相關圖法。首先提出相關圖法。n周期圖法。周期圖法。n1898年，年，Schuster在尋找太陽黑子數(shù)據(jù)中隱藏的周期性

2、的在尋找太陽黑子數(shù)據(jù)中隱藏的周期性的研究工作中，提出了周期圖法，但直到研究工作中，提出了周期圖法，但直到1965年提出年提出FFT以以后，周期圖法才受到人們的重視。后，周期圖法才受到人們的重視。 3簡介n相關圖法中，先由有限個觀測數(shù)據(jù)估計相關圖法中，先由有限個觀測數(shù)據(jù)估計自相關函數(shù)，然后計算自相關序列的傅自相關函數(shù)，然后計算自相關序列的傅里葉變換得到功率譜。里葉變換得到功率譜。n周期圖法直接對觀測數(shù)據(jù)進行傅里葉變周期圖法直接對觀測數(shù)據(jù)進行傅里葉變換，取模的平方，再除以換，取模的平方，再除以N得到功率譜。得到功率譜。周期圖法比相關圖法簡單，可用周期圖法比相關圖法簡單，可用FFT進行進行計算，得到

3、了廣泛的應用。計算，得到了廣泛的應用。4簡介n在周期圖譜估計中，我們?nèi)∫欢斡邢揲L在周期圖譜估計中，我們?nèi)∫欢斡邢揲L的數(shù)據(jù)進行傅里葉變換，相當于對原始的數(shù)據(jù)進行傅里葉變換，相當于對原始信號作了矩形窗運算。輸入數(shù)據(jù)通過一信號作了矩形窗運算。輸入數(shù)據(jù)通過一個窗函數(shù)相當于原始數(shù)據(jù)的頻譜與窗函個窗函數(shù)相當于原始數(shù)據(jù)的頻譜與窗函數(shù)頻譜的卷積。窗函數(shù)的頻譜由一個主數(shù)頻譜的卷積。窗函數(shù)的頻譜由一個主瓣和幾個旁瓣組成，主瓣以時域信號的瓣和幾個旁瓣組成，主瓣以時域信號的每個頻率成份為中心。旁瓣在主瓣的兩每個頻率成份為中心。旁瓣在主瓣的兩側(cè)以一定的間隔衰減至零。側(cè)以一定的間隔衰減至零。5簡介nFFT 產(chǎn)生離散的頻譜

4、，出現(xiàn)在產(chǎn)生離散的頻譜，出現(xiàn)在FFT 每個每個譜線的是在每個譜線上的連續(xù)卷積頻譜。譜線的是在每個譜線上的連續(xù)卷積頻譜。如果原始信號的頻譜成份與如果原始信號的頻譜成份與FFT 中的譜中的譜線完全一致，這種情況下采樣數(shù)據(jù)的長線完全一致，這種情況下采樣數(shù)據(jù)的長度為信號周期的整數(shù)倍，頻譜中只有主度為信號周期的整數(shù)倍，頻譜中只有主瓣。沒有出現(xiàn)旁瓣的原因是旁瓣正處在瓣。沒有出現(xiàn)旁瓣的原因是旁瓣正處在窗函數(shù)主瓣兩側(cè)采樣頻率間隔處的零分窗函數(shù)主瓣兩側(cè)采樣頻率間隔處的零分量點。量點。第4章 非參數(shù)譜估計6簡介n如果時間序列的長度不是周期的整數(shù)倍，如果時間序列的長度不是周期的整數(shù)倍，窗函數(shù)的連續(xù)頻譜將偏離主瓣的中

5、心，窗函數(shù)的連續(xù)頻譜將偏離主瓣的中心，頻率偏移量對應著信號頻率和頻率偏移量對應著信號頻率和FFT 頻率頻率分辨率的差異，這個偏移導致了頻譜中分辨率的差異，這個偏移導致了頻譜中出現(xiàn)旁瓣，所以，窗函數(shù)的旁瓣特性直出現(xiàn)旁瓣，所以，窗函數(shù)的旁瓣特性直接影響著各頻譜分量向相鄰頻譜的泄漏接影響著各頻譜分量向相鄰頻譜的泄漏寬度。寬度。7簡介n因此，周期圖法功率譜估計中頻譜泄露因此，周期圖法功率譜估計中頻譜泄露的改進辦法既是選擇合適的窗函數(shù)，盡的改進辦法既是選擇合適的窗函數(shù)，盡量減少頻譜泄露，改善功率譜估計的性量減少頻譜泄露，改善功率譜估計的性能。能。8常用窗函數(shù)特性020406000.20.40.60.81

6、nRectangle窗 函 數(shù)00.51-300-200-1000w/pi幅度(dB)Rectangle窗 函 數(shù) 的 頻 譜020406000.20.40.60.81nTriangular窗 函 數(shù)00.51-300-200-1000w/pi幅度(dB)Triangular窗 函 數(shù) 的 頻 譜020406000.20.40.60.81nHanning窗 函 數(shù)00.51-300-200-1000w/pi幅度(dB)Hanning窗 函 數(shù) 的 頻 譜020406000.20.40.60.81nHanming窗 函 數(shù)00.51-300-200-1000w/pi幅度(dB)Hanming窗

7、函 數(shù) 的 頻 譜020406000.20.40.60.81nblackman窗 函 數(shù)00.51-300-200-1000w/pi幅度(dB)blackman窗 函 數(shù) 的 頻 譜020406000.20.40.60.81nKaiser窗 函 數(shù)00.51-300-200-1000w/pi幅度(dB)Kaiser窗 函 數(shù) 的 頻 譜第4章 非參數(shù)譜估計9周期圖法n定義定義n在實際應用中，周期圖譜估計的計算式為：在實際應用中，周期圖譜估計的計算式為：n其中，其中，w(n)為窗函數(shù)。為窗函數(shù)。212()011( )()( )NPerjj nxnSX ex n eNN221()021( )( )

8、 ( )1 ( ) ( ) , 0,1,.,1NjknPWNxnSkx n w n eNDFT x n w nkNN10周期圖法n估計性能估計性能n漸近無偏性漸近無偏性n周期圖的方差周期圖的方差(當當N較大時較大時)n周期圖譜估計的方差不隨數(shù)據(jù)記錄長度周期圖譜估計的方差不隨數(shù)據(jù)記錄長度 N 的增大而減小，的增大而減小，而是近似于功率譜理論值的平方。而是近似于功率譜理論值的平方。n周期圖譜估計不是一致，這是一個令人失望的結(jié)果。周期圖譜估計不是一致，這是一個令人失望的結(jié)果。 ()lim( )( )PerxxNE SS2()2( ) 0( )2( ) 0,PerxxxSVar SS;第4章 非參數(shù)譜

9、估計11周期圖法改進多個周期圖求平均n把數(shù)據(jù)記錄切分為把數(shù)據(jù)記錄切分為K個分段，分別求周期個分段，分別求周期圖，然后求平均。圖，然后求平均。nBartlett方法：方法：D=L。Welch方法：方法： D=L/2( )() ( )01,01ix nx iDn w nnLiK 212,011( )()( )Ljj nx iiinSX ex n eLL112(),0011( )( )()KKPAjxx iiiiSSX eKKL12多個周期圖求平均n設設K 個數(shù)據(jù)分段之間互不相關，則個數(shù)據(jù)分段之間互不相關，則n為一個漸近無偏估計和一致性估計。為一個漸近無偏估計和一致性估計。n如果如果N 固定，且固定

10、，且 N = KL，為了降低方差而增，為了降低方差而增加加K，會導致，會導致L的減少，也就是分辨率的下降。的減少，也就是分辨率的下降。n在實際應用中，用在實際應用中，用DFT/FFT計算計算DTFT，則，則()()211( )( )()PAPerjxxxVar SVar SSeKK112()(),0011( )()()( )2/, ( ) ( ), 0,1,.,KKPAPAxxkx ikiiikFFTiiFFTSkSSX kKKLk NX kDFT x nkN第4章 非參數(shù)譜估計13多個周期圖求平均14多個周期圖求平均n用用Welch方法估計隨機過程的功率譜及演示方法估計隨機過程的功率譜及演示

11、穿函數(shù)對估計性能的影響穿函數(shù)對估計性能的影響123( )cos(0.35)2cos(0.4)0.5cos(0.8)( )x nnnnv n1500.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Rectangle Welch Estimate,N=512,K=4,D=102,L=20400.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Triangular Welch Estimate,N=512,K=4,D=1

12、02,L=20400.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Hanning Welch Estimate,N=512,K=4,D=102,L=20400.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Hamming Welch Estimate,N=512,K=4,D=102,L=20400.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency

13、 (radian/pi)Blackman Welch Estimate,N=512,K=4,D=102,L=20400.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Kaiser Welch Estimate,N=512,K=4,D=102,L=204,Beta=201600.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Rectangle Welch Estimate,N=4096,K=4,D=819,L=16

14、3800.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Triangular Welch Estimate,N=4096,K=4,D=819,L=163800.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Hanning Welch Estimate,N=4096,K=4,D=819,L=163800.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Hamming Welch Estimate,N=4096,K=4,D=819,L=163800.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Blackman Welch Estimate,N=4096,K=4,D=819,L=163800.10.20.30.40.50.60.70.80.91-1001020Magnitude (dB)Frequency (radian/pi)Kaiser Welch