點擊上面"腦機接口社區"關注咱們node
更多技術乾貨第一時間送達python
網友Aeo[2]提供了下面的算法過程分析。git
算法過程分析github
-
篩選（Sifting） -
求極值點 經過Find Peaks算法獲取信號序列的所有 極大值和極小值 -
擬合包絡曲線 經過信號序列的 極大值和極小值組，通過三次樣條插值法得到兩條光滑的波峯/波谷擬合曲線，即信號的上包絡線與下包絡線 -
均值包絡線 將兩條極值曲線平均得到 平均包絡線 -
中間信號 原始信號減均值包絡線，獲得 中間信號 -
判斷本徵模函數（IMF） IMF須要符合兩個條件：1）在整個數據段內，極值點的個數和過零點的個數必須相等或相差最多不能超過一個。2）在任意時刻，由局部極大值點造成的上包絡線和由局部極小值點造成的下包絡線的平均值爲零，即上、下包絡線相對於時間軸局部對稱。 -
IMF 1 得到的第一個知足IMF條件的中間信號即爲原始信號的第一個本徵模函數份量 IMF 1（由原數據減去包絡平均後的新數據，若還存在負的局部極大值和正的局部極小值，說明這還不是一個本徵模函數，須要繼續進行「篩選」。） -
使用上述方法獲得第一個IMF後，用原始信號減IMF1，做爲新的原始信號，再經過上述的篩選分析，能夠獲得IMF2，以此類推，完成EMD分解。
下面利用公式來講明上面的分析過程。web
EMD算法步驟算法
任何複雜的信號都可視爲多個不一樣的固有模態函數疊加之和，任何模態函數能夠是線性的或非線性的，而且任意兩個模態之間都是相互獨立的。在這個假設 基礎上，複雜信號
的EMD分解步驟以下：
步驟1：
尋找信號 所有極值點，經過三次樣條曲線將局部極大值點連成上包絡線，將局部極小值點連成下包絡線。上、下包絡線包含全部的數據點。安全
步驟2：
由上包絡和下包絡線的平均值
，得出微信
若 知足IMF的條件，則可認爲 是 的第一個IMF份量。markdown
步驟3：
若
不符合IMF條件，則將
做爲原始數據，重複步驟一、步驟2，獲得上、下包絡的均值
，經過計算
是否適合IMF份量的必備條件，若不知足，重複如上兩步
次，直到知足前提下獲得。第1個IMF表示以下:網絡
步驟4：
將
從信號
中分離獲得：
將 做爲原始信號重複上述三個步驟，循環 次，獲得第二個IMF份量 直到第 個IMF份量 ，則會得出：
步驟5：
當
變成單調函數後，剩餘的
成爲殘餘份量。全部IMF份量和殘餘份量之和爲原始信號
：
案例1---Python實現EMD案例
結合上面的算法分析過程，從代碼角度來看看這個算法。
1.求極大值點和極小值點
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")
2. 擬合包絡函數
這一步是EMD的核心步驟，也是分解出本徵模函數IMFs的前提。
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")
用原信號減去平均包絡線即爲所得到的新信號，若新信號中還存在負的局部極大值和正的局部極小值，說明這還不是一個本徵模函數，須要繼續進行「篩選」。
獲取本徵模函數（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) #樣本點導入，生成參數 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')
案例2---利用PyEMD工具來實現EMD
# 導入工具庫import numpy as npfrom PyEMD import EMD, Visualisation
構建信號
時間t: 爲0到1s,採樣頻率爲100Hz，S爲合成信號
# 構建信號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()
參考
[1] 基於穩態視覺誘發電位的腦-機接口系統研究
[2]案例1 代碼來源於 https://github.com/tianyagk
文章來源於網絡，僅用於學術交流，不用於商業行爲
如有侵權及疑問，請後臺留言，管理員即時刪侵！
更多閱讀
腦機接口拼寫器是否真的安全？華中科技大學研究團隊對此作了相關研究
腦機接口BCI學習交流QQ羣：515148456
微信羣請掃碼添加，Rose拉你進羣
(請務必填寫備註，eg. 姓名+單位+專業/領域/行業)
長按關注咱們
歡迎點個在看鼓勵一下
本文分享自微信公衆號 - 腦機接口社區（Brain_Computer）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。