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

下面利用公式來說明上面的分析過程。

EMD算法步驟

任何復雜的信號均可視為多個不同的固有模態函數疊加之和，任何模態函數可以是線性的或非線性的，并且任意兩個模態之間都是相互獨立的。在這個假設 基礎上，復雜信號的EMD分解步驟如下：
步驟1：
尋找信號 全部極值點，通過三次樣條曲線將局部極大值點連成上包絡線，將局部極小值點連成下包絡線。上、下包絡線包含所有的數據點。

步驟2：
由上包絡和下包絡線的平均值 ，得出

若滿足IMF的條件，則可認為是的第一個IMF分量。

步驟3：
若不符合IMF條件，則將作為原始數據，重復步驟1、步驟2，得到上、下包絡的均值，通過計算是否適合IMF分量的必備條件，若不滿足，重復如上兩步次，直到滿足前提下得到。第1個IMF表示如下:

步驟4：
將從信號中分離得到：

將作為原始信號重復上述三個步驟，循環次，得到第二個IMF分量直到第個IMF分量 ，則會得出：

步驟5：
當變成單調函數后，剩余的成為殘余分量。所有IMF分量和殘余分量之和為原始信號：

from scipy.signal import argrelextrema
"""通過Scipy的argrelextrema函數獲取信號序列的極值點"""# 構建100個隨機數data = np.random.random(100)# 獲取極大值max_peaks = argrelextrema(data, np.greater)#獲取極小值min_peaks = argrelextrema(data, np.less)
# 繪制極值點圖像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")

from scipy.signal import argrelextrema
#進行樣條差值import scipy.interpolate as spi
data = np.random.random(100)-0.5index = 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) #樣本點導入，生成參數iy3_max = spi.splev(index, ipo3_max) #根據觀測點和樣條參數，生成插值
ipo3_min = spi.splrep(min_peaks, data[min_peaks],k=3) #樣本點導入，生成參數iy3_min = spi.splev(index, ipo3_min) #根據觀測點和樣條參數，生成插值
# 計算平均包絡線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")

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) #樣本點導入，生成參數    iy3_max = spi.splev(index, ipo3_max) #根據觀測點和樣條參數，生成插值
    ipo3_min = spi.splrep(min_peaks, data[min_peaks],k=3) #樣本點導入，生成參數    iy3_min = spi.splev(index, ipo3_min) #根據觀測點和樣條參數，生成插值
    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

# 判斷IMFsdef 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函數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.5IMFs = 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')

# 導入工具庫import numpy as npfrom PyEMD import EMD, Visualisation

# 構建信號t = np.arange(0,1, 0.01)S = 2*np.sin(2*np.pi*15*t) +4*np.sin(2*np.pi*10*t)*np.sin(2*np.pi*t*0.1)+np.sin(2*np.pi*5*t)

# 提取imfs和剩余emd = EMD()emd.emd(S)imfs, res = emd.get_imfs_and_residue()# 繪制 IMFvis = Visualisation()vis.plot_imfs(imfs=imfs, residue=res, t=t, include_residue=True)# 繪制并顯示所有提供的IMF的瞬時頻率vis.plot_instant_freq(t, imfs=imfs)vis.show()