李岳，韓賓，魯云

　?。ㄎ髂峡萍即髮W(xué) 信息工程學(xué)院，四川 綿陽(yáng) 621010）

       摘要：介紹了短時(shí)傅里葉變換、Cohen類時(shí)頻分布、小波變換、Hilbert-Huang變換四種典型的時(shí)頻分析方法，分析對(duì)比結(jié)果顯示了用Hilbert-Huang變換對(duì)聲音信號(hào)進(jìn)行時(shí)頻分析的優(yōu)越性，結(jié)合LabVIEW在數(shù)據(jù)采集和儀器控制領(lǐng)域的強(qiáng)大功能，提出以聲卡作為采集硬件、LabVIEW作為軟件編程、Hilbert-Huang變換作為時(shí)頻分析方法的一種聲音信號(hào)采集分析系統(tǒng)。實(shí)驗(yàn)結(jié)果表明，這種采集分析系統(tǒng)非常適合頻率在音頻范圍之內(nèi)（20 Hz~20 kHz）的聲音信號(hào)的時(shí)頻譜分析。

　　關(guān)鍵詞：聲卡；LabVIEW；時(shí)頻分析；EMD

0引言

　　聲音信號(hào)如同其他自然界的信號(hào)和人工合成的信號(hào)一樣，都是典型的非平穩(wěn)信號(hào)，其明顯的特點(diǎn)是信號(hào)是時(shí)變的且信號(hào)持續(xù)時(shí)間是有限的。想要從真實(shí)信號(hào)中提取出不同組成成分的時(shí)變信息，一般的做法是通過(guò)時(shí)頻分析方法將低維的一維時(shí)間信號(hào)映射到二維的時(shí)間-頻率函數(shù)空間，其目的是揭示信號(hào)包含了多少頻率分量以及各個(gè)頻率分量是如何隨時(shí)間變化的。

　　聲卡作為一個(gè)常見(jiàn)的計(jì)算機(jī)配置，其本身就是一個(gè)非常優(yōu)秀的數(shù)據(jù)采集系統(tǒng)，它搭載的A/D和D/A轉(zhuǎn)換器可以很方便地實(shí)現(xiàn)模擬信號(hào)和數(shù)字信號(hào)的相互轉(zhuǎn)換。如果被測(cè)對(duì)象的頻率在音頻范圍之內(nèi)（20 Hz~20 kHz），而且對(duì)采樣頻率沒(méi)有特別高的需求，則可以用計(jì)算機(jī)自帶的聲卡來(lái)構(gòu)建一個(gè)數(shù)據(jù)采集系統(tǒng)。而LabVIEW由于其直觀的編程方式和強(qiáng)大的功能函數(shù)庫(kù)等特點(diǎn)，已經(jīng)廣泛地被各界科研工作者和工程師們所采用，也被視為標(biāo)準(zhǔn)的儀器控制和數(shù)據(jù)采集軟件。通過(guò)聲卡和LabVIEW的聯(lián)合開(kāi)發(fā)，能夠?qū)β曇粜盘?hào)進(jìn)行數(shù)據(jù)采集、數(shù)據(jù)保存、波形顯示、信號(hào)分析等功能。

　　本文分析對(duì)比了四種典型的時(shí)頻分析方法：短時(shí)傅里葉變換、Cohen類時(shí)頻分布、小波變換、Hilbert-Huang變換，通過(guò)調(diào)頻合成信號(hào)的理論計(jì)算，闡述了四種方法的優(yōu)缺點(diǎn)。并且以計(jì)算機(jī)自帶的聲卡為硬件平臺(tái)，以LabVIEW作為軟件平臺(tái)，設(shè)計(jì)開(kāi)發(fā)了一種高效低成本的聲音信號(hào)的經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）時(shí)頻分析的系統(tǒng)。通過(guò)實(shí)驗(yàn)表明，把信號(hào)進(jìn)行EMD分解得到信號(hào)的時(shí)頻分布，對(duì)信號(hào)有很強(qiáng)的自適應(yīng)性，能夠很好地處理非線性、非平穩(wěn)的信號(hào)。

1時(shí)頻分析理論基礎(chǔ)

　　1.1短時(shí)傅里葉變換

　　短時(shí)傅里葉變換（ShortTime Fourier Transform，STFT）是1946年Gabor提出來(lái)的，基本思想是在信號(hào)進(jìn)行傅里葉變換之前先乘以一個(gè)時(shí)間有限的窗函數(shù)，并假定信號(hào)在窗內(nèi)是平穩(wěn)的，以此來(lái)確定窗內(nèi)存在的頻率成分，然后通過(guò)窗在時(shí)間軸上的移動(dòng)逐段分析信號(hào)，最后得到所需的時(shí)頻分布［1］。信號(hào)x(t)的短時(shí)傅里葉變換為：

　　式（1）中，x(t)為被分析的信號(hào)，g(t)為窗函數(shù)，當(dāng)窗函數(shù)g(t)取值為常數(shù)1時(shí)，短時(shí)傅里葉變換退化為傳統(tǒng)的傅里葉變換。

　　短時(shí)傅里葉的反變換為：

　　1.2Cohen類時(shí)頻分布

　　Cohen在20世紀(jì)60年代中期發(fā)現(xiàn)眾多的時(shí)頻分布只是WignerVille時(shí)頻分布的變形，可以用一個(gè)統(tǒng)一的形式表示，習(xí)慣稱之為Cohen類時(shí)頻分布［2］，其表達(dá)式為：

　　式（3）中,φ(τ,v)為核函數(shù)，當(dāng)其取值為1時(shí)，Cohen類時(shí)頻分布退化為WignerVille時(shí)頻分布。

　　常用的加上核函數(shù)之后的Cohen類時(shí)頻分布有BornJordan分布（BJD）、偽WignerVille時(shí)頻分布（PWD）、平滑WignerVille時(shí)頻分布（SPWD）。

　　1.3小波變換

　　小波變換是一個(gè)較新的應(yīng)用數(shù)學(xué)分支，在20世紀(jì)80年代后期工程應(yīng)用的需求促使其迅速發(fā)展起來(lái)。法國(guó)學(xué)者I.Daubechies和S.Mallat把這一理論應(yīng)用到工程應(yīng)用領(lǐng)域。其在信號(hào)處理領(lǐng)域起到非常重要的作用［3］。

　　信號(hào)x(t)的連續(xù)小波變換為：

　　式（4）中，φ(t)是母小波函數(shù)，φ(t)是φ(t)的共軛函數(shù)，τ是平移因子，a是尺度因子［3］。連續(xù)小波變換的逆變換為：

　　式（5）中，Cφ=∫+∞0|ψ(aw)|2ada<∞，是φ(t)需要滿足的容許條件。

　　1.4HilbertHuang變換

　　HilbertHuang變換是N.E.Huang等人在1998年首次提出的，它通過(guò)提取信號(hào)上下包絡(luò)的均值，反復(fù)篩選，自適應(yīng)地得到一系列的時(shí)域局部對(duì)稱且瞬時(shí)頻率具有明確物理意義的IMF（Intrinsic Mode Function）信號(hào)，能解決Hilbert變換不能處理多值頻率的信號(hào)的問(wèn)題［4］。

　　HilbertHuang變換包含2個(gè)主要步驟:

　　(1) 對(duì)原始數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）的預(yù)處理，把數(shù)據(jù)分解為滿足Hilbert變換所要求的n階固有模式函數(shù)(IMF)和殘余函數(shù)之和。

　　(2) 對(duì)每一階IMF進(jìn)行Hilbert變換，得到瞬時(shí)頻率，從而求得時(shí)頻分布［4］。

2幾種時(shí)頻分析比較

　　本文用一個(gè)正向和一個(gè)反向的高斯型調(diào)頻信號(hào)疊加成一個(gè)合成信號(hào)來(lái)檢驗(yàn)幾種時(shí)頻分析方法的應(yīng)用效果。在0~2 000 ms時(shí)間段由中心點(diǎn)1 000 ms、頻率為180 Hz的一個(gè)高斯型調(diào)頻信號(hào)疊加另一個(gè)中心點(diǎn)1 000 ms、頻率為50 Hz的高斯型調(diào)頻信號(hào)，兩個(gè)高斯型調(diào)頻信號(hào)的頻率都是在50 Hz~180 Hz之間變化，并且在287 ms~318 ms和1 679 ms~1 710 ms時(shí)間段疊加了兩個(gè)頻率為常數(shù)156 Hz的信號(hào)（圖1、圖2所示）。

　　分別對(duì)圖1中的信號(hào)進(jìn)行短時(shí)傅里葉變換、Cohen類時(shí)頻分布、小波變換和HilbertHuang變換，其時(shí)頻分析結(jié)果如圖3~圖6。圖3短時(shí)傅里葉變換圖3中，其窗函數(shù)選用的是65點(diǎn)的hanning窗，可以從時(shí)頻變換結(jié)果看出，短時(shí)傅里葉變換在一定程度上彌補(bǔ)了傳統(tǒng)的傅里葉變換不具備局部頻率分析的能力。但是，短時(shí)傅里葉變換還是通過(guò)滑動(dòng)時(shí)間窗來(lái)計(jì)算信號(hào)的頻譜，其必然會(huì)受到Heisenberg測(cè)不準(zhǔn)原理的約束［5］，也就是長(zhǎng)窗口有高的頻率分辨率和低的時(shí)間分辨率而短窗口有高的時(shí)間分辨率和低的頻率分辨率，即選用的窗函數(shù)確定之后，對(duì)應(yīng)的時(shí)頻分辨率就是固定的，而窗函數(shù)的時(shí)間分辨率和頻率分辨率并沒(méi)有隨著信號(hào)頻率的變化而自適應(yīng)地變化。所以，在信號(hào)的大概頻率未知的情況下，選擇合適的窗函數(shù)是非常關(guān)鍵的步驟。

　　圖4是對(duì)信號(hào)進(jìn)行Cohen類時(shí)頻分布分析的結(jié)果，可以看出分析結(jié)果有很好的時(shí)頻聚焦特性。但是，由于Cohen類時(shí)頻分布本身有交叉項(xiàng)的影響，會(huì)降低時(shí)頻分辨率。如圖4所示，在287 ms~318 ms和1 679 ms~1 710 ms時(shí)間段疊加的兩個(gè)常數(shù)頻率就被掩蓋了圖4Cohen類時(shí)頻分布。所以，雖然Cohen類時(shí)頻分布有非常好的時(shí)頻聚焦性，但是由于其固有交叉項(xiàng)的干擾，其應(yīng)用效果受到影響。

　　圖5是對(duì)信號(hào)進(jìn)行連續(xù)小波變換分析的結(jié)果，連續(xù)小波變換能夠隨頻率的變化自適應(yīng)變化其時(shí)頻分辨率，即在低頻部分有很好的頻率分辨率，在高頻部分有很好的時(shí)間分辨率，解決了短時(shí)傅里葉變換在時(shí)域和頻域上不能同時(shí)自圖5連續(xù)小波變換適應(yīng)變化的問(wèn)題。但是，如圖5所示，連續(xù)小波變換時(shí)頻分析結(jié)果的時(shí)頻聚焦性不好。而且，圖6HilbertHuang變換連續(xù)小波變換的小波基一旦確定，在整個(gè)分析過(guò)程中就無(wú)法被替換，所確定小波基的類型直接影響到信號(hào)分析的效果，如何判斷和選用合適的小波基來(lái)分析信號(hào)是一個(gè)難點(diǎn)。

　　圖6是對(duì)信號(hào)進(jìn)行HilbertHuang變換分析的結(jié)果，可以看出HilbertHuang變換不再受傅里葉變換的限制，不需要預(yù)先設(shè)置基函數(shù)，能夠根據(jù)信號(hào)自身的特點(diǎn)自適應(yīng)地選擇頻帶，其分析結(jié)果有良好的時(shí)頻聚焦特性，在一定程度上改善了小波分析結(jié)果的模糊性，適用于分析非線性、非平穩(wěn)的信號(hào)。HilbertHuang變換作為一種較新的信號(hào)分析方法，其基本理論還不是很完善，還需要一些準(zhǔn)確的定義上的證明［6］。但是，HilbertHuang變換良好的自適應(yīng)性和更精確的時(shí)頻分辨率使其在處理非平穩(wěn)信號(hào)時(shí)具有出色的表現(xiàn)，已經(jīng)引起了工程師和科研工作者的廣泛關(guān)注［7］。

3聲卡和LabVIEW結(jié)合實(shí)現(xiàn)EMD時(shí)頻分析

　　常用的聲卡能夠?qū)β曇粜盘?hào)進(jìn)行雙聲道16位的數(shù)據(jù)采集，而且采集到的數(shù)據(jù)是高保真的，最高采樣率能夠達(dá)到176.4 kHz，這樣一個(gè)較高的采樣精度和采樣率，能夠滿足大多數(shù)科研和工程測(cè)量的需求。值得注意的是，聲卡的輸入端電壓不要超過(guò)1 V。

　　LabVIEW提供了許多便于對(duì)聲卡操作的函數(shù)，通過(guò)在LabVIEW后面板調(diào)用底層的聲卡操作函數(shù)來(lái)實(shí)現(xiàn)對(duì)聲音信號(hào)的采集。LabVIEW還提供了MATLAB程序調(diào)用接口MATLAB Script，立足于LabVIEW自身數(shù)據(jù)采集、數(shù)據(jù)保存、波形顯示的優(yōu)勢(shì)，結(jié)合MATLAB強(qiáng)大的數(shù)值計(jì)算能力，就可以快速地開(kāi)發(fā)出一個(gè)聲音信號(hào)采集分析系統(tǒng)。

　　圖7是LabVIEW聲音采集分析程序，包含了配置聲音參數(shù)、寫入聲音文件、聲音文件波形、調(diào)用MATLAB Script進(jìn)行EMD時(shí)頻分析、聲音文件存儲(chǔ)等功能。在實(shí)驗(yàn)室環(huán)境下，對(duì)兩個(gè)周期的規(guī)律性變化的聲音信號(hào)進(jìn)行采集和EMD時(shí)頻分析。

　　圖8是實(shí)驗(yàn)室環(huán)境下基于聲卡和LabVIEW采集到的聲音信號(hào)的時(shí)域波形，圖9是EMD時(shí)頻分析的結(jié)果。從實(shí)驗(yàn)結(jié)果可以看到，在采集到的0~6 s時(shí)間段一直存在一個(gè)20 Hz左右的低頻信號(hào)，采集到的高頻噪聲由濾波器濾除，信號(hào)經(jīng)過(guò)大約0.5 s的時(shí)間，幅值會(huì)有一個(gè)較大的變化，圖8采集的聲音信號(hào)時(shí)域波形會(huì)有一個(gè)強(qiáng)信號(hào)的輸入，頻率大概在100 Hz~450 Hz之間，而且在這個(gè)0~6 s時(shí)間段大約有兩個(gè)周期的信號(hào)在波動(dòng)。

4結(jié)論

　　本文提出了一種基于聲卡和LabVIEW對(duì)聲音信號(hào)進(jìn)行圖9EMD時(shí)頻分析結(jié)果EMD時(shí)頻分析的系統(tǒng)，這種系統(tǒng)非常適合對(duì)頻率在音頻范圍之內(nèi)（20 Hz~20 kHz）的聲音信號(hào)進(jìn)行頻譜分析，分析結(jié)果有良好的時(shí)頻聚焦性，對(duì)信號(hào)有非常好的自適應(yīng)能力。但是，HilbertHuang變換是一種較新的信號(hào)分析方法，其自身理論還在發(fā)展構(gòu)建當(dāng)中，一些相關(guān)的理論完善工作還有待進(jìn)一步開(kāi)展。

參考文獻(xiàn)

　?。?］ DURAK L. Shiftinvariance of shorttime Fourier transform in fractional Fourier domains ［J］. Journal of the Franklin InstituteEngineering and Applied Mathematics, 2009, 346(2): 136146.

　?。?］ Wang Yan, Wu Xi, Li Wenzao, et al. Analysis of microdoppler signatures of vibration targets using EMD and SPWVD ［J］. Neurocomputing, 2016(171): 4856.［3］ IGLEWSKANOWAK I. Continuous wavelet transforms on ndimensional spheres ［J］. Applied and Computational Harmonic Analysis, 2014,39(2): 248276.

　?。?］ HUANG N E, SHEN Z, LONG S R, et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］. Proc. R. Soc. London. A,1998, 454（1971）: 903995.

　?。?］ 楊建國(guó)．小波分析及其工程應(yīng)用［M］ ．北京：機(jī)械工業(yè)出版社，2005.

　　［6］ 苗晟，王威廉，姚紹文．HilbertHuang變換發(fā)展歷程及其應(yīng)用［J］. 電子測(cè)量與儀器學(xué)報(bào)，2014,28(8)：812818．

　?。?］ 胡愛(ài)軍．HilbertHuang變換在旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)分析中的應(yīng)用研究［D］ 北京：華北電力大學(xué)，2008.