周期圖譜估計(jì)——頻譜泄露及改進(jìn)

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

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

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

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

5、心，窗函數(shù)的連續(xù)頻譜將偏離主瓣的中心，頻率偏移量對(duì)應(yīng)著信號(hào)頻率和頻率偏移量對(duì)應(yīng)著信號(hào)頻率和FFT 頻率頻率分辨率的差異，這個(gè)偏移導(dǎo)致了頻譜中分辨率的差異，這個(gè)偏移導(dǎo)致了頻譜中出現(xiàn)旁瓣，所以，窗函數(shù)的旁瓣特性直出現(xiàn)旁瓣，所以，窗函數(shù)的旁瓣特性直接影響著各頻譜分量向相鄰頻譜的泄漏接影響著各頻譜分量向相鄰頻譜的泄漏寬度。寬度。7簡(jiǎn)介n因此，周期圖法功率譜估計(jì)中頻譜泄露因此，周期圖法功率譜估計(jì)中頻譜泄露的改進(jìn)辦法既是選擇合適的窗函數(shù)，盡的改進(jìn)辦法既是選擇合適的窗函數(shù)，盡量減少頻譜泄露，改善功率譜估計(jì)的性量減少頻譜泄露，改善功率譜估計(jì)的性能。能。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ù)譜估計(jì)9周期圖法n定義定義n在實(shí)際應(yīng)用中，周期圖譜估計(jì)的計(jì)算式為：在實(shí)際應(yīng)用中，周期圖譜估計(jì)的計(jì)算式為：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估計(jì)性能估計(jì)性能n漸近無(wú)偏性漸近無(wú)偏性n周期圖的方差周期圖的方差(當(dāng)當(dāng)N較大時(shí)較大時(shí))n周期圖譜估計(jì)的方差不隨數(shù)據(jù)記錄長(zhǎng)度周期圖譜估計(jì)的方差不隨數(shù)據(jù)記錄長(zhǎng)度 N 的增大而減小，的增大而減小，而是近似于功率譜理論值的平方。而是近似于功率譜理論值的平方。n周期圖譜估計(jì)不是一致，這是一個(gè)令人失望的結(jié)果。周期圖譜估計(jì)不是一致，這是一個(gè)令人失望的結(jié)果。 ()lim( )( )PerxxNE SS2()2( ) 0( )2( ) 0,PerxxxSVar SS;第4章 非參數(shù)譜

9、估計(jì)11周期圖法改進(jìn)多個(gè)周期圖求平均n把數(shù)據(jù)記錄切分為把數(shù)據(jù)記錄切分為K個(gè)分段，分別求周期個(gè)分段，分別求周期圖，然后求平均。圖，然后求平均。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多個(gè)周期圖求平均n設(shè)設(shè)K 個(gè)數(shù)據(jù)分段之間互不相關(guān)，則個(gè)數(shù)據(jù)分段之間互不相關(guān)，則n為一個(gè)漸近無(wú)偏估計(jì)和一致性估計(jì)。為一個(gè)漸近無(wú)偏估計(jì)和一致性估計(jì)。n如果如果N 固定，且固定

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

11、穿函數(shù)對(duì)估計(jì)性能的影響穿函數(shù)對(duì)估計(jì)性能的影響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