基于聽覺原理的諧波頻率自動提取方法

李允公， 葉利麗， 鄭國威， 毛怒濤， 應遠軍

(1.東北大學機械工程與自動化學院 沈陽，110819)(2.松下家電研究開發(杭州)有限公司 杭州，310018)

基于聽覺原理的諧波頻率自動提取方法

李允公1， 葉利麗1， 鄭國威2， 毛怒濤2， 應遠軍2

(1.東北大學機械工程與自動化學院 沈陽，110819)(2.松下家電研究開發(杭州)有限公司 杭州，310018)

基于聽覺系統的運行機制，提出一種可自動識別和提取信號中諧波頻率的計算方法。首先，對信號進行基底膜帶通濾波；其次，模擬聽皮層時頻感受野現象對不同時頻域內信號進行二維傅里葉變換并提取幅值和頻率信息，基于幅值信息的側抑制結果，提出了一種諧波成分頻率點的判據方法；最后，對提取得到的諧波頻率點進行重組，從而可揭示信號中的主要頻率成分以及各成分的時變情況。數值仿真和試驗數據驗證結果表明，所提方法可準確提取信號中的各諧波頻率，并具有一定的抗噪聲干擾能力，對于微弱諧波成分的頻率提取亦有較好效果。

聽覺模型；頻率提取；帶通濾波；側抑制；故障診斷

引 言

在設備故障診斷中，通常需明確振動信號所包含的主要頻率分量及其時變情況，有時信號的頻率結構比幅值信息更為重要[1-2]。實際中常遇到的一個問題是，當信號結構較為復雜、諧波分量幅值微弱、存在隨機噪聲干擾或缺少必要的輔助信息時(如鍵相位或編碼器信號)，除了頻譜中幅值突出的離散譜線之外，僅通過人為觀察往往難以準確判斷頻譜中哪些譜線對應諧波成分及哪些譜線對應隨機噪聲，雖然結合時頻譜圖可提高判斷的準確性，但仍可能遺漏幅值微弱的頻率成分。同時，對于信號數據量過大或現場連續監測的場合，則需由計算機對每一頻率成分的性質進行自動判斷。因此，有必要研究一種能夠自動識別、提取與跟蹤諧波頻率成分的信號分析方法。

目前，在缺少鍵相位信號或設備參數等輔助信息時，人工的信號分析思路多是搜索頻譜中幅值較高的譜線所對應的頻率值；進而尋找倍頻及邊頻值;若進行自動識別，則多是基于圖像處理的方法從時頻譜圖中提取頻率值等信號特征[3-5]。顯然，當存在噪聲干擾和諧波幅值微弱時，兩類方法均難以獲得良好的效果，且均未充分利用信號的內部信息。

人類聽覺系統可有效區分不同頻率和不同頻率結構的聲音，并能夠跟蹤同一聲源的時變頻率[6]。同時，在從外耳至聽覺中樞這一完整的信息處理通道上，頻率始終是一個的重要參量而被不斷使用[7]。例如，聽覺皮層中不同的神經元具有不同的時頻感受野[8]，即只對特定的頻率成分產生最大的響應，從而可實現頻率的識別和提取。因此，筆者借鑒聽覺系統的運行機制，通過模擬耳蝸基底膜頻率分解、側抑制和時頻感受野等生物學現象，提出了一種可實現諧波頻率識別、提取與跟蹤功能的計算方法，數值仿真和試驗數據檢驗均獲得較好效果。

1 方法的基本原理

所提方法的實現過程如圖1所示。其中，帶通濾波模擬耳蝸基底膜的頻率分解功能，二維傅里葉變換(discrete Fourier transform, 簡稱DFT)用于提取各局部時頻區域的頻率和幅值信息，其中的幅值信息由側抑制計算實現峰值銳化。頻率點提取工作將篩選出諧波成分頻率點，再將同一諧波對應的頻率點歸為一組，并略去與隨機噪聲對應的頻率點，從而完成頻率點重組。所得結果可描述信號的主要頻率結構和各頻率成分的時變情況。

圖1 所提方法的實現過程

Fig.1 Schematic diagram of the proposed method

2 方法的實現

2.1 帶通濾波

設被分析的信號為x(t)。使用聽覺模型中常用的Gammatone基底膜濾波器組[9]進行帶通濾波，設濾波器數目為M，第m個濾波器的表達式為

h(m,t)=B4t3e-2πBtcos(2πfmt+φm)

(1)

其中：fm為第m個濾波器的中心頻率；相位φm通常取為零；參數B用于控制帶寬。

B的計算公式為

B=1.019(24.7+0.108fm)

(2)

為避免同一頻率成分在不同濾波通道中出現相位差異，在基于傅里葉變換進行卷積運算時忽略濾波器h(m,t)的相頻特性，即

y(m,t)=F-1[X(f)|H(m,f)|]

(3)

其中：F-1為傅里葉逆變換；X(f)和H(m,f)分別為x(t)和h(m,t)的頻域表達。

設y(m,t)的離散形式為y(m,n)(n=1,2,…,N)。

2.2 二維DFT

大腦聽皮層中不同的神經元具有不同的時頻感受野，即只對特定的頻率和波形產生最大響應，目前，模擬這一功能的常用方法是對聽覺譜進行二維DFT[10]或二維濾波[11]，但聽覺譜具有頻率模糊性，筆者直接對y(m,n)的各時頻區域做二維DFT計算。

使用二維高斯窗wij(m,n)進行時頻區域選取，窗函數表達式為

(4)

其中：Wr和Wc為wij(m,n)在m和n方向上給定的寬度。

二維DFT為

(5)

其中：2L1+1和2L2+1為兩個方向上的二維DFT點數。

設采樣頻率為fs，高斯窗覆蓋的頻率范圍為Δf，令

Gij(p,q)實現了對y(m,n)的維數擴張，有助于提取信號中的更深層次信息。結合本研究目的，首先搜索|Gij(p,q)|中的最大值，即

A(i,j)=max{|Gij(p,q)|}

(8)

設A(i,j)對應的坐標為(Ωij,Γij)，其中Ωij為wij(m,n)覆蓋的時頻區域內能量最大的諧波頻率。

2.3 側抑制

側抑制[12]是視覺、聽覺及觸覺等感官系統中的一種生物現象，其作用是銳化邊界或波峰，可提高輸入信號的對比度。在本研究中，由于中心頻率相近的Gammatone濾波器之間會具有較大的頻率重合度，所以一頻率成分往往會使A(i,j)中存在覆蓋范圍較大的凸起曲面，為易于實現后續的頻率點提取，首先對A(i,j)進行側抑制處理。

考慮到A(i,j)是離散的，令A(0,j)=0和A(M+1,j)=0,首先計算

AΔ(i,j)=[A(i,j)-A(i-1,j)]

[A(i,j)-A(i+1,j)]

(9)

(i∈{1,2,…,M}，j∈{1,2,…,N})

對A(i,j)的側抑制處理為

(10)

利用A1(i,j)生成頻率函數ζ(i,j)，即

(11)

2.4 頻率點提取

頻率函數ζ(i,j)是時間-濾波通道-頻率空間內的離散點，但有一部分對應干擾噪聲，筆者綜合j時刻下ζ(i,j)中各相鄰非零點的間隔與A1(i,j)的幅值信息完成諧波頻率點的提取。

2.4.1 諧波頻率點判據條件

信號x(t)中的某一頻率分量為

x0(t)=asin(2πf0t)

(12 )

設x(t)中的白噪聲為no(t)，經帶通濾波后，在中心頻率與f0相近的濾波通道中會存在幅值明顯的濾波信號，設第k個濾波器的中心頻率與f0最接近，則該濾波信號為

y(k,t)=|H(k,f0)|asin(2πf0t)+nok(t) (13)其中：nok(t)為no(t)的濾波信號。

nok(t)的幅頻曲線以fk為中心，其包絡與|H(k,f)|相似，考慮到式(8)中只提取了|Gij(p,q)|的最大值信息，所以y(k,t)可近似表示為

y(k,t)≈|H(k,f0)|asin(2πf0t)+γksin(2πfkt)

(14)

其中：γk為nok(t)在fk處的幅值。

同理，與第k個通道相鄰的濾波信號可表示為

y(k+η,t)=|H(k+η,f0)|asin(2πf0t)+

γk+ηsin(2πfk+ηt)

(15)

其中：η∈[-η1,η2]，η1和η2均為正整數，理論上可取值至另一頻率分量出現為止。

同時，x(t)中的其他信號分量的濾波信號可近似設為

y(m,t)≈γmsin(2πfmt)

(16)

由式(5)可知，當i∈[k-η1-Wr/2,k+η2+Wr/2]時，窗函數wij(m,n)的覆蓋范圍內都會包含x0(t)和干擾噪聲的濾波信號，將兩種濾波信號對應的幅值表達為集合的形式，即

(17)

因為γi±l對應的fi±l互不相等，如果

(18)

即

其中：i′為max{Θi}對應的濾波通道數，可得Ωij=f0。

進一步考慮到高斯窗的不斷平移，若一直可滿足式(12)成立，則有

(i∈[k-η1-Wr/2,k+η2+Wr/2])

(19)

所以對于實測信號，若噪聲水平滿足

(20)

或

(21)

則可知ζ(k,j)=f0，但與其鄰近的η3或η4個點處的值則均為零。

另一方面，當窗函數wij(m,n)覆蓋的區域內只有干擾噪聲的濾波信號時，即

A(i,j)=γi

(22)

由于噪聲的隨機性，γi關于i值也會具有隨機性，所以經側抑制后得到的頻率函數會較為密集，很難出現較大的幅值為零的區間。

綜上可得諧波頻率點的判據條件：若ζ(i,j)中一非零點在i方向上與兩側相鄰的非零點間的間隔都大于σ1Wr，即可認為該點對應一諧波分量，則ζ(i,j)保留原值；否則令ζ(i,j)=0。其中，σ1為給定的系數。考慮到實測信號中的干擾噪聲難以避免，本研究令σ1=0.45。

2.4.2 補充條件

當信號中的兩相鄰諧波分量的間隔小于σ1Wr時，上述判據方法會發生漏判。當兩分量的能量均大于噪聲能量時，會在A1(i,j)的相應位置處存在明顯幅值。因此，設j時刻下A1(i,j)的最大值為Λj=max{A1(i,j),i=1,2，…,M}，若A1(i,j)≥σ2Λj,可保留ζ(i,j)原值，否則令ζ(i,j)=0。其中，σ2為給定的系數，本研究令σ2=0.1。

2.4.3 頻率點提取過程示意

圖2為對一混有高斯白噪聲的正弦信號的頻率提取過程。由圖可見，A1(i,j)中只保留了A(i,j)的峰值信息，同時，在正弦信號頻率附近存在寬度約為Wr的空白區間，該區間外的頻率點則較為密集，對應噪聲分量。由這一過程可以發現，二維高斯窗的不斷平移與二維DFT可使諧波成分在時頻域內具有自身的形態特征，當噪聲水平一定時，其形態特征與諧波幅值關系不大，從而可保證頻率點提取的準確性，還有助于提取幅值微弱的諧波頻率點。

圖2 頻率提取過程示意Fig.2 Illustration of frequency extraction

2.5 頻率點重組

頻率函數ζ(i,j)是一組離散點，但其中大多數點的值為零，為便于后續工作，首先提取ζ(i,j)中的非零點。設ζ(iu,j)為ζ(i,j)中第j個時刻下的第u個非零點，u∈{1,2,…,U},U為ζ(i,j)中各時刻下非零點數目中的最大值。生成一函數ξ(u,j)，令ξ(u,j)=ζ(iu,j)。

由于實測振動信號中的諧波分量數目往往較為穩定，且頻率值不會在短時間內產生劇烈波動，加之前期計算中已有效降低了噪聲的影響，所以同一諧波分量的頻率點多會在ξ(u,j)中具有一定的連續性，即聯接成段，且每段中的點對應的u坐標值相同。因此，首先從ξ(u,j)中分離出各頻率點段，并計算ξ(u,j)中相鄰點間的頻率變化量，即

Δξ(u,j)=ξ(u,j+1)-ξ(u,j)

(23)

(24)

本研究取σ2=50 Hz。因同一信號分量的頻率點可能會被劃分成多個數據段，故以不同數據段間首尾點在時頻網格中的距離最近為主要依據進行數據段組合，設兩數據段D(k)和D(l)首尾點距離為

(25)

3 方法驗證

3.1 仿真驗證

設濾波器數目M=200，高斯窗函數參數Wr和Wc分別取為40和1 000，窗函數在時間方向的移動步長為100個點，在頻率方向為1個點，其他參數在前面已給出。

生成仿真信號x(t)如下

x(t)=sin(2π1 000t3)+sin(2π3 000t2)+no(t)

其中：no(t)是方差為0.01的高斯白噪聲。

計算得到的各中間結果如圖3所示，包括基底膜帶通濾波結果y(m,n)、二維DFT后得到的A(i,j)、側抑制結果A1(i,j)以及頻率函數ζ(i,j)。其中，圖3(d)為頻率值投影圖。由圖3可知，使用二維DFT可以有效提取各時頻區域內的主要頻率成分，側抑制處理可使各頻率成分在時頻平面上占據的寬度變窄，從而簡化了時頻信息。由圖3(d)可發現，有用信號與噪聲信號間存在空白區域。

圖3 仿真信號分析所得中間結果Fig.3 Intermediate results in simulation analysis

圖4 函數ξ(u,j)中頻率點分布情況Fig.4 Distribution of frequency points in ξ(u,j)

圖5 頻率點重組結果Fig.5 Result of frequency recombination

在圖4、圖5中的低頻部分存在頻率點丟失現象，這與窗函數的寬度Wr有關。如取Wr=10而其他參數不變，所得結果如圖7所示，可見，頻率點丟失現象已明顯改善。同時在本研究驗證中也發現，過小的Wr會導致過多的噪聲頻率點被漏判，因此，宜使用Wr略大的窗函數。

圖6 頻率重組后各數據段內點數Fig.6 Points count of data segments after frequency recombination

圖7 Wr=10時的頻率點重組結果Fig.7 Result of frequency recombination while Wr=10

3.2 試驗驗證

現結合兩實例驗證所提方法性能：a.某型渦輪泵啟動過程振動信號的頻率提?。籦.吸塵器振動信號的諧波成分分析。

3.2.1 試驗驗證1

渦輪泵是一種應用廣泛的旋轉設備，由于結構緊湊且轉子剛度與支撐剛度均較大，其振動信號中必然存在諧波成分。某型渦輪泵常出現轉子振幅超標現象，使用電渦流位移傳感器測試得到啟動過程轉子位移振動信號，如圖8所示，測試時采樣頻率為51 200 Hz，測點位于近渦輪端。取濾波器數目M=300，Wc=3 000，窗函數在時間軸方向移動步長為150個點。帶通濾波結果y(m,t)如圖9所示，圖中主要有渦輪泵旋轉頻率及其2倍頻成分隨時間的演變過程，但其他成分十分微弱。

圖8 渦輪泵啟動過程振動信號Fig.8 Vibration signal of turbopump startk-up procedure

圖9 帶通濾波結果y(m,t)Fig.9 Band-pass filtering result y(m,t)

使用本研究方法做進一步分析，頻率點重組后得到85個數據段，各數據段內的頻率點數目如圖10所示，其中12個數據段內的點數大于15，各段內的頻率值如圖11所示。由圖11可以發現：a.所提方法可以準確地提取設備的旋轉頻率及其各高次倍頻，且時變特征也與實際情況相符；b.信號中的微弱成分得到了有效表征。除了旋轉頻率的高次倍頻成分外，當設備穩速運行時，在低于旋轉頻率的頻段內存在一斷續出現的頻率成分，雖然這一成分在圖9中亦隱約可見，但無法明確其具體性質，而根據圖11則可判斷為諧波成分。在其出現的時段內，旋轉頻率為1 331Hz，低頻成分的頻率值為512Hz左右，為旋轉頻率的0.38倍，考慮與軸承保持架振動有關。

圖10 頻率重組后各數據段內點數Fig.10 Points count of data segments after frequency recombination

圖11 頻率點重組結果Fig.11 Result of frequency recombination

綜合上述計算結果，推測渦輪泵轉子軸系和軸承兩部分的裝配質量欠佳，改進相關裝配工藝后，振動超標問題得以改善。

3.2.2 試驗驗證2

圖12為某型臥式吸塵器，在對該產品的噪聲控制方法研究過程中進行了機體的振動測試。該設備的核心部件為串勵電機和風機，電機工作轉速約為3.9×104r/min，額定功率為1.5kW，風機葉片數為9。吸塵器的主要激勵源是電機和與之相聯的風機，而電機的轉速單一且穩定，所以，振動信號中必然存在電機的旋轉頻率成分和風機的通過頻率成分。

圖12 吸塵器與振動測點布置Fig.12 Cleaner and vibration measurement points arrangement

在吸塵器工作狀態下進行振動測試，采樣頻率為19 692Hz，測得的部分信號波形如圖13所示。

圖13 吸塵器振動信號Fig.13 Vibration signal of cleaner

使用本研究方法分析所測信號，相關參數與試驗1相同。帶通濾波結果如圖14所示，圖中的低頻部分可見一穩定的諧波分量，但難以明確其他諧波成分。圖15為頻率重組后各數據段內的頻率點數目情況。頻率重組結果如圖16所示，圖中在98，662，1 325和5 947 Hz頻率處存在十分穩定的信號分量。結合電機與風機參數并考慮到二維DFT的頻率誤差可知：98 Hz對應電機的電磁振動；662和1 325 Hz分別為電機的旋轉頻率及其二倍頻；5 947 Hz則對因風機葉片的通過頻率；而在圖14中只能明確98 Hz這一頻率成分。

做為對比，圖17為信號的幅值譜及其局部放大，圖中，除了98Hz處可見明顯的離散譜線，其他3個頻率處的譜線均不突出或被淹沒在由氣流導致的隨機振動頻率成分中，如1 325Hz處的譜線無法確定對應的是周期成分還是隨機振動成分，而5 947 Hz處的譜線幅值又極其微弱，約為98Hz頻率處幅值的0.005倍。由于本研究方法在提取頻率點過程中以頻率點間的間隔為主要判斷依據，所以能夠較好地提取微弱的頻率成分，并具有一定的抗噪聲干擾性能。

圖14 帶通濾波結果y(m,t)Fig.14 Band-pass filtering result y(m,t)

圖15 頻率重組后各數據段內點數Fig.15 Points count of data segments after frequency recombination

圖16 頻率點重組結果Fig.16 Result of frequency recombination

圖17 吸塵器振動信號幅值譜Fig.17 Frequency spectrum of vibration signal of cleaner

4 結 論

1) 基于聽覺系統的時頻感受野功能特性，對信號的不同時頻區域進行二維DFT計算，可提取各區域內的主要諧波頻率，并可精簡信號的時頻信息。

2) 經側抑制計算后，諧波與隨機成分對應的頻率點呈現不同的分布狀態，可根據頻率點間的頻率間隔進行諧波成分頻率點的識別和提取，從而具有一定的抗噪聲干擾能力，并有可能提取到幅值微弱的諧波成分頻率點。

3) 由于所提方法可自動實現對信號頻率結構的描述，對于時變信號又具有頻率跟蹤功能，從而可為設備狀態監測提供輔助信息。

4) 由計算過程可知，對于具有密集頻譜特點的信號，該方法尚無法保證準確提取所有諧波成分頻率點，這也將在后續工作中進行研究。

[1] 楊凱，丁開平，楊云鶴. 基于信號時頻分析理論識別時變模態參數實驗[J].振動、測試與診斷，2015,35(5):880-884.

Yang Kai, Ding Kaiping, Yang Yunhe. Experimental investigation on estimation of time-varying modal parameters using time-frequency analysis[J].Journal of Vibration,Measurement & Diagnosis, 2015,35(5):880-884.(in Chinese)

[2] 肖良瑜,李建偉,宋大鳳,等.立式屏蔽電機半速渦動異常振動試驗分析[J].振動、測試與診斷, 2015,35(2):316-321.

Xiao Liangyu,Li Jianwei,Song Dafeng,et al. Analysis of abnormal vibration of half speed eddy for vertical canned motor[J].Journal of Vibration, Measurement & Diagnosis, 2015,35(2):316-321. (in Chinese)

[3] 甄莉,彭真明. 基于廣義S變換的圖像局部時頻分析[J].航空學報，2008,29(4):1013-1019.

Zhen Li, Peng Zhenming. Local time-frequency analysis of images based on generalized S-transform[J]. Acta Aeronautica et Astronautica Sinica，2008, 29(4):1013-1019. (in Chinese)

[4] 黃朝明，于洪亮，關德林，等. 摩擦振動時頻圖像特征提取[J].振動與沖擊，2012,31(7):46-49,62.

Huang Chaoming, Yu Hongliang, Guang Delin, et al. Feature extraction of frictional vibration based on time-frequency image[J]. Journal of Vibration and Shock, 2012,31(7):46-49,62. (in Chinese)

[5] 權貝貝. WinCE平臺下信號采集與處理系統的開發[D].沈陽:東北大學,2014.

[6] Du Yi, Kong Lingzhi, Wang Qian, et al.Auditory frequency-following response: a neurophysiological mea-sure for studying the “cocktail-party problem”[J]. Neuroscience and Biobehavioral Reviews, 2011, 35 (10):2046-2057.

[7] 李洋，唐佳，付子英，等. 聽覺中樞神經元對聲信號的識別和處理[J]. 生物化學與生物物理進展，2011,38(6):499-505.

Li Yang, Tang Jia, Fu Ziying, et al. Sound signal recognition and processing of central auditory neurons[J]. Progress in Biochemistry and Biophysics，2011,38(6):499-505.(in Chinese)

[8] 劉冬云，趙巖，肖中舉. 短聲改變聽覺中樞下丘神經元對純音的反應特性[J].華南國防醫學雜志，2013,27(8):536-539.

Liu Dongyun, Zhao Yan, Xiao Zhongju. Click change the response property of inferior colliculers neurons to tone of mouse[J]. Military Medical Journal of South China，2013,27(8):536-539. (in Chinese)

[9] 李允公，張金萍，戴麗. 基于極值點概率密度和聽覺模型的瞬態信號提取方法研究[J].振動與沖擊，2015,34(21):37-44,53.

Li Yungong, Zhang Jinping, Dai Li. A method for extracting transient signal based on the probability density of extreme points and auditory model[J]. Journal of Vibration and Shock, 2015,34(21):37-44,53. (in Chinese)

[10] Ezzat T, Bouvrie J, Poggio T. Max-Gabor analysis and synthesis of spectrograms[C]∥Ninthinternational Conference on Spoken Language Processing(ICASLP 2006). New Your: Gurrant Associates, 2006:17-21.

[11] Coensel B D, Botteldooren D. A model of saliency-based auditory attention to environmental sounds[C]∥Proceedings of 20th International Congress on Acoustics, (ICA 2010).New York: Gurrant Associates,2010:23-27.

[12] 徐新洲，羅昕煒，方世良，等. 基于聽覺感知機理的水下目標識別研究進展[J].聲學技術,2013, 32(2): 150-158．

Xu Xinzhou, Luo Xinwei,Fang Shiliang,et al. Research progress of underwater target recognition based on auditory perception mechanism[J]. Technical Acoustics, 2013, 32(2):150-158.(in Chinese)

10.16450/j.cnki.issn.1004-6801.2017.06.008

國家自然科學基金資助項目(51275080)

2016-03-02；

2016-06-04

TH17

李允公，男，1976年6月生，博士、副教授。主要研究方向為機械故障診斷、工程信號分析等。曾發表《Auditory-model-based feature extraction method for mechanical faults diagnosis》(《Chinese Journal of Mechanical Engineering》2010，Vol.21,No.3 )等論文。

E-mail: ygli@mail.neu.edu.cn