CN111580654A - 一种基于emd的脑电信号的短时特征提取方法 - Google Patents

一种基于emd的脑电信号的短时特征提取方法 Download PDF

Info

Publication number
CN111580654A
CN111580654A CN202010378897.7A CN202010378897A CN111580654A CN 111580654 A CN111580654 A CN 111580654A CN 202010378897 A CN202010378897 A CN 202010378897A CN 111580654 A CN111580654 A CN 111580654A
Authority
CN
China
Prior art keywords
signal
time
maximum
time point
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202010378897.7A
Other languages
English (en)
Inventor
罗元
罗飞
张毅
刘鹏飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing University of Post and Telecommunications
Original Assignee
Chongqing University of Post and Telecommunications
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing University of Post and Telecommunications filed Critical Chongqing University of Post and Telecommunications
Priority to CN202010378897.7A priority Critical patent/CN111580654A/zh
Publication of CN111580654A publication Critical patent/CN111580654A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F3/00—Input arrangements for transferring data to be processed into a form capable of being handled by the computer; Output arrangements for transferring data from processing unit to output unit, e.g. interface arrangements
    • G06F3/01—Input arrangements or combined input and output arrangements for interaction between user and computer
    • G06F3/011—Arrangements for interaction with the human body, e.g. for user immersion in virtual reality
    • G06F3/015—Input arrangements based on nervous system activity detection, e.g. brain waves [EEG] detection, electromyograms [EMG] detection, electrodermal response detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Dermatology (AREA)
  • General Health & Medical Sciences (AREA)
  • Neurology (AREA)
  • Neurosurgery (AREA)
  • Human Computer Interaction (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明请求保护一种基于EMD的脑电信号的短时特征提取方法，该方法包括步骤：S1，通过窗函数对脑电信号进行分段截取；S2，对截取信号的边界极值进行延拓，以缓解EMD在分解短时脑电信号时产生的端点效应；S3，通过EMD将延拓后的信号自适应地分解成多个固有模态函数；S4，对分解的IMF进行Hilbert变换，获得表征时域特征的瞬时能量和表征频域特征的边际能量。本发明在处理短时间序列的脑电信号时，能有效地提取脑电信号特征，保证一定识别率的同时又有效地降低系统的处理时间和延迟时间。

Description

一种基于EMD的脑电信号的短时特征提取方法
技术领域
本发明属于脑机接口中脑电信号处理领域，特别是一种基于EMD的脑电信号的短时特征提取方法。
背景技术
脑-机接口(Brain-Computer Interface，BCI)通过解析输入的脑电信号，将用户的意图解码为控制指令来控制输出设备，实现人脑与外部设备的交互。脑-机接口技术的核心是脑电信号的识别，但脑电信号具有非线性和非平稳性等特点，如何有效地提取脑电信号特征成为识别脑电信号的关键。
当前，脑电信号特征提取主要采用时频域特征分析方法。基于时频域分析的脑电特征提取方法主要有STFT、WT和WPT，然而以上三种算法的本质都是傅里叶变换，都会受到测不准原理的影响，无法同时在时域和频域获得较高的分辨率。经验模态分解(EmpiricalMode Decomposition,EMD)算法能将一段原始的脑电信号自适应地分解成一系列的固有模态函数(IMF)，分解出来的各IMF都包含有原始脑电信号在不同时间尺度上的局部特征，并且能将非平稳的数据平稳化，是一种自适应的时频域分析方法，被广泛地应用于非线性和非平稳信号的分析中。
在EMD分解过程中，每一个IMF的筛选都需要根据信号的局部极值点构造出信号数据的上、下包络线。然而，信号左右两端的端点不可能同时是局部极值点。因此，数据序列两端的上、下包络线会发散，并且随着分解层数的增加，这种发散会逐渐向内扩散，导致出现的误差越来越严重，这种现象被称为端点效应。并且信号数据段越短，端点效应越严重，限制了EMD在短时信号分析中的应用，使得EMD无法适用于要求快速响应的系统中。
发明内容
本发明旨在解决现有技术中存在的端点效应问题，提出一种基于EMD的脑电信号短时特征提取方法。
鉴于此，本发明采用的技术方案，包括以下步骤：
S1，使用窗函数对脑电信号进行分段截取；
S2，对S1中截取信号的边界极值进行延拓，以缓解EMD在分解短时脑电信号时产生的端点效应；
S3，通过EMD将延拓后的信号自适应地分解成多个固有模态函数(IMF)；
S4，对分解的固有模态函数进行希尔伯特变换，获得表征时域特征的瞬时能量和表征频域特征的边际能量。
在上述方案中，步骤S2所述边界极值进行延拓包括以信号端点的一个特征波为依据分别向左和向右延拓两个极大值和两个极小值。
其中，所述向左延拓两个极大值和极小值具体包括：
信号左端第一个特征波包含的信号点数k1：
Figure BDA0002481294660000021
Pm(1)为第1个极大值在序列中对应的位置，Pm(2)为第2个极大值在序列中对应的位置，Pn(1)为第1个极小值在序列中对应的位置，Pn(2)为第2个极小值在序列中对应的位置，分段截取后的信号具有M个极大值和N个极小值；
向外延拓的两个极大值和极小值的时间位置(Tm,Tn)和对应时间的函数值(G,Y)为：
Tm(0)＝Tm(1)-k1Δt,G(0)＝G(1)
Tn(0)＝Tn(1)-k1Δt,Y(0)＝Y(1)
Tm(-1)＝Tm(1)-2k1Δt,G(-1)＝G(1)
Tn(-1)＝Tn(1)-2k1Δt,Y(-1)＝Y(1)
Tm(0)表示向左延拓的第1个极大值对应的时刻点，Tm(1)表示第1个极大值对应的时刻点，G(0)表示向左延拓的第1个极大值对应的函数值，G(1)表示左边第1个极大值对应的函数值，Tn(0)表示向左延拓的第1个极小值对应的时刻点，Tn(1)表示第1个极小值对应的时刻点，Y(0)表示向左延拓的第1个极小值对应的函数值，Y(1)表示第1个极小值对应的函数值，Tm(-1)表示向左延拓的第2个极大值对应的时刻点，G(-1)表示向左延拓的第2个极大值对应的函数值，Tn(-1)表示向左延拓的第2个极小值对应的时刻点，Y(-1)表示向左延拓的第2个极小值对应的函数值，Δt表示采样步长；
其中，所述向右延拓两个极大值和极小值具体包括：
信号右端第一个特征波包含的信号点数k2：
Figure BDA0002481294660000022
Pm(M)表示第M个极大值在序列中对应的位置，Pm(M-1)表示第M-1个极大值在序列中对应的位置，Pn(N)表示第N个极小值在序列中对应的位置，Pn(N-1)表示第N-1个极小值在序列中对应的位置；
向外延拓的两个极大值和极小值的位置为(Tm,Tn)和对应时间的函数值(G,Y)为：
Tm(M+1)＝Tm(M)+k2Δt,G(M+1)＝G(M)
Tn(N+1)＝Tn(N)-k2Δt,Y(N+1)＝Y(N)
Tm(M+2)＝Tm(M)+2k2Δt,G(M+2)＝G(M)
Tn(N+2)＝Tn(N)-2k2Δt,Y(N+2)＝Y(N)
Tm(M+1)表示向右延拓的第1个极大值对应的时刻点，G(M+1)表示向右延拓的第1个极大值对应的函数值，G(M)表示第M个极大值对应的函数值，Tn(N+1)表示向右延拓的第一个极小值对应的时刻点，Tn(N)表示第N个极小值对应的时刻点，Y(N+1)表示向右延拓的的1个极小值点对应的函数值，Y(N)表示第N个极小值对应的函数值，Tm(M+2)表示向右延拓的第2个极大值对应的时刻点，G(M+2)表示向右延拓的第2个极大值对应的函数值，Tn(N+2)表示向右延拓的第2个极小值对应的时刻点，Y(N+2)表示向右延拓的第2个极小值对应的函数值。
进一步地，当端点的数值小于相邻的极小值或大于相邻的极大值时，还包括进行以下处理：
Tm(0)＝t1,G(0)＝x1,x1＞G(1)
Tn(0)＝t1,Y(0)＝x1,x1＜Y(1)
Tm(M+1)＝tn,G(M+1)＝xn,xn＞G(M)
Tn(N+1)＝tn,Y(N+1)＝xn,xn＜Y(N)
t1表示分割的离散信号的第1个时间点，x1表示分割的离散信号的第1个时间点对应的函数值，tn表示分割的离散信号的最后一个时间点，xn表示分割的离散信号的最后一个时间点对应的函数值。
更进一步地，所述EMD分解的具体过程如下：
(1)根据延拓后的信号x(t)，采用三次样条插值拟合出曲线的上下包络线：
g(t)＝f(Tm,Gm,t),y(t)＝f(Tn,Yn,t)
Tm表示极大值对应的时间点，Gm表示极大值对应时间点的函数值，Tn表示极小值对应的时间点，Yn表示极小值对应时间点的函数值，t表示时间点。
(2)计算出上包络线g(t)和下包络线y(t)的平均值：
Figure BDA0002481294660000031
(3)计算出延拓后的信号x(t)和m(t)的差值：
c(t)＝x(t)-m(t)
若c(t)满足IMF的定义截止条件，则c(t)即为分离出来的第一个IMF；否则，x(t)被c(t)取代，重复步骤(1)到(3)，直到满足IMF的定义；
分离出第一个IMF后，计算出剩余信号r(t)：
r(t)＝x(t)-c(t)
将剩余信号r(t)作为原始信号，并重复步骤(1)到(3)。
具体地，所述希尔伯特变换包括对前三阶固有模态函数分量进行希尔伯特变换，以此构造出解析信号，得出希尔伯特谱。包括以下步骤：
对前三阶固有模态函数分量进行希尔伯特变换，具体如下：
Figure BDA0002481294660000041
Yi(t)表示希尔伯特变换结果，π为圆周率，ci(τ)表示第i个IMF分量，t表示时间点，τ表示积分值。
由此构造出解析信号为Zi(t)：
Figure BDA0002481294660000042
其中，Ai(t)为瞬时幅值，
Figure BDA0002481294660000043
为瞬时相位；
根据Ai(t)和
Figure BDA0002481294660000044
进一步求取瞬时频率：
Figure BDA0002481294660000045
则计算出描述信号幅度随时间和瞬时频率的变化规律，即希尔伯特谱：
Figure BDA0002481294660000046
W表示IMF的个数。
根据希尔伯特谱进一步求取希尔伯特瞬时能量谱IES和边际能量谱MES：
Figure BDA0002481294660000047
Figure BDA0002481294660000048
式中，[ω1,ω2]和[t1,t2]分别表示信号的频率范围和时间范围。
本发明的优点及有益效果如下：
针对EMD方法在提取短时间序列的脑电信号特征时，容易受到端点效应的影响，导致分解的结果失真的问题，本发明提供了一种基于EMD的脑电信号短时特征提取方法。本发明的方法相较于其它时频分析方法，不受测不准原理的影响，根据信号自身的特征时间尺度判别内蕴振荡模式，通过三次样条插值将信号自适应地分解成多个具有物理意义的固有模态函数IMF，从而能根据信号自身的特征将信号自适应地分解，在提取短时间序列的脑电信号时。通过将脑电信号截取成短时间序列的脑电信号，并对截取信号的边界极值进行延拓，一方面缓解了EMD在分解短时间序列的脑电信号时产生的端点效应，减小了分解误差，从而使得本方法在分解短时间序列的脑电信号时也能获得较高的识别率；另一方面减少了单次计算的数据量，加快了分解速度，降低了信号的处理时间；同时，对短时间序列的脑电信号进行处理，信号产生到指令执行的延迟时间也随之降低。与传统EMD方法相比，本方法能更好地应用于短时间序列的脑电信号的处理，可以保证识别率的同时降低处理时间的延迟时间，保证一定识别率的同时又能有效地降低系统的处理时间和延迟时间。将该方法集成应用于脑机接口控制系统中，能有效地降低脑电指令的决策时间，提高系统的实时性。
附图说明
图1为本发明提供的一种基于EMD的脑电信号短时特征提取方法的流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、详细地描述。所描述的实施例仅仅是本发明的一部分实施例。
参见图1，本发明具体包括如下步骤：
S1，通过窗函数对原始脑电信号X(t)进行分段截取；计算公式为：
xi(t)＝X(t)wN(t)
式中，X(t)为原始脑电信号，wN(t)为矩形窗函数，xi(t)为分段截取后的信号。
S2，对截取信号xi(t)的边界极值进行延拓，以缓解EMD在分解短时脑电信号时产生的端点效应。具体计算过程如下：
对于截取信号xi(t)，其采样步长Δt＝1/fs，fs表示采样频率，那么xi(t)就可以表示为xi(t)＝[x(t1),x(t2),...,x(tn)]＝[x1,x2,...xn。x(tn)简写为xn，tn，xn。tn表示分割的离散信号的最后一个时间点，xn表示分割的离散信号的最后一个时间点对应的函数值。假设它有M个极大值和N个极小值，那么极值点在序列中对应的位置(Pm,Pn)，该位置对应的时刻点(Tm(i),Tn(j))以及对应的函数值(G(i),Y(j))可分别表示为：
Pm＝[Pm(1),Pm(2),...,Pm(M)]
Pn＝[Pn(1),Pn(2),...,Pn(N)]
Figure BDA0002481294660000051
Figure BDA0002481294660000052
Pm表示极大值在序列中对应的位置，Pn表示极小值在序列中对应的位置，Tm(i)表示极大值Pm对应的时刻点，Tn(j)表示极小值Pn对应的时刻点，G(i)表示极大值Pm对应的函数值，Y(j)表示极小值Pn对应的函数值，Pm(M)表示第M个极大值在序列中对应的位置，Pn(N)表示第N个极小值在序列号中对应的位置，
Figure BDA0002481294660000053
表示极大值Pm对应的时刻点，
Figure BDA0002481294660000054
表示极小值Pn对应的时刻点，
Figure BDA0002481294660000055
表示极大值Pm对应的函数值，
Figure BDA0002481294660000056
表示极小值点Pn对应的函数值。
极值延拓只需要对边界的极值进行延拓，根据极大值和极小值的分布特性以及样条插值的要求，本发明以信号端点的一个特征波为依据分别向左右两端延拓两个极大值和两个极小值。
(1)左端点处理
信号左端第一个特征波包含的信号点数k1：
Figure BDA0002481294660000057
Pm(1)为第1个极大值在序列中对应的位置，Pm(2)为第2个极大值在序列中对应的位置，Pn(1)为第1个极小值在序列中对应的位置，Pn(2)为第2个极小值在序列中对应的位置，分段截取后的信号具有M个极大值和N个极小值；
向外延拓的两个极值的时间位置(Tm,Tn)和函数值(G,Y)为：
Tm(0)＝Tm(1)-k1Δt,G(0)＝G(1)
Tn(0)＝Tn(1)-k1Δt,Y(0)＝Y(1)
Tm(-1)＝Tm(1)-2k1Δt,G(-1)＝G(1)
Tn(-1)＝Tn(1)-2k1Δt,Y(-1)＝Y(1)
Tm(0)表示向左延拓的第1个极大值对应的时刻点，Tm(1)表示第1个极大值对应的时刻点，G(0)表示向左延拓的第1个极大值对应的函数值，G(1)表示左边第1个极大值对应的函数值，Tn(0)表示向左延拓的第1个极小值对应的时刻点，Tn(1)表示第1个极小值对应的时刻点，Y(0)表示向左延拓的第1个极小值对应的函数值，Y(1)表示第1个极小值对应的函数值，Tm(-1)表示向左延拓的第2个极大值对应的时刻点，G(-1)表示向左延拓的第2个极大值对应的函数值，Tn(-1)表示向左延拓的第2个极小值对应的时刻点，Y(-1)表示向左延拓的第2个极小值对应的函数值，Δt表示采样步长；
(2)右端点处理
信号右端第一个特征波包含的信号点数k2：
Figure BDA0002481294660000061
Pm(M)，Pm(M-1)，Pn(N)，Pn(N-1)分别表示什么？(Pm(M)表示第M个极大值在序列中对应的位置，Pm(M-1)表示第M-1个极大值在序列中对应的位置，Pn(N)表示第N个极小值在序列中对应的位置，Pn(N-1)表示第N-1个极小值在序列中对应的位置)
向外延拓的两个极值的位置为(Tm,Tn)和函数值(G,Y)为：
Tm(M+1)＝Tm(M)+k2Δt,G(M+1)＝G(M)
Tn(N+1)＝Tn(N)-k2Δt,Y(N+1)＝Y(N)
Tm(M+2)＝Tm(M)+2k2Δt,G(M+2)＝G(M)
Tn(N+2)＝Tn(N)-2k2Δt,Y(N+2)＝Y(N)
Tm(M+1)表示向右延拓的第1个极大值对应的时刻点，G(M+1)表示向右延拓的第1个极大值对应的函数值，G(M)表示第M个极大值对应的函数值，Tn(N+1)表示向右延拓的第一个极小值对应的时刻点，Tn(N)表示第N个极小值对应的时刻点，Y(N+1)表示向右延拓的的1个极小值点对应的函数值，Y(N)表示第N个极小值对应的函数值，Tm(M+2)表示向右延拓的第2个极大值对应的时刻点，G(M+2)表示向右延拓的第2个极大值对应的函数值，Tn(N+2)表示向右延拓的第2个极小值对应的时刻点，Y(N+2)表示向右延拓的第2个极小值对应的函数值。
当端点的数值小于相邻的极小值或大于相邻的极大值时，需要对其进行额外的处理：
Tm(0)＝t1,G(0)＝x1,x1＞G(1)
Tn(0)＝t1,Y(0)＝x1,x1＜Y(1)
Tm(M+1)＝tn,G(M+1)＝xn,xn＞G(M)
Tn(N+1)＝tn,Y(N+1)＝xn,xn＜Y(N)
t1表示分割的离散信号的第1个时间点，x1表示分割的离散信号的第1个时间点对应的函数值，tn表示分割的离散信号的最后一个时间点，xn表示分割的离散信号的最后一个时间点对应的函数值。
根据延拓的极值点，由三次样条插值求得曲线的上下包络线：
g(t)＝f(Tm,Gm,t),y(t)＝f(Tn,Yn,t)
S3，通过EMD将延拓后的信号将信号自适应地分解成多个固有模态函数(IMF)；具体计算过程如下：
(1)对于延拓后的信号x(t)，用三次样条曲线拟合出极大值组成的上包络线g(t)和极小值组成的下包络线y(t)。
(2)计算出g(t)和y(t)的平均值：
Figure BDA0002481294660000071
(3)计算出x(t)和m(t)的差值：
c(t)＝x(t)-m(t)
若c(t)满足IMF的定义截止条件，则c(t)即为分离出来的第一个IMF；否则，x(t)被c(t)取代，重复计算步骤(1)到步骤(3)，直到满足IMF的定义。
分离出第一个IMF后，计算出剩余信号r(t)：
r(t)＝x(t)-c(t)
将剩余信号r(t)作为原始信号，并重复步骤(1)到步骤(3)。EMD过程结束后，原始信号会被分解成n个IMF分量和一个剩余分量rn(t)，则xi(t)可表示为：
Figure BDA0002481294660000072
S4，对分解的IMF进行Hilbert变换，获得表征时域特征的瞬时能量和表征频域特征的边际能量。具体计算过程如下：
Figure BDA0002481294660000073
Yi(t)表示希尔伯特变换结果，π为圆周率，ci(τ)表示第i个IMF分量，t表示时间点，τ表示积分值。
由此便可以构造出解析信号：
Figure BDA0002481294660000081
其中，Ai(t)为瞬时幅值，
Figure BDA0002481294660000082
为瞬时相位。
根据Ai(t)和
Figure BDA0002481294660000083
可以进一步求取瞬时频率：
Figure BDA0002481294660000084
则可计算出描述信号幅度随时间和瞬时频率的变化规律，即Hilbert谱：
Figure BDA0002481294660000085
W表示IMF个数，H(ω,t)表示Hilbert谱。
根据Hilbert谱可进一步求取Hilbert瞬时能量谱IES和边际能量谱MES：
Figure BDA0002481294660000086
Figure BDA0002481294660000087
式中，[ω1,ω2]和[t1,t2]分别表示信号的频率范围和时间范围。
IES和MES反映了信号在时域和频域上的能量特征，将IES定义为运动想象脑电信号的时域特征，将MES定义为频域特征。
获得以上的IES和MES后，可将其作为脑机接口控制系统中的脑电信号特征，进一步使用支持向量机对其进行分类，识别出脑电信号对应的控制指令，最终控制系统完成相应的运动指令。
以上这些实施例应理解为仅用于说明本发明而不用于限制本发明的保护范围。在阅读了本发明的记载的内容之后，技术人员可以对本发明作各种改动或修改，这些等效变化和修饰同样落入本发明权利要求所限定的范围。

Claims (7)

1.一种基于EMD的脑电信号的短时特征提取方法，其特征在于，包括以下步骤：
S1，使用窗函数对脑电信号进行分段截取；
S2，对S1中截取信号的边界极值进行延拓；
S3，通过EMD将延拓后的信号自适应地分解成多个固有模态函数；
S4，对分解的固有模态函数进行希尔伯特变换，获得表征时域特征的瞬时能量和表征频域特征的边际能量。
2.根据权利要求1所述一种基于EMD的脑电信号的短时特征提取方法，其特征在于：步骤S2所述边界极值进行延拓包括以信号端点的一个特征波为依据分别向左和向右延拓两个极大值和两个极小值。
3.根据权利要求2所述一种基于EMD的脑电信号的短时特征提取方法，其特征在于：所述向左延拓两个极大值和两个极小值具体包括：
信号左端第一个特征波包含的信号点数k1：
Figure FDA0002481294650000011
Pm(1)为第1个极大值在序列中对应的位置，Pm(2)为第2个极大值在序列中对应的位置，Pn(1)为第1个极小值在序列中对应的位置，Pn(2)为第2个极小值在序列中对应的位置，分段截取后的信号具有M个极大值和N个极小值；
向外延拓的两个极大值和极小值的时间位置(Tm,Tn)和对应时间的函数值(G,Y)为：
Tm(0)＝Tm(1)-k1Δt,G(0)＝G(1)
Tn(0)＝Tn(1)-k1Δt,Y(0)＝Y(1)
Tm(-1)＝Tm(1)-2k1Δt,G(-1)＝G(1)
Tn(-1)＝Tn(1)-2k1Δt,Y(-1)＝Y(1)
Tm(0)表示向左延拓的第1个极大值对应的时刻点，Tm(1)表示第1个极大值对应的时刻点，G(0)表示向左延拓的第1个极大值对应的函数值，G(1)表示第1个极大值对应的函数值，Tn(0)表示向左延拓的第1个极小值对应的时刻点，Tn(1)表示第1个极小值对应的时刻点，Y(0)表示向左延拓的第1个极小值对应的函数值，Y(1)第1个极小值对应的函数值，Tm(-1)表示向左延拓的第2个极大值对应的时刻点，G(-1)表示向左延拓的第2个极大值对应的函数值，Tn(-1)表示向左延拓的第2个极小值对应的时刻点，Y(-1)表示向左延拓的第2个极小值对应的函数值，Δt表示采样步长；
所述向右延拓两个极大值和极小值具体包括：
信号右端第一个特征波包含的信号点数k2：
Figure FDA0002481294650000021
Pm(M)表示第M个极大值在序列中对应的位置，Pm(M-1)表示第M-1个极大值在序列中对应的位置，Pn(N)表示第N个极小值在序列中对应的位置，Pn(N-1)表示第N-1个极小值在序列中对应的位置；
向外延拓的两个极大值和极小值的时间位置为(Tm,Tn)和对应时间的函数值(G,Y)为：
Tm(M+1)＝Tm(M)+k2Δt,G(M+1)＝G(M)
Tn(N+1)＝Tn(N)-k2Δt,Y(N+1)＝Y(N)
Tm(M+2)＝Tm(M)+2k2Δt,G(M+2)＝G(M)
Tn(N+2)＝Tn(N)-2k2Δt,Y(N+2)＝Y(N)
Tm(M+1)表示向右延拓的第1个极大值对应的时刻点，G(M+1)表示向右延拓的第1个极大值对应的函数值，G(M)表示第M个极大值对应的函数值，Tn(N+1)表示向右延拓的第1个极小值对应的时刻点，Tn(N)表示第N个极小值对应的时刻点，Y(N+1)表示向右延拓的的1个极小值点对应的函数值，Y(N)表示第N个极小值对应的函数值，Tm(M+2)表示向右延拓的第2个极大值对应的时刻点，G(M+2)表示向右延拓的第2个极大值对应的函数值，Tn(N+2)表示向右延拓的第2个极小值对应的时刻点，Y(N+2)表示向右延拓的第2个极小值对应的函数值。
4.根据权利要求3所述一种基于EMD的脑电信号的短时特征提取方法，其特征在于：当端点的数值小于相邻的极小值或大于相邻的极大值时，还包括进行以下处理：
Tm(0)＝t1,G(0)＝x1,x1＞G(1)
Tn(0)＝t1,Y(0)＝x1,x1＜Y(1)
Tm(M+1)＝tn,G(M+1)＝xn,xn＞G(M)
Tn(N+1)＝tn,Y(N+1)＝xn,xn＜Y(N)
t1表示分割的离散信号的第1个时间点，x1表示分割的离散信号的第1个时间点对应的函数值，tn表示分割的离散信号的最后一个时间点，xn表示分割的离散信号的最后一个时间点对应的函数值。
5.根据权利要求3或4所述一种基于EMD的脑电信号的短时特征提取方法，其特征在于：所述EMD分解的具体过程如下：
(1)根据延拓后的信号x(t)，采用三次样条插值拟合出曲线的上下包络线：
g(t)＝f(Tm,Gm,t),y(t)＝f(Tn,Yn,t)
Tm表示极大值对应的时间点，Gm表示极大值对应时间点的函数值，Tn表示极小值对应的时间点，Yn表示极小值对应时间点的函数值，t表示时间点；
(2)计算出上包络线g(t)和下包络线y(t)的平均值：
Figure FDA0002481294650000031
(3)计算出延拓后的信号x(t)和m(t)的差值：
c(t)＝x(t)-m(t)
若c(t)满足IMF的定义截止条件，则c(t)即为分离出来的第一个IMF；否则，x(t)被c(t)取代，重复步骤(1)到(3)，直到满足IMF的定义；
分离出第一个IMF后，计算出剩余信号r(t)：
r(t)＝x(t)-c(t)
将剩余信号r(t)作为原始信号，并重复步骤(1)到(3)。
6.根据权利要求5所述一种基于EMD的脑电信号的短时特征提取方法，其特征在于：所述希尔伯特变换包括对前三阶固有模态函数分量进行希尔伯特变换，以此构造出解析信号，得出希尔伯特谱。
7.根据权利要求6所述一种基于EMD的脑电信号的短时特征提取方法，其特征在于：所述希尔伯特变换的具体步骤为：
对前三阶固有模态函数分量进行希尔伯特变换，具体如下：
Figure FDA0002481294650000032
Yi(t)表示希尔伯特变换结果，π为圆周率，ci(τ)表示第i个IMF分量，t表示时间点，τ表示积分值；
由此构造出解析信号为Zi(t)：
Figure FDA0002481294650000033
其中，Ai(t)为瞬时幅值，
Figure FDA0002481294650000034
为瞬时相位；
根据Ai(t)和
Figure FDA0002481294650000035
进一步求取瞬时频率：
Figure FDA0002481294650000036
则计算出描述信号幅度随时间和瞬时频率的变化规律，即希尔伯特谱：
Figure FDA0002481294650000037
W表示IMF的个数；
根据希尔伯特谱进一步求取希尔伯特瞬时能量谱IES和边际能量谱MES：
Figure FDA0002481294650000038
Figure FDA0002481294650000039
式中，[ω1,ω2]和[t1,t2]分别表示信号的频率范围和时间范围。
CN202010378897.7A 2020-05-07 2020-05-07 一种基于emd的脑电信号的短时特征提取方法 Pending CN111580654A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010378897.7A CN111580654A (zh) 2020-05-07 2020-05-07 一种基于emd的脑电信号的短时特征提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010378897.7A CN111580654A (zh) 2020-05-07 2020-05-07 一种基于emd的脑电信号的短时特征提取方法

Publications (1)

Publication Number Publication Date
CN111580654A true CN111580654A (zh) 2020-08-25

Family

ID=72114509

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010378897.7A Pending CN111580654A (zh) 2020-05-07 2020-05-07 一种基于emd的脑电信号的短时特征提取方法

Country Status (1)

Country Link
CN (1) CN111580654A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112614505A (zh) * 2020-11-27 2021-04-06 江苏爱谛科技研究院有限公司 一种并行超快速emd信号处理系统
CN113191317A (zh) * 2021-05-21 2021-07-30 江西理工大学 一种基于极点构造低通滤波器的信号包络提取方法和装置
CN113849067A (zh) * 2021-09-26 2021-12-28 华东理工大学 基于经验模态分解的运动想象人工数据生成方法及装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103941091A (zh) * 2014-04-25 2014-07-23 福州大学 基于改进emd端点效应的电力系统hht谐波检测方法
CN105760347A (zh) * 2016-02-04 2016-07-13 福建工程学院 一种基于数据/极值联合对称延拓的hht端点效应抑制方法
CN106371002A (zh) * 2015-07-24 2017-02-01 国网四川省电力公司眉山供电公司 一种基于希尔伯特黄变换算法对断路器故障诊断的方法
CN110514921A (zh) * 2019-07-22 2019-11-29 华南理工大学 一种电力电子变换器非平稳信号中非线性现象的识别方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103941091A (zh) * 2014-04-25 2014-07-23 福州大学 基于改进emd端点效应的电力系统hht谐波检测方法
CN106371002A (zh) * 2015-07-24 2017-02-01 国网四川省电力公司眉山供电公司 一种基于希尔伯特黄变换算法对断路器故障诊断的方法
CN105760347A (zh) * 2016-02-04 2016-07-13 福建工程学院 一种基于数据/极值联合对称延拓的hht端点效应抑制方法
CN110514921A (zh) * 2019-07-22 2019-11-29 华南理工大学 一种电力电子变换器非平稳信号中非线性现象的识别方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SNEHAL R. KHARVATKAR等: "Detection of Pitch Frequency of Indian Classical Music Based on Hilbert-Huang Transform for Automatic Note Transcription", 《2018 FOURTH INTERNATIONAL CONFERENCE ON COMPUTING COMMUNICATION CONTROL AND AUTOMATION (ICCUBEA)》 *
罗飞等: "多特征融合的运动想象脑电特征提取方法", 《计算机应用》 *
黄大吉等: "希尔伯特-黄变换的端点延拓", 《海洋学报》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112614505A (zh) * 2020-11-27 2021-04-06 江苏爱谛科技研究院有限公司 一种并行超快速emd信号处理系统
CN113191317A (zh) * 2021-05-21 2021-07-30 江西理工大学 一种基于极点构造低通滤波器的信号包络提取方法和装置
CN113849067A (zh) * 2021-09-26 2021-12-28 华东理工大学 基于经验模态分解的运动想象人工数据生成方法及装置

Similar Documents

Publication Publication Date Title
CN109307862B (zh) 一种目标辐射源个体识别方法
CN111580654A (zh) 一种基于emd的脑电信号的短时特征提取方法
Lin et al. Elimination of end effects in empirical mode decomposition by mirror image coupled with support vector regression
CN108335694B (zh) 远场环境噪声处理方法、装置、设备和存储介质
CN109785854A (zh) 一种经验模态分解和小波阈值去噪相结合的语音增强方法
CN108009122B (zh) 一种改进的hht方法
CN112560699A (zh) 基于密度和压缩感知的齿轮振动信源欠定盲源分离方法
He et al. Super-gaussian bss using fast-ica with chebyshev–pade approximant
CN116706876A (zh) 一种双高电力系统宽频振荡识别方法、装置及设备
CN117347043A (zh) 一种旋转机械瞬时转速提取方法、装置、设备及存储介质
Xie et al. Adaptive symbolic transfer entropy and its applications in modeling for complex industrial systems
Das et al. Review of adaptive decomposition-based data preprocessing for renewable generation rich power system applications
JPH0573093A (ja) 信号特徴点の抽出方法
Ye et al. Fingerprint image enhancement algorithm based on two dimension emd and gabor filter
CN115267548A (zh) 锂电池电压采样方法、系统和可读存储介质
CN112580701B (zh) 一种基于分类变换扰动机制的均值估计方法及装置
CN113948088A (zh) 基于波形模拟的语音识别方法及装置
CN105513587B (zh) 一种mfcc提取方法及装置
CN107632963A (zh) 长度为复合数的基于fft的电力系统采样信号希尔伯特变换方法
Sumarno On The Performace of Segment Averaging of Discrete Cosine Transform Coefficients on Musical Instruments Tone Recognition
Juang et al. A non-feedback-loop and low-computation-complexity algorithm design for a novel 2-D sliding DFT computation
CN113361819B (zh) 一种线性预测方法及装置
Gensler et al. Fast feature extraction for time series analysis using least-squares approximations with orthogonal basis functions
US20130191430A1 (en) Method for Fast Wavelet Transform of a Signal
CN111354372A (zh) 一种基于前后端联合训练的音频场景分类方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20200825

RJ01 Rejection of invention patent application after publication