CN105279376A - 一种基于gpgpu的非线性非稳态复杂信号自适应分解方法 - Google Patents

Abstract

本发明公开了一种基于GPGPU的非线性非稳态复杂信号自适应分解方法，属于信号分析领域。本发明针对传统的EEMD算法在处理大规模信号数据时，受限于算法本身的密集型的计算，导致传统的EEMD算法不能满足实际应用中实时性的需求。本发明通过分析EEMD算法中包含大量的可高度并行计算步骤，通过应用基于CUDA的GPGPU方法对EEMD算法进行并行化设计，使该算法在数据精度和时间消耗上达到一个优化的状态，并结合希尔伯特黄变换，应用希尔伯特变换与香农熵的概念得到希尔伯特-黄谱熵对分解信号进一步研究。实验证明，本方法在实际的信号分解分析中具有更好的效率和可用性。

Description

一种基于GPGPU的非线性非稳态复杂信号自适应分解方法

技术领域

本发明属于信号分析技术领域，涉及一种复杂信号分解方法，尤其涉及一种基于GPGPU的非线性非稳态复杂信号自适应分解方法。

背景技术

人脑是复杂的非线性系统，脑电信号的研究是当今生命科学的重要前沿领域之一。脑电信号处理对于脑部相关疾病的检测、诊断和治疗至关重要，然而脑电信号的研究涉及到大量的脑电信号数据的采集和计算，大量神经数据的存储、管理和利用是一个巨大的挑战。信号分析在实际应用具有大脑是一个高度复杂的非线性、非平稳性系统，脑电信号是大量神经元活动产生的，也具有非线性、非平稳性的特点。传统的线性时频分析方法，如短时傅立叶变换、Winger-Ville分布和小波变换等，已经广泛应用于脑电信号的频谱分析中，但是这些方法的缺点之一就是不能准确估计信号的瞬时振幅和相位/频率(文献[1])。

经验模式分解(EmpiricalModeDecomposition，EMD)作为一种可替代的线性方法，弥补了传统信号分析方法的一些不足，被用于非平稳、非线性的脑电信号分析中。EMD是一种完全数据驱动方法，不同于小波变换需要利用小波函数来分解复杂的信号，而是可以直接分解得到一组包含瞬时特征的固有模态函数(IntrinsicModeFunction，IMF)。EMD已经广泛应用于各种神经信号处理中，例如癫痫信号的检测和预测(文献[2，3])，麻醉状态(文献[4])和睡眠中(文献[5，6])的脑电信号监测。虽然EMD在处理神经信号方面表现优异，但当信号出现间歇性中断时，IMF会包含不同程度的大幅震荡，容易产生模态混叠现象。这也是EMD的最大问题所在，为了克服EMD中出现的混频现象，Huang提出了加噪声进行辅助信号分析(NoiseAssistedDataAnalysis，NADA)的方法——总体平均经验模式分解法(EnsembleEmpiriealModeDecomposition，EEMD)(文献[7])。

EEMD是一种噪声辅助分析方法，算法要求大量的噪声添加，输出的精度取决于进行白噪声添加试验次数是否足够多。因此EEMD的计算是密集型的，大量的、耗时的计算使它受限于实时性的应用，因此亟待发明一种快速有效的信号分解方法。

[文献1]赵治栋,唐向宏,赵知劲,等.基于Hilbert-HuangTransform的心音信号谱分析[J].传感技术学报,2005,18(1):18-22.

[文献2]FineAS,NichollsDP,MogulDJ.AssessingInstantaneousSynchronyofNonlinearNonstationaryOscillatorsintheBrain[J].JournalofNeuroscienceMethods,2010.186(1):42-51.

[文献3]LiXL,JefferysJGR,FoxJ,etal..NeuronalPopulationOscillationsofRatHippocampusduringEpilepticSeizures[J].NeuralNetworks,2008,21(8):1105-1111.

[文献4]PachoriRB,DiscriminationbetweenIctalandSeizure-freeEEGSignalsusingEmpiricalModeDecomposition[J].ResearchLettersinSignalProcessing,2008,2008:1-5.

[文献5]LiXL,LiD,LiangZH,etal.,AnalysisofDepthofAnesthesiawithHilbert-HuangSpectralEntropy[J].ClinicalNeurophysiology,2008,119(11):2465-2475.

[文献6]CausaL,HeldCM,CausaJ,etal.,AutomatedSleep-spindleDetectioninHealthyChildrenPolysomnograms[J].IEEETransactiononBiomedicalEngineering,2010.57(9):2135-2146.

[文献7]WuZH,HuangNE.EnsembleEmpiricalModeDecomposition:ANoiseAssistedDataAnalysisMethod[M].WorldScientific,2008,1(1):1-41.

发明内容

针对现有方法的不足，本发明应用CPU多线程线程池技术和基于CUDA的GPGPU方法对EEMD算法进行并行化设计，使该算法在数据精度和时间消耗上达到一个优化的状态，并结合希尔伯特黄变换(Hilbert-HuangTransform，HHT)，应用希尔伯特变换与香农熵的概念得到希尔伯特-黄谱熵(Hilbert-HuangSpectralEntropy，HHSE)来更提取信号更多的物理特征，HHSE是一个更精确和近连续信号的能量分布，各个IMF和HHSE精确地描述了原始信号的物理本质，应用EEMD和HHT算法对脑电信号进行更深入的研究和分析。

本发明所采用的技术方案是：一种基于GPGPU的非线性非稳态复杂信号自适应分解方法，其特征在于，包括以下步骤：

步骤1：初始化输入信号x(t)，计算本征模态函数IMF(intrinsicmodefunctions，本征模态函数)的数目NI：

NI＝log₂(length)-1

其中，imf_i(i＝1,2,…,NI)表示NI个本征模态函数，length是初始信号x(t)的长度；

步骤2：设置输入信号x(t)的白噪声信号幅度，设置整体实验次数trial为NTE，初始化trial次数k＝1；

步骤3：在第k次trial时，令x_k(t)＝x(t)+n_k(t)，其中n_k(t)表示通过设定幅度添加的白噪声，并且令i＝1；初始化剩余信号为r_i(t)＝x_k(t)；

步骤4：在提取第k次trial的第i个IMF时，令j＝1，h_j-1(t)＝r_i(t)，h_j-1(t)表示第j-1次trial的IMF；

步骤5：计算h_j-1(t)的局部最小值h_min(t)和最大值h_max(t)；

步骤6：对最小值h_min(t)和最大值h_max(t)用三次曲线进行插值，用来提取h_j-1(t)上部包络线和下部包络线，然后计算上部包络线和下部包络线的均值m_j-1(t)，得到第j次的h_j(t)，即：

h_j(t)＝h_j-1(t)-m_j-1(t)；

步骤7：检验h_j(t)是否为IMF；

若否，则j＝j+1，并返回执行步骤5；

若是，第k次trial的第i个IMF值便可以得到，即imf_i[k](t)＝h_j(t)，然后令r_i+1(t)＝r_i(t)-imf_i[k](t)作为新的残余信号，并顺序执行下述步骤8；

步骤8：判断r_i+1(t)是否有2个以上的极值；

若是，则令i＝i+1，并返回执行步骤4；

若否，则分解的过程第k次trial便终止了，因此，这个添加了白噪声的信号x_k(t)能被分解为其中r_NI(t)是x_k(t)的残余信号；并顺序执行下述步骤9；

步骤9：判断k是否大于阈值NTE；

若否，则令k＝k+1，并返回执行步骤3；

若否，则预先设置的全部实验次数结束，并且原始信号x(t)第i个IMF信号便可以是所有次数trial的平均值，即：

${\stackrel{\‾}{imf}}_i\left(t\right)=(\underset{k=1}{\overset{NTE}{\Σ}}{\mathrm{imf}}_i\[k\](t\left)\right)/NTE;$

其中，i＝1，2，…，NTE；

步骤10：在得到IMF分量后，对IMF分量用希尔伯特变换进行处理：

$Z\left(t\right)=imf\left(t\right)+iH\[imf\left(t\right)\]=a\left(t\right)e^{i\∫\mathrm{\ω}\left(t\right)dt};$

其中，IMF分量的瞬时角频率w(t)和幅度a(t)为：

$a\left(t\right)=\sqrt{{\mathrm{imf}}^2\left(t\right)+H^2\[imf\left(t\right)\]}$

$\mathrm{\ω}\left(t\right)=\frac{d}{dt}\[arctan(H\[imf(t)\]/imf\left(t\right))\];$

步骤11：振幅的时频分布被定义为振幅谱，记为H(ω,t)，称为Hilbert谱，用函数(f)来表示(w,t)，即：h(f)＝∫H(w,t)dt；

Hilbert边际谱标准化后得到：

$\widehat{h\left(f\right)}=\frac{h\left(f\right)}{\mathrm{\Σ}h\left(f\right)};$

则熵值表示为：

$H=-\underset{f}{\mathrm{\Σ}}\widehat{h\left(f\right)}\log \left(\widehat{h\left(f\right)}\right);$

通过对熵值的进一步log(N)分解，熵值被归一化为：

$HHSE=\frac{H}{log\left(N\right)};$

式中：HHSE表示希尔伯特黄谱熵，N表示频率的个数。

作为优选，步骤4中对NTE次trial做并行化处理，对于NTE个Trial的处理，主机调用NTE个线程块，每个线程块分配128个线程，每个线程用于计算8个数据点的极值；步骤5中主机调用NTE个线程块，每个线程块分配1个线程用于处理包络线的计算，包络线的计算采用三次样条插值法；主机再次初始化NTE个线程块，每个线程块分配128个线程，调用128×NT个CUDA线程并行完成NTE个Trial的IMF提取运算。

与现有的传统的基于总体经验均值模式方法(EEMD)相比，本发明具有以下优点和有益效果：

(1)本发明通过对传统的EEMD算法进行改进，使得算法能够快速的处理大规模高维数据分析问题；

(2)应用基于CPU多线程的线程池技术和基于GPU多线程的CUDA技术，利用两种技术都达到了EEMD算法并行化的目的，而且并行化后的算法比之前在效率和精度上都有了提高；

(3)同时本发明结合HHT方法对信号进行分析，提取出信号更多的物理特性，为信号研究提供了高效的解决方案。

附图说明

附图1：本发明的流程图。

附图2：本发明实施例的流程图。

附图3：本发明实施例原始信号图。

附图4：本发明实施例中间结果生成图，其中图(a)为癫痫发作间期(0～4s)，图(b)为癫痫发作期(64～68s)。

附图5:本发明实施例信号分解特征提取结果图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本实施例采用的是CUDA4.2(统一计算设备架构),而传统的小波相干性方法是基于C语言实现的。实验的硬件环境是NVIDIAGeForceGTX295显卡，处理器是Intel(R)Core(TM)i7-2600。实验数据取自癫痫病患者的脑电信号，利用6通道(F3，F4，C3，C4，O1和O2)头皮电极进行采集，电极按照10-20系统放置。原始脑电信号序列如图2所示，数据长度80s，采样率为256Hz，在63s处癫痫开始发作。本文主要通过计算HHSE来研究分析癫痫脑电信号。

请见图1和图2，本发明提供的一种基于GPGPU的非线性非稳态复杂信号自适应分解方法，包括以下步骤：

步骤1：初始化输入信号6通道x(t)，计算IMF(intrinsicmodefunctions，本征模态函数)的数目，imf_i(i＝1,2,…,NI)。其计算方式如下：

NI＝log₂(length)-1

length是指初始信号x(t)的长度。

步骤2：在计算时采用滑动窗方法，时间窗大小为1024，重叠750，对于每一个Trial添加噪声的幅度为0.1×(对应时间窗的每一个单元Epoch的标准偏差)，重复添加100次，每个Epoch得到100个Trial，初始化trial次数k＝1；

步骤3：对于包含1024个点的Epoch，初始化之后，主机调用8个线程块，每个线程块分配128个线程来计算1024个点的标准偏差STD，进一步计算每个Epoch与它本身STD的商，即Epoch[i]/STD，对其重复添加100次白噪声得到对应的Trial[j](j＝1,2,…,100)；

步骤4：对于100个Trial的处理，主机调用100个线程块，每个线程块分配128个线程，每个线程用于计算8个数据点的极值；

步骤5：主机调用NT个线程块，每个线程块分配1个线程用于处理上部包络线和下部包络线的计算，上部包络线和下部包络线的计算采用三次样条插值法；

步骤6：主机再次初始化NT个线程块，每个线程块分配128个线程，调用128×NT个CUDA线程并行完成NT个Trial的IMF提取运算。重复循环10次步骤2至4，直到完成所有IMF的提取；

步骤7：在得到IMF分量后，对IMF分量用希尔伯特变换进行处理：

$Z\left(t\right)=imf\left(t\right)+iH\[imf\left(t\right)\]=a\left(t\right)e^{i\∫\mathrm{\ω}\left(t\right)dt}$

其中，IMF分量的瞬时角频率w(t)和幅度a(t)为：

$a\left(t\right)=\sqrt{{\mathrm{imf}}^2\left(t\right)+H^2\[imf\left(t\right)\]}$

$\mathrm{\α}t)=\frac{d}{dt}\[arctan(H\[imf(t)\]/imf(t\left)\right)\];$

步骤8：振幅的时频分布被定义为振幅谱，记为H(ω,t)，称为Hilbert谱，用函数(f)来表示(w,t)，即：h(f)＝∫H(w,t)dt；

Hilbert边际谱标准化后得到：

$\widehat{h\left(f\right)}=\frac{h\left(f\right)}{\mathrm{\Σ}h\left(f\right)};$

则熵值表示为：

$H=-\underset{f}{\mathrm{\Σ}}\widehat{h\left(f\right)}\log \left(\widehat{h\left(f\right)}\right);$

通过对熵值的进一步log(N)分解，熵值被归一化为：

$HHSE=\frac{H}{log\left(N\right)};$

式中：HHSE表示希尔伯特黄谱熵，N表示频率的个数。

经计算后第一通道的希尔伯特黄谱如图3所示，分别选取癫痫发作间期和癫痫发作期两段数据进行分析，图(a)为癫痫发作间期(0～4s)，图(b)为癫痫发作期(64～68s)，相应的希尔伯特黄谱和边际谱如图4，两段数据谱分析的差异可以明显地从图中观察到，在癫痫发作期，可以观察到3-4Hz的活动增加，0-2Hz和30Hz以上的频段都被滤除，中/低频的脑电活动在癫痫发作时增强。图5为6通道数据的HHSE，可以看到熵值在癫痫发作期明显降低。可以说明本发明的对信号特征提取的有效性。

本实例通过比较不同数据大小下，运算速度时间的对比来展示本发明的方法的效率。从表1中本发明可以看出，相比于传统的EEMD方法，本发明的基于GPU加速的方法能够明显提升方法的运行速度。

表1：在不同处理次数不同方法的处理时间(s)

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (2)

1.一种基于GPGPU的非线性非稳态复杂信号自适应分解方法，其特征在于，包括以下步骤：

步骤1：初始化输入信号x(t)，计算本征模态函数IMF的数目NI：

NI＝log₂(length)-1

其中，imf_i(i＝1,2,…,NI)表示表示得到的NI个本征模态函数，length是初始信号x(t)的长度；

步骤2：设置输入信号x(t)的白噪声信号幅度，设置整体实验次数trial为NTE，初始化trial次数k＝1；

步骤3：在第k次trial时，令x_k(t)＝x(t)+n_k(t)，其中n_k(t)表示通过设定幅度添加的白噪声，并且令i＝1；初始化剩余信号为r_i(t)＝x_k(t)；

步骤4：在提取第k次trial的第i个IMF时，令j＝1，h_j-1(t)＝r_i(t)，h_j-1(t)表示第j-1次trial的IMF；

步骤5：计算h_j-1(t)的局部最小值h_min(t)和最大值h_max(t)；

步骤6：对最小值h_min(t)和最大值h_max(t)用三次曲线进行插值，用来提取h_j-1(t)上部包络线和下部包络线，然后计算上部包络线和下部包络线的均值m_j-1(t)，得到第j次的h_j(t)，即：

h_j(t)＝h_j-1(t)-m_j-1(t)；

步骤7：检验h_j(t)是否为IMF；

若否，则j＝j+1，并返回执行步骤5；

若是，第k次trial的第i个IMF值便可以得到，即imf_i[k](t)＝h_j(t)，然后令r_i+1(t)＝r_i(t)-imf_i[k](t)作为新的残余信号，并顺序执行下述步骤8；

步骤8：判断r_i+1(t)是否有2个以上的极值；

若是，则令i＝i+1，并返回执行步骤4；

若否，则分解的过程第k次trial便终止了，因此，这个添加了白噪声的信号x_k(t)能被分解为其中r_NI(t)是x_k(t)的残余信号；并顺序执行下述步骤9；

步骤9：判断k是否大于阈值NTE；

若否，则令k＝k+1，并返回执行步骤3；

若否，则预先设置的全部实验次数结束，并且原始信号x(t)第i个IMF信号便可以是所有次数trial的平均值，即：

${\stackrel{\‾}{imf}}_i\left(t\right)=(\underset{k=1}{\overset{NTE}{\Σ}}{\mathrm{imf}}_i\[k\](t\left)\right)/NTE;$

其中，i＝1，2，…，NTE；

步骤10：在得到IMF分量后，对IMF分量用希尔伯特变换进行处理：

$Z\left(t\right)=imf\left(t\right)+iH\[imf\left(t\right)\]=a\left(t\right)e^{i\∫\mathrm{\ω}\left(t\right)dt};$

其中，IMF分量的瞬时角频率w(t)和幅度a(t)为：

$a\left(t\right)=\sqrt{{\mathrm{imf}}^2\left(t\right)+H^2\[imf\left(t\right)\]}$

$\mathrm{\ω}\left(t\right)=\frac{d}{dt}\[arctan\left(H\[imf\left(t\right)\]/imf\left(t\right)\right)\];$

步骤11：振幅的时频分布被定义为振幅谱，记为H(ω,t)，称为Hilbert谱，用函数(f)来表示(w,t)，即：h(f)＝∫H(w,t)dt；

Hilbert边际谱标准化后得到：

$\widehat{h\left(f\right)}=\frac{h\left(f\right)}{\mathrm{\Σ}h\left(f\right)};$

则熵值表示为：

$H=-\underset{f}{\Σ}\widehat{h\left(f\right)}log\left(\widehat{h\left(f\right)}\right);$

通过对熵值的进一步log(N)分解，熵值被归一化为：

$HHSE=\frac{H}{log\left(N\right)};$

式中：HHSE表示希尔伯特黄谱熵，N表示频率的个数。

2.根据权利要求1所述的基于GPGPU的非线性非稳态复杂信号自适应分解方法，其特征在于：步骤4中对NTE次trial做并行化处理，对于NTE个Trial的处理，主机调用NTE个线程块，每个线程块分配128个线程，每个线程用于计算8个数据点的极值；步骤5中主机调用NTE个线程块，每个线程块分配1个线程用于处理包络线的计算，包络线的计算采用三次样条插值法；主机再次初始化NTE个线程块，每个线程块分配128个线程，调用128×NT个CUDA线程并行完成NTE个Trial的IMF提取运算。