法律状态公告日
法律状态信息
法律状态
2018-06-12
授权
授权
2016-02-24
实质审查的生效 IPC(主分类):G06F19/00 申请日:20151021
实质审查的生效
2016-01-27
公开
公开
技术领域
本发明属于信号分析技术领域,涉及一种复杂信号分解方法,尤其涉及一种基于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=log2(length)-1
其中,imfi(i=1,2,…,NI)表示NI个本征模态函数,length是初始信号x(t)的长度;
步骤2:设置输入信号x(t)的白噪声信号幅度,设置整体实验次数trial为NTE,初始化trial次数k=1;
步骤3:在第k次trial时,令xk(t)=x(t)+nk(t),其中nk(t)表示通过设定幅度添加的白噪声,并且令i=1;初始化剩余信号为ri(t)=xk(t);
步骤4:在提取第k次trial的第i个IMF时,令j=1,hj-1(t)=ri(t),hj-1(t)表示第j-1次trial的IMF;
步骤5:计算hj-1(t)的局部最小值hmin(t)和最大值hmax(t);
步骤6:对最小值hmin(t)和最大值hmax(t)用三次曲线进行插值,用来提取hj-1(t)上部包络线和下部包络线,然后计算上部包络线和下部包络线的均值mj-1(t),得到第j次的hj(t),即:
hj(t)=hj-1(t)-mj-1(t);
步骤7:检验hj(t)是否为IMF;
若否,则j=j+1,并返回执行步骤5;
若是,第k次trial的第i个IMF值便可以得到,即imfi[k](t)=hj(t),然后令ri+1(t)=ri(t)-imfi[k](t)作为新的残余信号,并顺序执行下述步骤8;
步骤8:判断ri+1(t)是否有2个以上的极值;
若是,则令i=i+1,并返回执行步骤4;
若否,则分解的过程第k次trial便终止了,因此,这个添加了白噪声的信号xk(t)能被分解为其中rNI(t)是xk(t)的残余信号;并顺序执行下述步骤9;
步骤9:判断k是否大于阈值NTE;
若否,则令k=k+1,并返回执行步骤3;
若否,则预先设置的全部实验次数结束,并且原始信号x(t)第i个IMF信号便可以是所有次数trial的平均值,即:
其中,i=1,2,…,NTE;
步骤10:在得到IMF分量后,对IMF分量用希尔伯特变换进行处理:
其中,IMF分量的瞬时角频率w(t)和幅度a(t)为:
步骤11:振幅的时频分布被定义为振幅谱,记为H(ω,t),称为Hilbert谱,用函数(f)来表示(w,t),即:h(f)=∫H(w,t)dt;
Hilbert边际谱标准化后得到:
则熵值表示为:
通过对熵值的进一步log(N)分解,熵值被归一化为:
式中: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,本征模态函数)的数目,imfi(i=1,2,…,NI)。其计算方式如下:
NI=log2(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分量用希尔伯特变换进行处理:
其中,IMF分量的瞬时角频率w(t)和幅度a(t)为:
步骤8:振幅的时频分布被定义为振幅谱,记为H(ω,t),称为Hilbert谱,用函数(f)来表示(w,t),即:h(f)=∫H(w,t)dt;
Hilbert边际谱标准化后得到:
则熵值表示为:
通过对熵值的进一步log(N)分解,熵值被归一化为:
式中:HHSE表示希尔伯特黄谱熵,N表示频率的个数。
经计算后第一通道的希尔伯特黄谱如图3所示,分别选取癫痫发作间期和癫痫发作期两段数据进行分析,图(a)为癫痫发作间期(0~4s),图(b)为癫痫发作期(64~68s),相应的希尔伯特黄谱和边际谱如图4,两段数据谱分析的差异可以明显地从图中观察到,在癫痫发作期,可以观察到3-4Hz的活动增加,0-2Hz和30Hz以上的频段都被滤除,中/低频的脑电活动在癫痫发作时增强。图5为6通道数据的HHSE,可以看到熵值在癫痫发作期明显降低。可以说明本发明的对信号特征提取的有效性。
本实例通过比较不同数据大小下,运算速度时间的对比来展示本发明的方法的效率。从表1中本发明可以看出,相比于传统的EEMD方法,本发明的基于GPU加速的方法能够明显提升方法的运行速度。
表1:在不同处理次数不同方法的处理时间(s)
应当理解的是,本说明书未详细阐述的部分均属于现有技术。
应当理解的是,上述针对较佳实施例的描述较为详细,并不能因此而认为是对本发明专利保护范围的限制,本领域的普通技术人员在本发明的启示下,在不脱离本发明权利要求所保护的范围情况下,还可以做出替换或变形,均落入本发明的保护范围之内,本发明的请求保护范围应以所附权利要求为准。
机译: 一种基于接收信号噪声比和反馈通道信息的多干扰中继去除系统,该信号具有零精度和最小均方误差技术中的复杂性,该方法使用相同的方法和一种方法录制中号
机译: 基于非线性能量收集的Swipt系统自适应时隙信号接收方法
机译: 基于时变非线性变换的视频信号自适应对比度增强方法及装置