基于ST 的EMD 算法在FPGA 上的設(shè)計與實(shí)現(xiàn)

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

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

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

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

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

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

1 EMD 算法的原理

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

1.1 IMF 的特性

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

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

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

1.2 EMD 算法的步驟

對于任意的信號x（t），EMD 算法可以分為以下幾個步驟：

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

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

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

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

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

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

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

整個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 分量。這種方法對函數(shù)進(jìn)行一次數(shù)據(jù)處理，得到唯一的IMF 分量，且不需要進(jìn)行EMD 算法的篩選過程[7]，速度更快，可預(yù)測性更強(qiáng)，且對硬件環(huán)境更為友好。

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

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

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

鋸齒變換可以分為以下幾個步驟：

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

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

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

其中，k+l=m。

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

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

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

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

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

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

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

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

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

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

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

圖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é)果。而鋸齒空間對函數(shù)的值（縱坐標(biāo)）沒有進(jìn)行改變，所以只需要在鋸齒空間中找到對應(yīng)橫坐標(biāo)即可，公式如下：

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

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

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

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

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

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

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

3.2 極值尋取模塊

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

圖6 極值尋取模塊

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

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

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

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

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

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

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

包絡(luò)模塊的輸入控制器（以下簡稱控制器A）的功能是將極值模塊的輸出加載到生成包絡(luò)模塊。因為使用了單獨(dú)的模塊尋找信號的極值，因此為了檢測出是否找到了極大值與極小值，并將其輸入到包絡(luò)模塊，設(shè)計了控制器A。

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

圖8 控制器A流程圖

3.4 輸出計算模塊

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

圖9 輸出計算模塊

3.5 門控時鐘模塊

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

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

圖10 門控時鐘電路

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

表1 FPGA 主要資源表

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

4.1 仿真的結(jié)果

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

圖11 利用Matlab 繪制EMD 輸入信號及IMF

圖12 Vivado 中IMF 的仿真波形

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

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

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

表3 使用ST與使用CSI的EMD算法的采樣率與計算速度對比

5 結(jié)論

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