目录
- 简介
- EMD算法原理
- python实现EMD案例
本教程为脑机学习者Rose发表 .QQ交流群：903290195
简介
SSVEP信号中含有自发脑电和大量外界干扰信号，属于典型的非线性非平稳信号。传统的滤波方法通常不满足对非线性非平稳分析的条件，1998年黄鄂提出希尔伯特黄变换(HHT)方法，其中包含经验模式分解(EMD)和希尔伯特变换(HT)两部分。EMD可以将原始信号分解成为一系列固有模态函数(IMF) [1]，IMF分量是具有时变频率的震荡函数，能够反映出非平稳信号的局部特征，用它对非线性非平稳的SSVEP信号进行分解比较合适。
EMD算法原理
任何复杂的信号均可视为多个不同的固有模态函数叠加之和，任何模态函数可以是线性的或非线性的，并且任意两个模态之间都是相互独立的。在这个假设基础上，复杂信号 x ( t ) x(t) x(t)的EMD分解步骤如下：
步骤1：
寻找信号 全部极值点，通过三次样条曲线将局部极大值点连成上包络线，将局部极小值点连成下包络线。上、下包络线包含所有的数据点。
步骤2：
由上包络和下包络线的平均值 m 1 ( t ) m_{1}(t) m1(t) ，得出
h 1 ( t ) = x ( t ) − m 1 ( t ) h_1(t)=x(t)-m_{1}(t) h1(t)=x(t)−m1(t)
若 h 1 ( t ) h_{1}(t) h1(t)满足IMF的条件，则可认为 h 1 ( t ) h_{1}(t) h1(t)是 x ( t ) x(t) x(t)的第一个IMF分量。
步骤3：
若 h 1 ( t ) h_{1}(t) h1(t)不符合IMF条件，则将 h 1 ( t ) h_{1}(t) h1(t)作为原始数据，重复步骤1、步骤2，得到上、下包络的均值 m 11 ( t ) m_{11}(t) m11(t)，通过计算 h 11 ( t ) = h 1 ( t ) − m 11 ( t ) h_{11}(t)=h_{1}(t)-m_{11}(t) h11(t)=h1(t)−m11(t)是否适合IMF分量的必备条件，若不满足，重复如上两步 k k k次，直到满足前提下得到 h 1 k ( t ) = h 1 ( k − 1 ) ( t ) − m 1 k ( t ) h_{1k}(t)=h_{1(k-1)}(t)-m_{1k}(t) h1k(t)=h1(k−1)(t)−m1k(t)。第1个IMF表示如下:
c 1 ( t ) = h 1 k ( t ) c_{1}(t)=h_{1k}(t) c1(t)=h1k(t)
步骤4：
将 c 1 ( t ) c_{1}(t) c1(t)从信号 x ( t ) x(t) x(t)中分离得到：
r 1 = x ( t ) − c 1 ( t ) r_{1}=x(t)-c_{1}(t) r1=x(t)−c1(t)
将 r 1 ( t ) r_{1}(t) r1(t)作为原始信号重复上述三个步骤，循环 n n n次，得到第二个IMF分量 c 2 ( t ) c_{2}(t) c2(t)直到第 n n n个IMF分量 ，则会得出：
{ r 2 ( t ) = r 1 ( t ) − c 2 ( t ) ⋯ r n ( t ) = r n − 1 ( t ) − c n ( t ) \begin{cases} r_{2}(t)=r_{1}(t)-c_{2}(t)\\ \cdots\\ r_{n}(t)=r_{n-1}(t)-c_{n}(t)\\ \end{cases} ⎩⎪⎨⎪⎧r2(t)=r1(t)−c2(t)⋯rn(t)=rn−1(t)−cn(t)
步骤5：
当 r n ( t ) r_{n}(t) rn(t)变成单调函数后，剩余的 r n ( t ) r_{n}(t) rn(t)成为残余分量。所有IMF分量和残余分量之和为原始信号 x ( t ) x(t) x(t)：
x ( t ) = ∑ n c i ( t ) + r n ( t ) x(t)=\sum^{n}c_{i}(t)+r_{n}(t) x(t)=∑nci(t)+rn(t)
用EMD进行滤波的基本思想是将原信号进行EMD分解后，只选取与特征信号相关的部分对信号进行重构。如下图中a部分为原始信号，b部分为将原始信号进行EMD分解获得的6个IMF分量和1个残余分量，c部分为将分解获得的6个IMF分量和1个残余分量进行重构后的信号，可以看出SSVEP信号用EMD分解后，基本上包含了原有信号的全部信息。
图片来源于[1]
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()
参考
[1] 基于稳态视觉诱发电位的脑-机接口系统研究
脑机学习者Rose笔记分享,QQ交流群：903290195
