Author: Aeo 轉(zhuǎn)載請注明出處

算法背景

經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition，縮寫EMD）是由黃鍔（N. E. Huang）在美國國家宇航局與其他人于1998年創(chuàng)造性地提出的一種新型自適應(yīng)信號時頻處理方法，特別適用于非線性非平穩(wěn)信號的分析處理。

算法過程分析

• 篩選（Sifting）
  1. 求極值點(diǎn) 通過Find Peaks算法獲取信號序列的全部極大值和極小值
  2. 擬合包絡(luò)曲線 通過信號序列的極大值和極小值組，經(jīng)過三次樣條插值法獲得兩條光滑的波峰/波谷擬合曲線，即信號的上包絡(luò)線與下包絡(luò)線
  3. 均值包絡(luò)線 將兩條極值曲線平均獲得平均包絡(luò)線
  4. 中間信號 原始信號減均值包絡(luò)線，得到中間信號
  5. 判斷本征模函數(shù)（IMF） IMF需要符合兩個條件：
     1）在整個數(shù)據(jù)段內(nèi)，極值點(diǎn)的個數(shù)和過零點(diǎn)的個數(shù)必須相等或相差最多不能超過一個。
     2）在任意時刻，由局部極大值點(diǎn)形成的上包絡(luò)線和由局部極小值點(diǎn)形成的下包絡(luò)線的平均值為零，即上、下包絡(luò)線相對于時間軸局部對稱。
• IMF 1 獲得的第一個滿足IMF條件的中間信號即為原始信號的第一個本征模函數(shù)分量IMF 1
  （由原數(shù)據(jù)減去包絡(luò)平均后的新數(shù)據(jù)，若還存在負(fù)的局部極大值和正的局部極小值，說明這還不是一個本征模函數(shù)，需要繼續(xù)進(jìn)行“篩選”。）
• 使用上述方法得到第一個IMF后，用原始信號減IMF1，作為新的原始信號，再通過上述的篩選分析，可以得到IMF2，以此類推，完成EMD分解。

1. 求極值點(diǎn)

from scipy.signal import argrelextrema

# 通過Scipy的argrelextrema函數(shù)獲取信號序列的極值點(diǎn)
data = np.random.random(100)
max_peaks = argrelextrema(data, np.greater)
min_peaks = argrelextrema(data, np.less)

# 繪制極值點(diǎn)圖像
plt.figure(figsize = (18,6))
plt.plot(data)
plt.scatter(max_peaks, data[max_peaks], c='r', label='Max Peaks')
plt.scatter(min_peaks, data[min_peaks], c='b', label='Max Peaks')
plt.legend()
plt.xlabel('time (s)')
plt.ylabel('Amplitude')
plt.title("Find Peaks")

極大值與極小值點(diǎn)

2. 擬合包絡(luò)函數(shù)

這一步是EMD的核心步驟，也是分解出本征模函數(shù)IMFs的前提。

from scipy.signal import argrelextrema

data = np.random.random(100)-0.5
index = list(range(len(data)))

# 獲取極值點(diǎn)
max_peaks = list(argrelextrema(data, np.greater)[0])
min_peaks = list(argrelextrema(data, np.less)[0])

# 將極值點(diǎn)擬合為曲線
ipo3_max = spi.splrep(max_peaks, data[max_peaks],k=3) #樣本點(diǎn)導(dǎo)入，生成參數(shù)
iy3_max = spi.splev(index, ipo3_max) #根據(jù)觀測點(diǎn)和樣條參數(shù)，生成插值

ipo3_min = spi.splrep(min_peaks, data[min_peaks],k=3) #樣本點(diǎn)導(dǎo)入，生成參數(shù)
iy3_min = spi.splev(index, ipo3_min) #根據(jù)觀測點(diǎn)和樣條參數(shù)，生成插值

# 計算平均包絡(luò)線
iy3_mean = (iy3_max+iy3_min)/2

# 繪制圖像
plt.figure(figsize = (18,6))
plt.plot(data, label='Original')
plt.plot(iy3_max, label='Maximun Peaks')
plt.plot(iy3_min, label='Minimun Peaks')
plt.plot(iy3_mean, label='Mean')
plt.legend()
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("Cubic Spline Interpolation")

三次樣條插值法擬合包絡(luò)函數(shù)

用原信號減去平均包絡(luò)線即為所獲得的新信號，若新信號中還存在負(fù)的局部極大值和正的局部極小值，說明這還不是一個本征模函數(shù)，需要繼續(xù)進(jìn)行“篩選”。

第一次篩選后的新信號與原始信號對比

獲取本征模函數(shù)（IMF）

def sifting(data):
    index = list(range(len(data)))

    max_peaks = list(argrelextrema(data, np.greater)[0])
    min_peaks = list(argrelextrema(data, np.less)[0])

    ipo3_max = spi.splrep(max_peaks, data[max_peaks],k=3) #樣本點(diǎn)導(dǎo)入，生成參數(shù)
    iy3_max = spi.splev(index, ipo3_max) #根據(jù)觀測點(diǎn)和樣條參數(shù)，生成插值

    ipo3_min = spi.splrep(min_peaks, data[min_peaks],k=3) #樣本點(diǎn)導(dǎo)入，生成參數(shù)
    iy3_min = spi.splev(index, ipo3_min) #根據(jù)觀測點(diǎn)和樣條參數(shù)，生成插值

    iy3_mean = (iy3_max+iy3_min)/2
    return data-iy3_mean


def hasPeaks(data):
    max_peaks = list(argrelextrema(data, np.greater)[0])
    min_peaks = list(argrelextrema(data, np.less)[0])
    
    if len(max_peaks)>3 and len(min_peaks)>3:
        return True
    else:
        return False


# 判斷IMFs
def isIMFs(data):
    max_peaks = list(argrelextrema(data, np.greater)[0])
    min_peaks = list(argrelextrema(data, np.less)[0])
    
    if min(data[max_peaks]) < 0 or max(data[min_peaks])>0:
        return False
    else:
        return True

    
def getIMFs(data):
    while(not isIMFs(data)):
        data = sifting(data)
    return data


# EMD函數(shù)
def EMD(data):
    IMFs = []
    while hasPeaks(data):
        data_imf = getIMFs(data)
        data = data-data_imf
        IMFs.append(data_imf)
    return IMFs


# 繪制對比圖
data = np.random.random(1000)-0.5
IMFs = EMD(data)
n = len(IMFs)+1

# 原始信號
plt.figure(figsize = (18,15))
plt.subplot(n, 1, 1)
plt.plot(data, label='Origin')
plt.title("Origin ")

# 若干條IMFs曲線
for i in range(0,len(IMFs)):
    plt.subplot(n, 1, i+2)
    plt.plot(IMFs[i])
    plt.ylabel('Amplitude')
    plt.title("IMFs "+str(i+1))

plt.legend()
plt.xlabel('time (s)')
plt.ylabel('Amplitude')

迭代分解結(jié)果

至此就完成了經(jīng)驗?zāi)B(tài)分解EMD的基本過程。

另代碼文件PyEMD.py附在Github上

Reference

【百度百科】經(jīng)驗?zāi)B(tài)分解

【知乎】這篇文章能讓你明白經(jīng)驗?zāi)B(tài)分解（EMD）——基礎(chǔ)理論篇

【知乎】這篇文章能讓你明白經(jīng)驗?zāi)B(tài)分解（EMD）——IMF的物理含義

【CSDN】Python實(shí)現(xiàn)線性插值和三次樣條插值