四種時頻分析方法的頻率分辨率研究

摘 要：時頻分析是分析時變譜的有力工具，其頻率分辨率是值得研究的關(guān)鍵問題之一。通過對基于短時傅里葉變換、小波變換、Choi-Williams分布和Hilbert-Huang變換的四種時頻分析方法的頻率分辨率的實驗比較， 說明Hilbert-Huang變換法具有最高的頻率分辨率，其次是Choi-Williams分布，小波緊隨其后，最差是短時傅里葉變換。



關(guān)鍵詞：短時傅里葉變換；小波變換；Choi-Williams分布；Hilbert-Huang變換；時頻分析；頻率分辨率

中圖分類號：TN911.72 文獻(xiàn)標(biāo)識碼：A

文章編號：1004-373X(2008)11-037-03

Experimental Research on Frequency Resolution of Time-Frequency

Distribution Using Four Methods



YANG Shenshen1，LI Ming1，GU Xuekang2

(1.School of Information Science and Technology，East China Normal University，Shanghai，200241，China;

2.China Ship Scientific Research Center，Wuxi，214082，China)



Abstract：Time-Frequency (TF) is a powerful tool for the analysis of time-varying spectrum.One of key issues in this regard is frequency resolution.Experimentally compared the frequency resolution of four types of TF analysis methods，namely，Short Time Fourier Transform (STFT)，Wavelet Transform (WT)，Choi-William Distribution (CWD) and Hilbert-Huang Transform (HHT).The results suggest that HHT has the highest frequency resolution，CWD is better than WT and STFT is the worst.



Keywords：short-time Fourier transform;wavelet transform;Choi-William distribution;Hilbert-Huang transform;time-frequency analysis;frequency resolution

時頻分析(TF)的基本思想是設(shè)計時間和頻率的聯(lián)合函數(shù)，用它同時描述信號在不同時間和頻率的能量密度或強度是分析非平穩(wěn)信號的有力工具。信號特征的提取與TF分布的頻率分辨率有很大關(guān)系。在頻率分辨率的研究中，取頻率分量較少且頻率之差盡量小的正弦疊加信號是較簡便的實驗做法^[1]。本文用4個頻率分量(f1，f2，f3，f4)的正弦疊加信號作為待分析信號，把STFT，WT，CWD和HHT四種時頻分析拿出來做比較，展現(xiàn)四者在頻率分辨率上的差異。

1 算法原理

1.1 Short Time Fourier變換(STFT)

為了研究信號在時間t的特性，一個自然的想法是加強在那個時間的信號，壓縮在其他時間的信號。這可以通過用中心在t的窗函數(shù)h(t)乘信號來實現(xiàn)，產(chǎn)生的被改變的信號為:



st(τ)=s(τ)h(τ－t)

(1)



這里，改變的信號是兩個時間即t(固定時間)和τ(執(zhí)行時間)的函數(shù)。窗函數(shù)決定留下的信號圍繞著時間t大體上不變，而在遠(yuǎn)離t的地方信號被壓縮了許多倍。也就是:



st(τ)s(τ) 對于τ接近t

0 對于τ遠(yuǎn)離t

(2)



因為改變的信號加強了圍繞著時間t的信號，所以，它的Fourier變換反映的是在時刻t的一個小的時間間隔內(nèi)的頻率分布，即:



St(ω)=12π∫e－jωτst(τ)dτ

=12π∫e－jωτst(τ)h(τ－t)dτ

(3)



因此，在時間t的能量密度頻譜是:



Psp(t，ω)=|St(ω)|2

=12π∫e－jωτ s(τ)h(τ－t)2

(4)



對于每一個不同的時間，都可以得到一個不同頻譜，這些頻譜的總體就構(gòu)成了一個時頻分布，即頻譜圖^[2]。

窗函數(shù)對STFT的性能有很大的影響：為了提高時間分辨率，需要選擇的窗函數(shù)h(t)盡可能短；但為了頻率分辨率，則要求選擇的窗函數(shù)h(t)時間寬度盡可能長，因此與時間分辨率的提高相矛盾。另一方面，對于要分析的非平穩(wěn)信號而言，也許某一小時段上是以高頻信息為主，我們希望用短時間窗進(jìn)行分析；而在某一場時間段上是一些低頻信息，我們希望用一個長時間窗進(jìn)行分析。因此對一個時變的非平穩(wěn)信號，利用STFT難以找到一個合適的時間窗口來適應(yīng)于不同的時間段，這也是它的最大不足之處^[3]。

1.2 Wavelet變換(WT)

小波分析^[4，5]本質(zhì)上是一種時間尺度分析。它克服了STFT的時頻窗不能改變形狀的缺陷，而可以有效地聚焦信號的瞬時結(jié)構(gòu)。

小波的一般定義如下:



W(a，b;X，Ψ)=|a|－12∫∞－∞X(t)Ψ*(t－ba)dt

(5)



其中，Ψ*(#8226;)為小波基函數(shù)，a為尺度因子，b表示與坐標(biāo)原點的平移。方程(5)的一個直觀的物理解釋是：W(a，b;X，Ψ)是X在時間t = b，尺度為a的能量。

在本文，我們使用Morlet小波計算信號x(t)的時頻分布。其計算表達(dá)式如下：



SCx(t，a;h)=|Tx(t，a;h)|2

=1|a|∫∞－∞x(s)h*s－tads2

(6)



其中a=v0v，v0是母小波函數(shù)h(t)的中心頻率^[6]。

1.3 Chio-Williamas分布(CWD)

20世紀(jì)60年代中期，Cohen發(fā)現(xiàn)眾多的時頻分布只是Wigner-Ville分布(WVD)的變形，可以用統(tǒng)一的形式表示，習(xí)慣稱之為Cohen類時頻分布，那么其表達(dá)式為^[2，7]：



Wx(t，ω)=∫∞－∞∫∞－∞∫∞－∞x(u+τ2)x(u－τ2)

φ(θ，τ)e－j(tτ+τω－uθ)dudθdτ

(7)



式中x(t)為解析信號，φ(θ，τ)為核函數(shù)。當(dāng)取φ(θ，τ)=1時，即表示W(wǎng)igner-Ville分布：



Wx(t，ω)=∫∞－∞x(t+τ2)x(t－τ2)e－jωπdτ

(8)



Choi和Williams^[7]提出了時頻分析核，其表現(xiàn)形式為：



φ(θ，τ)=e－θ2τ2/σ

(9)



Choi和Williams核函數(shù)與時間和頻率無關(guān)，是時延θ和頻延τ的函數(shù)，具有時頻移不變性。σ為衰減系數(shù)，它與交叉項的幅值成比例關(guān)系。

相應(yīng)的分布為:



Wx(t，ω)=∫∞－∞∫∞－∞∫∞－∞x(u+τ2)x(u+τ2)#8226;

e－θ2τ2/σ2e－j(tτ+τω－ωθ)dudθdτ

(10)



此分布稱為Choi-Williams分布。注意到當(dāng)σ→∞時，我們便得到WVD，因為此時核趨近于1。相反地，σ越小，相干項的衰減就越大。此分布具有以下數(shù)學(xué)性質(zhì)^[6]：

P1 y(t)=x(t－t0)Wy(t，ω)=Wx(t－t0，ω)

y(t)=x(t)ej2πω0tWy(t，ω)=Wx(t，ω－ω0)



時移和頻移不變性:

P2 y(t)=kx(kt)

k>0Wy(t，ω)=Wx(kt，ωk) 時頻伸縮性。

P3 fx(t)=∫+∞－∞ωWxa(t，ω)dω∫+∞－∞Wxa(t，ω)dω(xa是信號x對應(yīng)的解析信號)瞬時頻率。

P4 tx(ω)=∫+∞－∞tWxa(t，ω)dω∫+∞－∞Wxa(t，ω)dω群延遲。

我們在Matlab軟件中用時頻工具箱中tfrcw.m來計算。

1.4 Hilbert-Huang變換

非平穩(wěn)信號的主要特征是其頻率是時間的函數(shù)。但研究表明：STFT，WT，CWD這些時頻分析法的最終理論依據(jù)都是傅里葉變換，因而不可避免地會受到傅里葉變換分析非平穩(wěn)信號的缺陷的影響，如出現(xiàn)虛假頻率和多余信號分量等。N.E.Huang等人提出的HHT從根本上突破了傅里葉變換理論的限制，首次建立了一種基于瞬時頻率的、較系統(tǒng)和完整的信號分析方法^[8]。

文獻(xiàn)［9］提出了內(nèi)蘊模式函數(shù)的概念和EMD篩選算法。固有模態(tài)函數(shù)(IMF)必須滿足兩個條件：信號極值點的數(shù)量與過零點的數(shù)量必須相等，或最多相差一個；在任一時間點上，信號極大值定義的上包絡(luò)和極小值定義的下包絡(luò)的局部均值為零。

EMD算法：對任一實信號s(t)，首先確定出s(t)上的所有極大值點和極小值點，然后將所有的極大值點和所有的極小值點分別用一條光滑的曲線聯(lián)接起來，使兩條曲線間包含所有的信號數(shù)據(jù)。將這兩條曲線分別作為s(t)的上下包絡(luò)線，計算出它們的平均值曲線m(t)，s(t)與m(t)的差記作h(t)，則^[8]:



s(t)-m(t)=h(t)



把h(t)作為原信號，重復(fù)上面介紹的步驟，直到當(dāng)h(t)滿足一定的條件時，記:



c1(t)=h(t)



將c1(t)視為一個IMF，再作:



s(t)-c1(t)=r(t)



將r(t)視為新的s(t)，重復(fù)以上過程，一次得第二個IMF c2(t)，第三個IMF c3(t)，…。當(dāng)cn(t)或r(t)滿足給定的終止條件時，分解過程終止，得分解式:



s(t)=∑Nn=1cn(t)+rn(t)

(11)



其中，r(t)稱為殘余函數(shù)，代表信號的平均趨勢。

對式(11)中的每個IMF分別作Hilbert變換后得:



s(t)=Re∑ni=1ai(t)ejΦi(t)=Re∑ni=1ai(t)ej∫ωi(t)dt

(12)



這里省略了殘余函數(shù)r(t)，Re表示取實部。稱展開式(12)為Hilbert幅值譜，簡稱Hilbert譜，記作:



H(ω，t)=Re∑ni=1ai(t)ej∫ωi(t)dt

(13)



2 實驗結(jié)果

我們選用四個頻率的余弦信號進(jìn)行分析，如下：



x(t)=［cos(2πf1t)+cos(2πf2t)］［u(t)－u(t－100)］+

［cos(2πf3t)+cos(2πf4t)］［u(t－200)－u(t－300)］



其中f1，f2，f3，f4皆是頻率。令f1=0.05，f2=0.20，f3=0.25，f4=0.45，則:Δf1=0.15，Δf2=0.20。由此可見，兩組頻率之間的差是很微小的，通過觀察4種時頻方法的等高圖和立體圖，可以比較出它們頻率分辨率的優(yōu)劣。

首先來看HHT中各IMF分量及其余量的Hilbert譜：

圖1 信號x(t)HHT中各IMF分量及其余量的Hilbert譜

STFT的窗函數(shù)選的是長度為11的hamming窗，進(jìn)行WT分析時，利用時頻工具箱中的函數(shù)［tfr，t，f，wt］=tfrscalo(x，t，wave，fmin，fmax，N，trace)，其中wave取3；進(jìn)行CWD分析時，利用時頻工具箱中的函數(shù)［tfr，t，f］=tfrcw(x，t，N，G，H，signal，trace)，其中時間窗G取的是長度為15的hamming窗，時間窗h取的是長度為11的hamming窗，衰減系數(shù)σ即函數(shù)中的signal取為1。

由圖2可以看到STFT的頻率分辨率很糟糕，兩組頻率交織在一起；WT的情況好一些，第一組頻率f1，f2尚能分辨，但第二組頻率f3，f4就分不清楚了；CWD顯然比上面兩種情況要改善了許多，兩組頻率很清晰地顯示出來，只是Cohen類的分布都屬于二次型的，因而產(chǎn)生的相干項不可避免的干擾了時頻表示的可讀性；HHT頻率分辨率最優(yōu)。

圖2 信號x(t)的時頻等高圖

由立體圖圖3可以進(jìn)一步證實上面的結(jié)論。另一方面我們研究的信號幅值都是相等的，而WT這種方法幅值的變化過于凸顯，與事實不符，其他三種方法都表現(xiàn)良好。

圖3 信號的立體圖比較

3 結(jié) 語

本文對基于STFT，WT，CWD和HHT的四種時頻分布的頻率分辨率進(jìn)行了實驗比較。結(jié)果表明，STFT的頻率分辨率最低，WT較前者有了明顯改善，CWD表現(xiàn)更好， HHT的頻率分辨率最優(yōu)。在這里只對頻率分辨率做了實驗比較，時間分辨率也是一個很好的研究方向，但本文還未涉及，作者將在以后的工作中繼續(xù)對這方面的工作進(jìn)行深入研究。

參 考 文 獻(xiàn)

［1］Li Ming.Experimental Research for Processing of Choi-Williams Distribution and Bessel Distribution［A］.IEEE Conf.Proc.，1st Int.Conf.on Information，Communication Signal Processing，Singapore，1997，3:1 593-1 597.



［2］［美］ L.科恩.時-頻分析：理論與應(yīng)用［M］.白居憲，譯.西安:西安交通大學(xué)出版社，1998.

［3］葛哲學(xué)，陳仲生.Matlab時頻分析技術(shù)及其應(yīng)用［M］.北京:人民郵電出版社，2006.



［4］Mallat S.A Theory for Multi-Resolution Signal Decomposition Wavelet Representation ［J］.IEEE Trans.on Pattern Analysis and Machine Intell，1989，11(7):674-693.



［5］Coilman R R，Wickerhauser M W.Entropy-based Algorithms for Basis-Selection ［J］.IEEE Tran.Information Theory，1992，38(2):713-718.



［6］胡昌華，周濤，夏啟兵，等.基于MATLAB的系統(tǒng)分析與設(shè)計——時頻分析［M］.西安:西安電子科技大學(xué)出版社，2001.



［7］劉希強，沈萍，山長侖，等.數(shù)字化地震波形資料的時頻分析方法及應(yīng)用［J］.西北地震學(xué)報，2004，26(2):118-125.



［8］鐘佑明，秦樹人.HHT的理論依據(jù)探討——Hilbert變換的局部乘積定理［J］.振動與沖擊.2006，25(2):12-19.



［9］Huang N E，Zheng Shen，Long S R，et al.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear Non-Stationary Time Series Analysis ［J］.Proc RSOC London，1998，454:903-995.



［10］Li Ming，Jia W，Zhao W.Analysis of Pattern Similarity between Choi-Williams Kernel and Born-Jordan Kernel［A］.Proc.of IEEE ICICS99，Singapore，1999.



作者簡介 楊哂哂 女，1984年出生，碩士研究生，從事信號處理研究。

李 明 男，1955年出生，PHD，華東師范大學(xué)教授，博士生導(dǎo)師。主要研究信號處理，計算機(jī)科學(xué)等。

顧學(xué)康 男，博士，中國船舶科學(xué)研究中心研究員，博士生導(dǎo)師。

注：本文中所涉及到的圖表、注解、公式等內(nèi)容請以PDF格式閱讀原文。