本教程為腦機學習者Rose原創(轉載請聯系作者授權)發表於公眾號：腦機接口社區(微信號：Brain_Computer).QQ交流群：903290195
簡介
SSVEP信號中含有自發腦電和大量外界干擾信號，屬於典型的非線性非平穩信號。傳統的濾波方法通常不滿足對非線性非平穩分析的條件，1998年黃鄂提出希爾伯特黃變換(HHT)方法，其中包含經驗模式分解(EMD)和希爾伯特變換(HT)兩部分。EMD可以將原始信號分解成為一系列固有模態函數(IMF) [1]，IMF分量是具有時變頻率的震盪函數，能夠反映出非平穩信號的局部特征，用它對非線性非平穩的SSVEP信號進行分解比較合適。
EMD算法原理
任何復雜的信號均可視為多個不同的固有模態函數疊加之和，任何模態函數可以是線性的或非線性的，並且任意兩個模態之間都是相互獨立的。在這個假設基礎上，復雜信號$x(t)$的EMD分解步驟如下：
步驟1：
尋找信號 全部極值點，通過三次樣條曲線將局部極大值點連成上包絡線，將局部極小值點連成下包絡線。上、下包絡線包含所有的數據點。
步驟2：
由上包絡和下包絡線的平均值$m_{1}(t)$ ，得出
$$h_1(t)=x(t)-m_{1}(t)$$
若$h_{1}(t)$滿足IMF的條件，則可認為$h_{1}(t)$是$x(t)$的第一個IMF分量。
步驟3：
若$h_{1}(t)$不符合IMF條件，則將$h_{1}(t)$作為原始數據，重復步驟1、步驟2，得到上、下包絡的均值$m_{11}(t)$，通過計算$h_{11}(t)=h_{1}(t)-m_{11}(t)$是否適合IMF分量的必備條件，若不滿足，重復如上兩步$k$次，直到滿足前提下得到$h_{1k}(t)=h_{1(k-1)}(t)-m_{1k}(t)$。第1個IMF表示如下:
$$c_{1}(t)=h_{1k}(t)$$
步驟4：
將$c_{1}(t)$從信號$x(t)$中分離得到：
$$r_{1}=x(t)-c_{1}(t)$$
將$r_{1}(t)$作為原始信號重復上述三個步驟，循環$n$次，得到第二個IMF分量$c_{2}(t)$直到第$n$個IMF分量 ，則會得出：
$$\begin{cases}
r_{2}(t)=r_{1}(t)-c_{2}(t)\
\cdots\
r_{n}(t)=r_{n-1}(t)-c_{n}(t)\
\end{cases}
$$
步驟5：
當$r_{n}(t)$變成單調函數后，剩余的$r_{n}(t)$成為殘余分量。所有IMF分量和殘余分量之和為原始信號$x(t)$：
$$x(t)=\sum^{n}c_{i}(t)+r_{n}(t)$$
用EMD進行濾波的基本思想是將原信號進行EMD分解后，只選取與特征信號相關的部分對信號進行重構。如下圖中a部分為原始信號，b部分為將原始信號進行EMD分解獲得的6個IMF分量和1個殘余分量，c部分為將分解獲得的6個IMF分量和1個殘余分量進行重構后的信號，可以看出SSVEP信號用EMD分解后，基本上包含了原有信號的全部信息。
python實現EMD案例
# 導入工具庫
import numpy as np
from 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和剩余信號res
emd = EMD()
emd.emd(S)
imfs, res = emd.get_imfs_and_residue()
# 繪制 IMF
vis = Visualisation()
vis.plot_imfs(imfs=imfs, residue=res, t=t, include_residue=True)
# 繪制並顯示所有提供的IMF的瞬時頻率
vis.plot_instant_freq(t, imfs=imfs)
vis.show()
腦機學習者Rose筆記分享,QQ交流群：903290195
更多分享，請關注公眾號