手把手教你EMD算法原理与Python实现(更新)

脑机接口社区

发布于 2020-07-15 10:37:13

5.9K0

发布于 2020-07-15 10:37:13

Rose今天主要介绍一下EMD算法原理与Python实现。关于EMD算法之前介绍过《EMD算法之Hilbert-Huang Transform原理详解和案例分析》,

SSVEP信号中含有自发脑电和大量外界干扰信号，属于典型的非线性非平稳信号。传统的滤波方法通常不满足对非线性非平稳分析的条件，1998年黄鄂提出希尔伯特黄变换(HHT)方法，其中包含经验模式分解(EMD)和希尔伯特变换(HT)两部分。EMD可以将原始信号分解成为一系列固有模态函数(IMF) [1]，IMF分量是具有时变频率的震荡函数，能够反映出非平稳信号的局部特征，用它对非线性非平稳的SSVEP信号进行分解比较合适。

网友Aeo[2]提供了下面的算法过程分析。

算法过程分析

筛选（Sifting）
1. 求极值点 通过Find Peaks算法获取信号序列的全部极大值和极小值
2. 拟合包络曲线 通过信号序列的极大值和极小值组，经过三次样条插值法获得两条光滑的波峰/波谷拟合曲线，即信号的上包络线与下包络线
3. 均值包络线 将两条极值曲线平均获得平均包络线
4. 中间信号 原始信号减均值包络线，得到中间信号
5. 判断本征模函数（IMF） IMF需要符合两个条件：1）在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个。2）在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零，即上、下包络线相对于时间轴局部对称。
IMF 1 获得的第一个满足IMF条件的中间信号即为原始信号的第一个本征模函数分量IMF 1（由原数据减去包络平均后的新数据，若还存在负的局部极大值和正的局部极小值，说明这还不是一个本征模函数，需要继续进行“筛选”。）
使用上述方法得到第一个IMF后，用原始信号减IMF1，作为新的原始信号，再通过上述的筛选分析，可以得到IMF2，以此类推，完成EMD分解。

下面是EMD算法思路.

下面利用公式来说明上面的分析过程。

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分量的必备条件，若不满足，重复如上两步

次，直到满足前提下得到

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)

作为原始信号重复上述三个步骤，循环

次，得到第二个IMF分量

c_{2}(t)

直到第

个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分解后，基本上包含了原有信号的全部信息。

图片来源于[1]

案例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.5
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

# 绘制图像
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


# 判断IMFs
def 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.5
IMFs = 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 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和剩余
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] 基于稳态视觉诱发电位的脑-机接口系统研究

[2]案例1 代码来源于 https://github.com/tianyagk

文章来源于网络，不用于商业行为，若有侵权及疑问，请后台留言！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-07-12，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自脑机接口社区微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度