点击上面"脑机接口社区"关注咱们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源创计划”，欢迎正在阅读的你也加入，一块儿分享。