康世勛，孔德杰，馮進(jìn)良，馬晨陽(yáng)

（長(zhǎng)春理工大學(xué) 光電工程學(xué)院，長(zhǎng)春 130022）

數(shù)字信號(hào)處理是許多領(lǐng)域理論和應(yīng)用的重要組成部分，一些傳統(tǒng)的信號(hào)處理方法如傅里葉變換是處理線性系統(tǒng)和平穩(wěn)信號(hào)的有效方法，但是現(xiàn)實(shí)世界中的大多數(shù)系統(tǒng)都是非線性，且信號(hào)是非平穩(wěn)的[1]。

經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）算法是希爾伯特-黃變換（Hilbert-Huang Transform，HHT）的主要組成部分。其在結(jié)構(gòu)動(dòng)力分析、醫(yī)學(xué)檢測(cè)、地震檢測(cè)以及設(shè)備故障診斷[2-5]等工程領(lǐng)域有著重要的意義。其篩選過(guò)程需要大量的計(jì)算，利用軟件實(shí)現(xiàn)EMD 的高速計(jì)算是困難的。但是在許多應(yīng)用中，特別是與安全相關(guān)的應(yīng)用，如軸承故障檢測(cè)和轉(zhuǎn)子條斷裂檢測(cè)，高速分析和反應(yīng)是必要的。基于上述問(wèn)題，提出一種提高EMD 計(jì)算速度的解決方案。

為了提高EMD 的計(jì)算速度，現(xiàn)場(chǎng)可編程邏輯門(mén)陣列（FPGA）憑借其低功耗、高能耗比、并行化處理和可重配置的優(yōu)點(diǎn)得到了研究人員的重視，并且其并行化與流水線處理方式能夠很好地處理EMD 算法的高計(jì)算量。例如周浩等人[6]提出了基于FPGA 的EMD 算法的實(shí)現(xiàn)，提高了EMD 算法的計(jì)算速度及其靈活性。

然而，在實(shí)現(xiàn)EMD 算法的實(shí)現(xiàn)步驟中，因?yàn)镃SI 需要求取二階導(dǎo)數(shù)，并且需要執(zhí)行大量的乘除法運(yùn)算，為了進(jìn)一步提高EMD 算法在FPGA 上的采樣率與計(jì)算速度，提出了使用計(jì)算更為簡(jiǎn)單的鋸齒變換（ST）來(lái)代替CSI 生成包絡(luò)。以ST為基礎(chǔ)的EMD 算法可以更好地適應(yīng)硬件環(huán)境并可以提高采樣率與計(jì)算速度，使得在FPGA 上實(shí)現(xiàn)的EMD 算法可以更快的速度處理各個(gè)頻率的信號(hào)。

1 EMD 算法的原理

EMD 算法是由Huang 等人[1]提出的，它可以將任何一種非線性和非平穩(wěn)信號(hào)分解成若干個(gè)本征模函數(shù)（Intrinsic Mode Function，IMF）。

1.1 IMF 的特性

由EMD 算法產(chǎn)生的IMF 需要滿足以下兩個(gè)條件：

（1）在函數(shù)的整個(gè)時(shí)間范圍內(nèi)，局部極值點(diǎn)與函數(shù)零點(diǎn)的個(gè)數(shù)必須相等或最多只差一個(gè)。

（2）在任意時(shí)間點(diǎn)，由局部極大值定義的上包絡(luò)線和局部極小值定義的下包絡(luò)線的均值為0。

1.2 EMD 算法的步驟

對(duì)于任意的信號(hào)x（t），EMD 算法可以分為以下幾個(gè)步驟：

（1）求取輸入信號(hào)x(t)的極大值與極小值；

（2）分別使用CSI 或ST（本文使用的是ST）連接輸入信號(hào)的極大值和極小值，生成上包絡(luò)線u(t)和下包絡(luò)線l(t)；

（3）計(jì)算局部均值，即mean(t) =u(t) +l(t)/2；

（4）計(jì)算hj(t)，hj(t)=x(t)- meant，其 中j為生成IMF 的索引；

（5）如果hj(t) 滿足成為IMF 的條件，則hj(t)是一個(gè)IMF，c1(t)=hj(t)；若hj(t) 不滿足，則令x(t)=hj(t)，重新從步驟（1）開(kāi)始；

（6）計(jì)算第一個(gè)殘差信號(hào)，將殘差信號(hào)記作r1(t)=x(t)-c1(t)，將r1(t) 作為新的輸入信號(hào)，重復(fù)步驟（1）～（4）得到下一個(gè)IMF。當(dāng)最終的IMF為單調(diào)分量即信號(hào)不可再分解時(shí)，無(wú)法得到進(jìn)一步的IMF，EMD 算法的分解結(jié)束，經(jīng)EMD 分解后，最初的信號(hào)x(t)可以寫(xiě)成：

其中，得到的IMF個(gè)數(shù)為n個(gè)；第k個(gè)IMF記為ck(t)；最終殘差記為rn(t)。

整個(gè)EMD 算法的流程圖如圖1 所示。

圖1 EMD 算法流程圖

2 ST的原理及其在EMD算法中的應(yīng)用

由Lu[7]提出了基于鋸齒變換（ST）的EMD 算法更適合硬件環(huán)境的實(shí)現(xiàn)，基于ST 的EMD 算法不同于三次樣條插值（CSI）在原始數(shù)據(jù)上進(jìn)行處理，它是將原始數(shù)據(jù)函數(shù)轉(zhuǎn)換到鋸齒空間，然后在鋸齒空間中，連接函數(shù)的極大值與極小值構(gòu)建上下包絡(luò)線。IMF 為變換后的函數(shù)與上下包絡(luò)均值的差值，最后再將IMF 變換到原始的函數(shù)空間，得到IMF 分量。這種方法對(duì)函數(shù)進(jìn)行一次數(shù)據(jù)處理，得到唯一的IMF 分量，且不需要進(jìn)行EMD 算法的篩選過(guò)程[7]，速度更快，可預(yù)測(cè)性更強(qiáng)，且對(duì)硬件環(huán)境更為友好。

2.1 將原始函數(shù)變換到鋸齒空間

將原始函數(shù)變換到鋸齒空間中需要其對(duì)應(yīng)的鋸齒函數(shù)，鋸齒函數(shù)可以通過(guò)用直線連接函數(shù)的極值點(diǎn)實(shí)現(xiàn)，如圖2 所示為原始函數(shù)及其對(duì)應(yīng)的鋸齒函數(shù)。在鋸齒函數(shù)的每一段上，其極點(diǎn)值是與原函數(shù)的極大極小值一一對(duì)應(yīng)的，而原始函數(shù)數(shù)據(jù)的變化是與鋸齒函數(shù)的直線段一一對(duì)應(yīng)的。

圖2 原始函數(shù)與對(duì)應(yīng)的鋸齒波函數(shù)

鋸齒變換可以分為以下幾個(gè)步驟：

（1）設(shè)原始信號(hào)有m個(gè)極值點(diǎn)，記為：

（2）設(shè)有k個(gè)極大值組成上包絡(luò)線U(t)，記作：

設(shè)有l(wèi)個(gè)極小值組成下包絡(luò)線L(t)，記作：

其中，k+l=m。

（3）將原函數(shù)的鋸齒變換定義為：

除了函數(shù)上的點(diǎn)，原函數(shù)坐標(biāo)空間也需要進(jìn)行變換為對(duì)應(yīng)的鋸齒空間，在原函數(shù)空間中，點(diǎn)的坐標(biāo)為(t,y)；在鋸齒空間中，點(diǎn)的坐標(biāo)為(u,s)，將坐標(biāo)上的鋸齒變換定義為：

從式（5）和式（6）中可以看出鋸齒變換并不對(duì)縱坐標(biāo)（函數(shù)的值）做變換，而是對(duì)橫坐標(biāo)（時(shí)間）進(jìn)行壓縮或者拓展，使函數(shù)變?yōu)榉侄蔚匿忼X函數(shù)。

2.2 在鋸齒空間中求函數(shù)上下包絡(luò)線及均值

鋸齒函數(shù)在極大值U(ti)和極小值L(ti)之間連續(xù)的線性變換。

求上包絡(luò)線可以通過(guò)連接連續(xù)的極大值取出來(lái)，記作：

下包絡(luò)線同樣可以通過(guò)連接連續(xù)的極小值取出來(lái)，記作：

極大值和極小值的和等于極值的個(gè)數(shù)；

殘差是包絡(luò)的平均值，記作：

IMF 是鋸齒函數(shù)與均值的差值，記作：

函數(shù)c(u)具有上下包絡(luò)對(duì)稱(chēng)性的，滿足IMF的兩個(gè)特性，上包絡(luò)線和下包絡(luò)線的均值為0，且在鋸齒空間中IMF 的極大值與極小值是關(guān)于零值對(duì)稱(chēng)的。

圖3 原始函數(shù)的鋸齒變換及在鋸齒空間內(nèi)的上下包絡(luò)線和均值

圖4 鋸齒空間中的IMF 及其上下包絡(luò)線

在鋸齒空間中，IMF 和上下包絡(luò)線都是分段線性的函數(shù)，只需要在極值點(diǎn)處操作，不需要重復(fù)的操作。

2.3 鋸齒空間到原函數(shù)空間的轉(zhuǎn)換

要將在鋸齒空間中得到的IMF、上下包絡(luò)線和殘差轉(zhuǎn)換到原始函數(shù)空間得到cdata(t)、rdata(t)、Udata(t)、Ldata(t)中，只需要將鋸齒空間中的結(jié)果轉(zhuǎn)換為原始空間中的結(jié)果。而鋸齒空間對(duì)函數(shù)的值（縱坐標(biāo)）沒(méi)有進(jìn)行改變，所以只需要在鋸齒空間中找到對(duì)應(yīng)橫坐標(biāo)即可，公式如下：

由于鋸齒反變換并不會(huì)破壞IMF 的包絡(luò)對(duì)稱(chēng)性，所以變換回原空間的函數(shù)依然滿足IMF 的性質(zhì)。

3 基于ST的EMD算法在FPGA上的實(shí)現(xiàn)

基于ST 的EMD 算法將被分解為四個(gè)模塊在FPGA 上實(shí)現(xiàn)，分別是極值尋取模塊、ST 生成包絡(luò)模塊、輸出計(jì)算模塊、門(mén)控時(shí)鐘模塊。

3.1 算法在FPGA 上的整體實(shí)現(xiàn)

本文提出的實(shí)現(xiàn)EMD 算法的方法采用自底向上方法實(shí)現(xiàn)。首先介紹了極值尋取模塊、包絡(luò)模塊和輸出計(jì)算模塊等基本模塊，并使用Verilog-HDL 在Vivado 軟件上實(shí)現(xiàn)，然后將這些模塊整合到一起，形成頂層模塊。

在本文的設(shè)計(jì)中，采樣后的輸入信號(hào)被發(fā)送到極值尋取模塊的輸入緩沖區(qū)中，用于臨時(shí)存儲(chǔ)輸入數(shù)據(jù)，之后用于計(jì)算IMF。如果出現(xiàn)極大值或者極小值的話，那么極值尋取模塊會(huì)把極大值或極小值輸出。此外，該模塊根據(jù)輸出的極大值和極小值情況，分別將極大值或極小值能使信號(hào)置為高電平。包絡(luò)模塊的輸入控制器加載上包絡(luò)線和下包絡(luò)線的輸入數(shù)據(jù)，然后通過(guò)包絡(luò)線生成模塊計(jì)算上下包絡(luò)線并輸出，之后借助輸出計(jì)模塊計(jì)算信號(hào)的IMF 和殘差信號(hào)。在所提出的設(shè)計(jì)中，極值尋取模塊、包絡(luò)線輸入控制器模塊和輸出計(jì)算模塊以比例時(shí)鐘頻率工作，而其余模塊，即包絡(luò)線生成模塊以輸入頻率工作。使用這兩種不同頻率的目的是使系統(tǒng)進(jìn)行實(shí)時(shí)處理。當(dāng)極值尋取模塊忙于掃描輸入時(shí)在求極大值和極小值的同時(shí)，包絡(luò)模塊利用循環(huán)緩沖區(qū)中已經(jīng)存儲(chǔ)的極值來(lái)計(jì)算上包絡(luò)和下包絡(luò)，從而使系統(tǒng)速度更快。整個(gè)算法的結(jié)構(gòu)框圖如圖5 所示。

圖5 算法的整體結(jié)構(gòu)圖

3.2 極值尋取模塊

EMD 算法的第一步是找到輸入信號(hào)的極大值與極小值，在本文提出的設(shè)計(jì)中，有一個(gè)單獨(dú)的模塊用來(lái)實(shí)現(xiàn)這一功能，其結(jié)構(gòu)框圖如圖6 所示。這一模塊主要由一個(gè)輸入緩沖區(qū)（由三個(gè)寄存器組成）、兩個(gè)比較器、兩個(gè)計(jì)數(shù)器、與門(mén)和三態(tài)緩沖器組成。程序開(kāi)始運(yùn)行時(shí)，在每個(gè)時(shí)鐘周期，將輸入信號(hào)發(fā)送到由寄存器X、Y、Z組成的輸入緩沖區(qū)中。在兩個(gè)比較器的作下，寄存器X中的數(shù)據(jù)與寄存器Y中的數(shù)據(jù)作比較，寄存器Z中的數(shù)據(jù)也與寄存器Y中的數(shù)據(jù)作比較。兩個(gè)比較器的輸出G_1 與L_1 通過(guò)一個(gè)與門(mén)（圖6 中A1），兩個(gè)比較器的輸出G_2 與L_2 通過(guò)另一個(gè)與門(mén)（圖6 中A2）。如果在寄存器Y中的數(shù)據(jù)是極大值，寄存器Y中的數(shù)據(jù)比在寄存器X和寄存器Z中的數(shù)據(jù)更大，與門(mén)A1 將會(huì)輸出高電平，輸出極大值的使能端將為高電平。如果寄存器Y中的數(shù)據(jù)小于寄存器X和寄存器Z中的數(shù)據(jù)，那么寄存器Y中的數(shù)據(jù)是極小值，與門(mén)A2將輸出高電平，極小值的使能端將為高電平。極值尋取模塊也有兩個(gè)獨(dú)立的計(jì)數(shù)器。當(dāng)找到極大值時(shí)，極大值計(jì)數(shù)器清除為0，開(kāi)始新的計(jì)數(shù)。當(dāng)找到極小值時(shí)，極小值計(jì)數(shù)器清除為0，開(kāi)始新的計(jì)數(shù)。

圖6 極值尋取模塊

3.3 使用ST 的生成包絡(luò)線模塊

使用ST 的生成包絡(luò)線模塊的輸入是極值尋取模塊所得到的極值。此模塊的實(shí)現(xiàn)借助了有限狀態(tài)機(jī)（FSM），模塊內(nèi)部包括減法、除法、乘法和加法器結(jié)構(gòu)，并且為了進(jìn)一步提高計(jì)算速度，模塊內(nèi)的乘除法器采用了流水線結(jié)構(gòu)。計(jì)算上包絡(luò)線和下包絡(luò)線的結(jié)構(gòu)是一樣的。

計(jì)算包絡(luò)線的模塊的算法原理是基于第二節(jié)提到的ST，如果A0、A1是連續(xù)的極大值或者極小值，C是這兩個(gè)值之間的時(shí)間跨度，那么根據(jù)ST 算法，包絡(luò)線的輸出是：

其中，j是計(jì)數(shù)器從0 到C- 1 的輸出。模塊的輸入與ready 信號(hào)是由包絡(luò)模塊的輸入控制器給出的，下文將詳細(xì)介紹。

當(dāng)ready 信號(hào)的值為高電平時(shí)，包絡(luò)模塊的FSM 控制器通過(guò)減法器獲取A0、A1的差值，然后把差值和C的值輸入到除法器上。FSM 控制器將除法器輸出的結(jié)果與j作為乘法器的輸入（j的初始值設(shè)為0），將乘法器的輸出結(jié)果與A0相加，并輸出結(jié)果。與此同時(shí)FSM 控制器將包絡(luò)線生成模塊的輸出使能信號(hào)置為高電平，表示模塊的輸出結(jié)果在輸出端可用。然后模塊內(nèi)的計(jì)數(shù)器將j的值加1，并將再次輸入到乘法器中，與除法器的輸出結(jié)果相乘，乘法器的結(jié)果再次與A0相加，然后將加法器的輸出結(jié)果放在模塊的輸出處。此后一直重復(fù)這個(gè)過(guò)程，直到j(luò)的值等于C- 1。當(dāng)j的值等于C- 1時(shí)，F(xiàn)SM 將完成信號(hào)置為高電平并將輸出使能信號(hào)置為低電平，表示所有點(diǎn)的計(jì)算已經(jīng)完成。

使用ST 計(jì)算包絡(luò)線模塊的框圖如圖7 所示。

圖7 使用ST 計(jì)算包絡(luò)線模塊

包絡(luò)模塊的輸入控制器（以下簡(jiǎn)稱(chēng)控制器A）的功能是將極值模塊的輸出加載到生成包絡(luò)模塊。因?yàn)槭褂昧藛为?dú)的模塊尋找信號(hào)的極值，因此為了檢測(cè)出是否找到了極大值與極小值，并將其輸入到包絡(luò)模塊，設(shè)計(jì)了控制器A。

控制器A監(jiān)測(cè)極值輸入的使能信號(hào)，如果信號(hào)的極大值被找到，則極大值輸入使能信號(hào)為高電平，極大值計(jì)數(shù)寄存器的值加1，控制器A將信號(hào)的極大值作為A0輸入到上包絡(luò)生成模塊，極大值使能信號(hào)恢復(fù)為低電平；緊接著當(dāng)連續(xù)的下一個(gè)極值被找到時(shí)，極大值輸入使能信號(hào)為高電平，控制器A將該極值作為A1輸入到上包絡(luò)生成模塊，并且極大值計(jì)數(shù)寄存器的值為2，ready 信號(hào)為高電平，表示上包絡(luò)生成模塊已準(zhǔn)備好計(jì)算，并將極值模塊中極大值的計(jì)數(shù)器的輸出作為C輸入到上包絡(luò)生成模塊，極大值輸入使能信號(hào)為低電平，上包絡(luò)輸入模塊開(kāi)始運(yùn)行，極值數(shù)量寄存器的值減1。同樣的，下包絡(luò)線生成模塊的輸入也由相同結(jié)構(gòu)的控制器給出，控制器A的功能如圖8 所示。

圖8 控制器A流程圖

3.4 輸出計(jì)算模塊

輸出計(jì)算模塊的主要任務(wù)是計(jì)算上下包絡(luò)線的平均值，并將均值從輸入信號(hào)中減去得到IMF 信號(hào)。算法模塊由加法器、除法器、減法器和控制器組成?？刂破饔脕?lái)監(jiān)測(cè)包絡(luò)生成模塊中完成信號(hào)的值，如果包絡(luò)生成模塊的完成信號(hào)為高電平，則表明包絡(luò)線計(jì)算完成，將包絡(luò)線的值加載到輸出計(jì)算模塊，計(jì)算模塊的框圖如圖9 所示。

圖9 輸出計(jì)算模塊

3.5 門(mén)控時(shí)鐘模塊

由于在本文的設(shè)計(jì)中，每個(gè)模塊不是同時(shí)工作的，所以使用門(mén)控時(shí)鐘來(lái)降低整個(gè)系統(tǒng)的功耗。在這個(gè)電路中，使能信號(hào)來(lái)自另一個(gè)模塊。由于使能信號(hào)的不同，門(mén)控信號(hào)可能有故障。因此，為了減少故障，門(mén)控時(shí)鐘使用了D鎖存器來(lái)輸入使能信號(hào)。在門(mén)控時(shí)鐘的幫助下可以降低各電路的動(dòng)態(tài)功耗。

在本文的設(shè)計(jì)中，輸入信號(hào)的一些點(diǎn)可能沒(méi)有極大值或極小值可以用來(lái)生成上、下包絡(luò)線。在這種情況下，包絡(luò)生成模塊中的FSM 控制器將等待下一次可以生成包絡(luò)線的點(diǎn)的到來(lái)，模塊中的其他部分的時(shí)鐘信號(hào)將停止翻轉(zhuǎn)，從而降低動(dòng)態(tài)功耗。門(mén)控時(shí)鐘的電路圖如圖10 所示。

圖10 門(mén)控時(shí)鐘電路

將上述各個(gè)模塊及總體設(shè)計(jì)在Vivado 上編譯、綜合、布局布線后，得到的FPGA 資源使用報(bào)告如表1 所示

表1 FPGA 主要資源表

4 仿真結(jié)果與對(duì)比

4.1 仿真的結(jié)果

將第三節(jié)的設(shè)計(jì)編譯、綜合、布線之后再在Vivado 上進(jìn)行波形仿真，仿真的輸入波形是由Matlab 產(chǎn)生的。與DSP 相關(guān)的EMD 算法大部分適用于低頻信號(hào)，本文提出的基于FPGA 的EMD算法可用于高頻信號(hào)的分解，仿真波形是對(duì)頻率為100 kHz、500 kHz、和5 000 kHz 的混合信號(hào)進(jìn)行驗(yàn)證，其仿真結(jié)果在Matlab 如圖11 所示，Vivado 中仿真波形如圖12 所示。經(jīng)仿真后發(fā)現(xiàn)本文提出的EMD 算法可以通過(guò)其不同的模量分解出不同的頻率分量。

圖11 利用Matlab 繪制EMD 輸入信號(hào)及IMF

圖12 Vivado 中IMF 的仿真波形

4.2 兩種設(shè)計(jì)之間的對(duì)比

在本文中，提出了一種不同于傳統(tǒng)的生成包絡(luò)線的方法，本文提出的設(shè)計(jì)使用ST 生成包絡(luò)線。該方法在準(zhǔn)確度方面相對(duì)于傳統(tǒng)的EMD 算法較低（ST 生成的包絡(luò)線與CSI 生成的包絡(luò)線的相關(guān)度為0.969），但是在復(fù)雜度方面大大降低復(fù)雜性的對(duì)比度，如表2 所示，且對(duì)硬件環(huán)境更為友好。在采樣率方面，相對(duì)于在161 kHz 下工作的傳統(tǒng)的EMD 算法，本文提出的設(shè)計(jì)可以在25 MHz 上進(jìn)行采樣，并且計(jì)算速度也有提高，其對(duì)比如表3 所示。

表2 ST 與CSI 的計(jì)算步驟復(fù)雜度對(duì)比

表3 使用ST與使用CSI的EMD算法的采樣率與計(jì)算速度對(duì)比

5 結(jié)論

為了提高EMD 算法的實(shí)時(shí)性能與計(jì)算速度，本文提出了基于ST 的EMD 算法，并在FPGA上進(jìn)行了實(shí)現(xiàn)。相對(duì)于傳統(tǒng)的基于CSI 的EMD算法，本文的設(shè)計(jì)在準(zhǔn)確度方面略微降低，但在采樣率和計(jì)算速度上大幅度提高，提高了算法的實(shí)時(shí)性，使得其在故障檢測(cè)等實(shí)時(shí)性要求較高的場(chǎng)合能夠得到更好的應(yīng)用。并且與其他與DSP 相關(guān)的設(shè)計(jì)大部分用于低頻信號(hào)的分解相比，本文的設(shè)計(jì)也可以用于處理高頻信號(hào)。在以后的工作中，考慮將鋸齒變換得到的結(jié)果進(jìn)行平滑處理，進(jìn)一步優(yōu)化，提高其準(zhǔn)確度。