CN115211869A - 一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统 - Google Patents
一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统 Download PDFInfo
- Publication number
- CN115211869A CN115211869A CN202211015305.0A CN202211015305A CN115211869A CN 115211869 A CN115211869 A CN 115211869A CN 202211015305 A CN202211015305 A CN 202211015305A CN 115211869 A CN115211869 A CN 115211869A
- Authority
- CN
- China
- Prior art keywords
- signal
- imf
- kalman filtering
- component
- noise
- 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
Links
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Pathology (AREA)
- Animal Behavior & Ethology (AREA)
- Psychiatry (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- Veterinary Medicine (AREA)
- Biophysics (AREA)
- Public Health (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- General Health & Medical Sciences (AREA)
- Physiology (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统。针对头部电极采集的信号,发挥了经验模态分解处理非线性非平稳信号的自适应能力,对脑电信号进行自适应分解得到本征模态分量,保证了局部信息的完整性;根据相关性准则对若干本征模态分量进行分类处理;使用卡尔曼滤波算法对本征模态分量信号中与原始脑电信号相似度次高的分量进行去噪处理,避免了卡尔曼滤波在非线性系统下发散问题的同时提升了计算速度,最终得到了去除噪声的脑电信号。本发明是基于数据的信号去噪方法,数据处理步骤完备,选取的算法兼顾了实际脑电信号特征与实际自适应信号模型特征,具有一定的理论价值与实际工程意义。
Description
技术领域
本发明涉及一种脑电信号去噪领域,特别涉及一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统。
背景技术
脑电图(Electroencephalogram,简称EEG)是神经信息学的重要组成部分,也是当今流行的无损伤探测大脑的高级技术。由于EEG信号中包括了丰富的神经信息和大量有用的生理、心理及疾病信息,可反应出大脑的病变,因此脑电图常被用作辅助诊断疾病。但由于大脑的神经系统特殊的互联结构,不同通道的EEG信号之间一定存在相互的噪声干扰,从而使得每一通道的EEG信号存在相互混迭和信息冗余的现象。另外在信号传输过程中,也必然会受到周围环境的干扰。
滤波是去除干扰的有效手段,根据信号的不同,可以选择不同的滤波器。传统的EEG去噪方法——小波变换,具有良好的时频特性,在去除含噪信号中随机噪声的同时,能最大限度保留信号中的突变特征成分。但该方法是以信号不变特性或统计特性平稳为前提,更适用于分析线性平稳信号;并且小波变换很大程度上取决于小波基的选择,其结果用于分析所有数据,使该方法的分解过程不能实现自适应。而脑电信号是一种随机性很强的非线性、非平稳信号,所以小波变换在脑电信号分析中会受到一定限制。
1960年卡尔曼提出了卡尔曼滤波算法,克服了维纳滤波器的缺点,首次将状态空间引入到信号处理中,更好地实现非平稳信号的滤波功能,因此卡尔曼滤波被广泛地应用。卡尔曼滤波是一种状态估计技术,它的操作是基于信号的状态空间表示。对于卡尔曼滤波来说不同的状态空间表示可能导致有不同的结构。卡尔曼滤波的一个很明显很重要的特点用是状态空间的方法来描述其数学推导公式。另一个比较创新的特点是,它的估计值是通过递归计算的方式得到的,即可应用于平稳系统也可以应用于非平稳的情况。
经验模态分解(Empirical Mode Decomposition,EMD)是Huang等人经过大量的研究分析,在1998年提出了一种新的自适应信号处理方法,不但适用于线性和平稳的序列,而且适用于非平稳的非线性的序列,从本身信号提取信息保证了局部信息的完整性。EMD方法在理论上可以实现对任意类型的非线性非平稳信号都进行分解处理,因此该方法一经提出后就迅速有效的应用在不同的工程领域,如应用于海洋、大气、以及各种机械故障诊断等方面。
针对传统去噪方法存在的问题,本发明提出一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法。该方法在有效信号得到最大程度保护的前提下消除尽可能多的噪声,并提高计算速度,以达到更好的去噪效果。
发明内容
本发明目的在于对现有研究和技术存在的不足加以完善与规范化,提出一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统。通过对脑电信号的经验模态分解可以保留脑电信号局部信息的完整性;通过卡尔曼滤波同时考虑对非平稳脑电信号的时变敏感最佳估计与节省数据计算的成本,具有一定的理论价值与实际应用意义。
本发明的目的通过以下的技术方案实现:第一方面,本发明提供了一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法,包括以下步骤:
(1)将电极在一段时间内采集的需去噪处理的脑电信号记作原始信号X(t),其中t=1,2,3,...,T;T为采样点数;
(2)初始化ri-1(t)=X(t),i=1;
(3)初始化Hj-1(t)=ri-1(t),j=1;
(4)对Hj-1(t)的极大极小值点进行三次样条插值,形成上下包络线,并计算出均值包络线Mj-1(t);
(5)令Hj(t)=Hj-1(t)-Mj-1(t);
(6)若Hj(t)是IMF分量,则求得本征模态函数IMFi=Hj(t),并进行步骤(7);否则,j=j+1,转到(4);
(7)令ri(t)=ri-1(t)-IMFi;
(8)如果ri(t)极值点仍多于2个,则i=i+1,转到(3);否则,分解结束,得到残余分量ri(t);令n=i,则原始信号X(t)用下式表示:
式中,X(t)表示原始信号,IMFi表示信号经过分解后的第i个IMF分量,且所有IMF分量频率按照从高到低的顺序进行排列,rn(t)表示残差分量;
(9)对于步骤(8)中的IMFi,通过计算与原信号相关系数进行筛选得到m个待去噪的信号IMFa,a=1,2,...,m;
(10)把一段长为t的信号视为一个t维系统,系统状态k表示为变量xk∈Rt,系统控制输入为uk,系统过程噪声为wk,则系统状态k由上一个状态k-1决定的随机差分方程为:
xk=Axk-1+Buk-1+wk-1
其中A表示在没有过程噪声下的状态转移矩阵,B是可选的控制输入uk的增益;wk-1为上一个状态的系统过程噪声;
定义观测变量z∈Rm,得到状态k下的观测方程为:
zk=Cxk+vk
其中C为观测矩阵,vk表示观测噪声,zk表示真实测量值;假设系统过程噪声和观测噪声是互相独立且正态分布的噪声,并且它们的协方差矩阵分别为Q和R;wk和vk的概率分布满足正态分布;
Gk是卡尔曼增益,由下式给出:
(12)对于步骤(9)中的m组IMF分量信号IMFa,令zka=IMFa,a=1,2,...,m,对于每个IMFa执行步骤(10)-步骤(12),最终得到向量X=[xk1 xk2 xk3 … xkm]H;xkm表示第m个IMF分量滤波后的结果;对X按列求和得到Yk0,作为经过卡尔曼滤波系统处理后的脑电信号;
(13)将(9)中未筛选的信号和Yk0加和,得到最终脑电信号Yk。
进一步地,步骤(4)中,对步骤(3)中的信号Hj-1(t)的局部上、下极值点,采用三次样条插值函数分别拟合形成原始信号的上、下包络线Uj-1(t)和Vj-1(t);求解上下包络线的均值Mj-1(t)=(Uj-1(t)+Vj-1(t))/2,得到该信号的均值包络线;从步骤(4)中的信号Hj-1(t)中减去均值Mj-1(t),得到中间信号Hj(t)=Hj-1(t)-Mj-1(t)。
进一步地,对步骤(1)中的原始信号X(k)基于信号本身特性进行经验模态分解时,分解结果应当是满足以下两个条件的本征模态函数(即IMF分量),即每个固有模态函数需满足:
1)整个信号段内的过零点数和极值点数相同或最多相差1;
2)由信号局部极大值点构造的上包络线和极小值点构造的下包络线的平均值在任意时刻均等于零,即上下包络线关于时间轴局部对称。
进一步地,对步骤(11)中的观测变量与预测得到的值之间的差值被称作测量过程中的革新或残余;实际值和预测值之间不一致程度是通过该差值来反映的;若该差值为0,则表明以上两者完全吻合;另外,卡尔曼滤波算法需要设置初始值和P(0);状态矢量的初始估计值设置为0或者是从先验信息得到的其他值,P(0)能够设置为单位矩阵。
进一步地,对于步骤(13)得到的最终处理后的脑电信号Yk通过计算和原信号X(k)的相关系数、计算信号Yk信噪比、计算信号Yk均方误差三种方法来衡量与原信号之间的去噪效果。
进一步地,对于步骤(9)中的各项IMFi的筛选方法采用计算与原信号X(k)相关系数的准则实现,具体实现方法如下:
1)给出任意两个离散信号相关系数R(X,Yi)的计算公式如下:
2)去掉残差以及相关系数较小的分量;按照手动设置的区间(A,B)来划分IMFi,具体参考如下的准则:
a.R<A,该分量相关性差,舍去该分量;
b.R∈[A,B],该分量相关性一般,将该分量加入IMFa中,进行卡尔曼滤波过程;
c.R>B,该分量相关性强,直接保留,不进行滤波处理;
3)根据实际信号分解结果,A为0.1,B为0.6,若处理结果不达标,在计算资源允许的情况下继续调整A和B的大小;
4)根据实际信号分类结果,IMFi的非线性与震荡性与经验模态分解次数S有关,具体表现为:
第二方面,本发明还提供了一种基于经验模态分解和卡尔曼滤波的脑电信号去噪装置,包括存储器和一个或多个处理器,所述存储器中存储有可执行代码,所述处理器执行所述可执行代码时,用于实现基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤。
第三方面,本发明还提供了一种计算机可读存储介质,其上存储有程序,该程序被处理器执行时,实现基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤。
第四方面,本发明还提供了一种基于经验模态分解和卡尔曼滤波的脑电信号去噪系统,该系统包括脑电信号采集模块、经验模态分解模块和卡尔曼滤波模块;
所述脑电信号采集模块用于采集脑电信号,将电极在一段时间内采集的需去噪处理的脑电信号记作原始信号X(t),其中t=1,2,3,...,T;T为采样点数;
所述经验模态分解模块用于基于脑电信号采集模块采集的原始信号,初始化ri-1(t)=X(t),i=1;初始化Hj-1(t)=ri-1(t),j=1;对Hj-1(t)的极大极小值点进行三次样条插值,形成上下包络线,并计算出均值包络线Mj-1(t);令Hj(t)=Hj-1(t)-Mj-1(t);并判断Hj(t)是否是IMF分量,若若Hj(t)不是IMF分量,令j=j+1,计算均值包络线,继续判断;若Hj(t)是IMF分量,则求得本征模态函数IMFi=Hj(t),令ri(t)=ri-1(t)-IMFi,判断ri(t)极值点是否仍多于2个,若是则i=i+1,继续经验模态分解;否则,分解结束,得到残余分量ri(t);令n=i,则原始信号X(t)用下式表示:
式中,X(t)表示原始信号,IMFi表示信号经过分解后的第i个IMF分量,且所有IMF分量频率按照从高到低的顺序进行排列,rn(t)表示残差分量;
所述卡尔曼滤波模块用于基于信号经过分解后的IMFi,通过计算与原信号相关系数进行筛选得到m个待去噪的信号IMFa,a=1,2,...,m;把一段长为t的信号视为一个t维系统,系统状态k表示为变量xk∈Rt,系统控制输入为uk,系统过程噪声为wk,则系统状态k由上一个状态k-1决定的随机差分方程为:
xk=Axk-1+Buk-1+wk-1
其中A表示在没有过程噪声下的状态转移矩阵,B是可选的控制输入uk的增益;wk-1为上一个状态的系统过程噪声;
定义观测变量z∈Rm,得到状态k下的观测方程为:
zk=Cxk+vk
其中C为观测矩阵,vk表示观测噪声,zk表示真实测量值;假设系统过程噪声和观测噪声是互相独立且正态分布的噪声,并且它们的协方差矩阵分别为Q和R;wk和vk的概率分布满足正态分布;
Gk是卡尔曼增益,由下式给出:
对于m组IMF分量信号IMFa,令zka=IMFa,a=1,2,...,m,对于每个IMFa,经过卡尔曼滤波后,最终得到向量X=[xk1 xk2 xk3 … xkm]H;xkm表示第m个IMF分量滤波后的结果;对X按列求和得到Yk0,作为经过卡尔曼滤波系统处理后的脑电信号;将未筛选的信号和Yk0加和,得到最终脑电信号Yk。
与现有技术相比,本发明具有以下创新优势及显著效果:
(1)针对脑电信号的非线性非平稳特性要求,选取经验模态分解方法对其进行自适应性分解,兼顾了脑电信号自身特征和数据运行特征,能尽可能保留了脑电信号的局部信息,也避免了小波分解的自适应问题。
(2)针对脑电信号的随机噪声形式要求,选取卡尔曼滤波方法对含有随机噪声的脑电信号进行理论上的最优递推最小方差估计,对随机噪声脑电信号有着较好的滤波效果。
(3)针对脑电信号处理过程中的计算量要求,提出了计算相关系数和IMF分量排序的双重筛选方法对脑电信号的各IMF分量进行筛选,减少了需要计算的信号组数目,极大降低了计算量。
(4)针对实际信号的非线性与震荡性与相关系数有关的特性和卡尔曼滤波系统要求,提出了计算相关系数和IMF分量排序的双重筛选方法对脑电信号的各IMF分量进行筛选,筛选出非线性较弱且相关系数较高的信号组继承给卡尔曼滤波处理,很好地克服了卡尔曼滤波在强非线性信号下的发散问题,提升了系统的全局收敛性。
附图说明
图1为待处理的脑电信号的图像。
图2和图3为脑电信号经过经验模态分解后各分量的图像。
图4为最终得到的去噪后的脑电信号图像。
图5为卡尔曼滤波流程图。
图6为一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的具体流程图。
图7为本发明提供的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪装置的结构图。
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图对本发明的具体实施方式做详细的说明。
在下面的描述中阐述了很多具体细节以便于充分理解本发明,但是本发明还可以采用其他不同于在此描述的其它方式来实施,本领域技术人员可以在不违背本发明内涵的情况下做类似推广,因此本发明不受下面公开的具体实施例的限制。
本发明实施例基于某医院某新生儿的的某段时间的一个电极接收到的连续脑电信号。该信号的采集时间在1min左右,共28000个离散点,波形如图1所示。
本实施例中需要去噪的脑电信号即为上述长度为28000个点的某新生儿的部分脑电数据,方法结果为得到最终处理后的脑电信号Yk。如图6所示,其详细实施步骤如下:
(1)根据脑电采集系统(电极)采集的相应需求周期中的脑电信号,将上述一个电极一段时间内记录的需去噪处理的脑电信号记作原始信号X(t),其中t=1,2,3,...,28000。对原始信号X(k)基于信号本身特性进行经验模态分解时,分解结果应当是满足以下两个条件的本征模态函数(即IMF分量),即每个固有模态函数需满足:1)整个信号段内的过零点数和极值点数相同或最多相差1;2)由信号局部极大值点构造的上包络线和极小值点构造的下包络线的平均值在任意时刻均等于零,即上下包络线关于时间轴局部对称。
(2)初始化ri-1(t)=X(t),i=1;
(3)初始化Hj-1(t)=ri-1(t),j=1;
(4)对Hj-1(t)的极大极小值点进行三次样条插值,形成上下包络线,并计算出均值包络线Mj-1(t);具体为:对步骤(3)中的信号Hj-1(t)的局部上、下极值点,采用三次样条插值函数分别拟合形成原始信号的上、下包络线Uj-1(t)和Vj-1(t);求解上下包络线的均值Mj-1(t)=(Uj-1(t)+Vj-1(t))/2,得到该信号的均值包络线;从步骤(4)中的信号Hj-1(t)中减去均值Mj-1(t),得到中间信号Hj(t)=Hj-1(t)-Mj-1(t)。
(5)令Hj(t)=Hj-1(t)-Mj-1(t);
(6)若Hj(t)是IMF分量,则求得本征模态函数IMFi=Hj(t),并进行步骤(7);否则,j=j+1,转到(4);
(7)令ri(t)=ri-1(t)-IMFi;
(8)如果ri(t)极值点仍多于2个,则i=i+1,转到(3);否则,分解结束,得到残余分量ri(t);令n=i,则原始信号X(t)用下式表示:
式中,X(t)表示原始信号,IMFi表示信号经过分解后的第i个IMF分量,且所有IMF分量频率按照从高到低的顺序进行排列,rn(t)表示残差分量;
(9)对于步骤(8)中的IMFi,通过计算与原信号相关系数进行筛选得到m个待去噪的信号IMFa,a=1,2,...,m;如图2和图3所示,为脑电信号经过经验模态分解后各分量的图像。各项IMFi的筛选方法采用计算与原信号X(k)相关系数的准则实现,具体实现方法如下:
1)给出任意两个离散信号相关系数R(X,Yi)的计算公式如下:
这里X(t)和Yi(t)表示待求取相关系数的两个信号,和代表信号X(t)和Yi(t)的平均值;在本实施例中,X(t)为原信号X(k),Yi(t)为各分量IMFi。计算R得到2×2实对称阵,选取R(1,2)作为指标R即可。具体计算数值如下:
从IMF1到IMF11,计算各IMF得到R分别如下所示:
第一个R的大小是:0.30810339792077146
第二个R的大小是:0.24811274559154803
第三个R的大小是:0.1943704693190584
第四个R的大小是:0.3098241917222679
第五个R的大小是:0.7012166374977818
第六个R的大小是:0.42366331808549323
第七个R的大小是:0.1915760853110749
第八个R的大小是:0.04258587245892567
第九个R的大小是:0.00063794930171923
第十个R的大小是:6.209173001462e^-05
第十一个R的大小是:0.0012534013607053
2)去掉残差以及相关系数较小的分量;按照手动设置的区间(A,B)来划分IMFi,具体参考如下的准则:
a.R<A,该分量相关性差,舍去该分量;
b.R∈[A,B],该分量相关性一般,将该分量加入IMFa中,进行卡尔曼滤波过程;
c.R>B,该分量相关性强,直接保留,不进行滤波处理;
3)根据实际信号分解结果,B一般不大于0.6;在本实施例中,A设置为0.1,B设置为0.6,若处理结果不达标,在计算资源允许的情况下继续调整A和B的大小;
4)根据实际信号分类结果,IMFi的非线性与震荡性与经验模态分解次数S有关,本实施例中,分解次数为11次,具体表现为:
该筛选方法有以下优点:
信号的非线性会影响卡尔曼滤波的效果,这是因为卡尔曼滤波最初是出自线性系统的一种方法,以贝叶斯概率模型为基础,并且要求信号的噪声满足高斯分布,在强非线性下会出现不收敛的情况,本方法选择非线性介于强弱之间的信号集合IMFi很好避免了这一点。
选取线性介于强弱之间的信号集合IMFi,只占所有信号IMFi的一部分,且信号波形简化,显著降低了计算量。
脑电噪声以随机噪声为主,使用卡尔曼滤波能很好适应滤波要求。经过筛选后得到IMFa,a=1,2,...,m。
(10)如图5所示,把一段长为t的信号视为一个t维系统,系统状态k表示为变量xk∈Rt,本发明中,x1设置为1,在下面迭代中,满足xk=xk-1,系统控制输入为uk,本发明方法处理的系统中无输入,故uk=0,系统过程噪声为wk,则系统状态k由上一个状态k-1决定的随机差分方程为:
xk=Axk-1+Buk-1+wk-1
其中A表示在没有过程噪声下的状态转移矩阵,B是可选的控制输入uk的增益;wk-1为上一个状态的系统过程噪声;
定义观测变量z∈Rm,得到状态k下的观测方程为:
zk=Cxk+vk
其中C为观测矩阵,它体现了状态xk和测量值zk的关系。vk表示观测噪声,zk表示真实测量值;假设系统过程噪声和观测噪声是互相独立且正态分布的噪声,并且它们的协方差矩阵分别为Q和R;wk和vk的概率分布P(W)和P(V)满足正态分布:
其中矩阵C和A与前述相同。观测变量与预测得到的值之间的差值被称作测量过程中的革新或残余;实际值和预测值之间不一致程度是通过该差值来反映的;若该差值为0,则表明以上两者完全吻合;另外,卡尔曼滤波算法需要设置初始值和P(0);状态矢量的初始估计值设置为0或者是从先验信息得到的其他值,P(0)能够设置为单位矩阵。
Gk是卡尔曼增益,由下式给出:
其中Q是过程激励噪声协方差矩阵,本发明中设置成长度与原脑电信号等长的随机噪声的协方差。是k-1状态下的先验误差协方差矩阵,Pk是k-1状态下后验误差协方差矩阵,并用来进一步计算k状态下的先验误差协方差矩阵I是单位矩阵。在以上各步的公式中H代表矩阵的转置;本发明中,测量值是一维信号,所以矩阵C和Gk是向量。R和Pk成为标量。本发明中,A设置为单位矩阵,P的初始值设置为10,在迭代过程中,满足:
P=P+Q
Gk=P/(P+R)
(12)对于步骤(9)中的m组IMF分量信号IMFa,令zka=IMFa,a=1,2,...,m,对于每个IMFa执行上述步骤(10)-步骤(12),最终得到向量X=[xk1 xk2 xk3 … xkm]H,xkm表示第m个IMF分量滤波后的结果;对X按列求和得到Yk0,作为经过卡尔曼滤波系统处理后的脑电信号。
(13)将(9)中未筛选的信号和Yk0加和,得到去噪后的最终脑电信号Yk,如图4所示。得到的最终处理后的脑电信号Yk通过计算和原信号X(k)的相关系数、计算信号Yk信噪比、计算信号Yk均方误差三种方法来衡量与原信号之间的去噪效果。
本发明整个实施例按照图6中所示的流程,对脑电信号数据进行处理并得到最终去噪后的信号。图1-图4为使用本发明基于经验模态分解和卡尔曼滤波的脑电信号去噪方法进行脑电信号去噪方法的实施与应用的各个环节结果。通过经验模态分解和卡尔曼滤波,保证了局部信息的完整性,也避免了非线性情况下的不稳定,还降低了计算量。因此,本方法顺利完成了对脑电信号的降噪流程,可以得到去噪后的脑电信号。本发明对后续的脑电信号工作等研究有指导性的意义。
另一方面,本发明还提供了一种基于经验模态分解和卡尔曼滤波的脑电信号去噪系统,该系统包括脑电信号采集模块、经验模态分解模块和卡尔曼滤波模块;
所述脑电信号采集模块用于采集脑电信号,将电极在一段时间内采集的需去噪处理的脑电信号记作原始信号X(t),其中t=1,2,3,...,T;T为采样点数;该模块具体实现过程参考本发明提供的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤描述。
所述经验模态分解模块用于基于脑电信号采集模块采集的原始信号,初始化ri-1(t)=X(t),i=1;初始化Hj-1(t)=ri-1(t),j=1;对Hj-1(t)的极大极小值点进行三次样条插值,形成上下包络线,并计算出均值包络线Mj-1(t);令Hj(t)=Hj-1(t)-Mj-1(t);并判断Hj(t)是否是IMF分量,若若Hj(t)不是IMF分量,令j=j+1,计算均值包络线,继续判断;若Hj(t)是IMF分量,则求得本征模态函数IMFi=Hj(t),令ri(t)=ri-1(t)-IMFi,判断ri(t)极值点是否仍多于2个,若是则i=i+1,继续经验模态分解;否则,分解结束,得到残余分量ri(t);令n=i,则原始信号X(t)用下式表示:
式中,X(t)表示原始信号,IMFi表示信号经过分解后的第i个IMF分量,且所有IMF分量频率按照从高到低的顺序进行排列,rn(t)表示残差分量;
经验模态分解模块具体实现过程参考本发明提供的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤描述。
所述卡尔曼滤波模块用于基于信号经过分解后的IMFi,通过计算与原信号相关系数进行筛选得到m个待去噪的信号IMFa,a=1,2,...,m;把一段长为t的信号视为一个t维系统,系统状态k表示为变量xk∈Rt,系统控制输入为uk,系统过程噪声为wk,则系统状态k由上一个状态k-1决定的随机差分方程为:
xk=Axk-1+Buk-1+wk-1
其中A表示在没有过程噪声下的状态转移矩阵,B是可选的控制输入uk的增益;wk-1为上一个状态的系统过程噪声;
定义观测变量z∈Rm,得到状态k下的观测方程为:
zk=Cxk+vk
其中C为观测矩阵,vk表示观测噪声,zk表示真实测量值;假设系统过程噪声和观测噪声是互相独立且正态分布的噪声,并且它们的协方差矩阵分别为Q和R;wk和vk的概率分布满足正态分布;
Gk是卡尔曼增益,由下式给出:
对于m组IMF分量信号IMFa,令zka=IMFa,a=1,2,...,m,对于每个IMFa,经过卡尔曼滤波后,最终得到向量X=[xk1 xk2 xk3 … xkm]H;xkm表示第m个IMF分量滤波后的结果;对X按列求和得到Yk0,作为经过卡尔曼滤波系统处理后的脑电信号;将未筛选的信号和Yk0加和,得到最终脑电信号Yk。
卡尔曼滤波模块具体实现过程参考本发明提供的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤描述。
与前述基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的实施例相对应,本发明还提供了基于经验模态分解和卡尔曼滤波的脑电信号去噪装置的实施例。
参见图7,本发明实施例提供的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪装置,包括存储器和一个或多个处理器,所述存储器中存储有可执行代码,所述处理器执行所述可执行代码时,用于实现上述实施例中的基于经验模态分解和卡尔曼滤波的脑电信号去噪方法。
本发明基于经验模态分解和卡尔曼滤波的脑电信号去噪装置的实施例可以应用在任意具备数据处理能力的设备上,该任意具备数据处理能力的设备可以为诸如计算机等设备或装置。装置实施例可以通过软件实现,也可以通过硬件或者软硬件结合的方式实现。以软件实现为例,作为一个逻辑意义上的装置,是通过其所在任意具备数据处理能力的设备的处理器将非易失性存储器中对应的计算机程序指令读取到内存中运行形成的。从硬件层面而言,如图7所示,为本发明基于经验模态分解和卡尔曼滤波的脑电信号去噪装置所在任意具备数据处理能力的设备的一种硬件结构图,除了图7所示的处理器、内存、网络接口、以及非易失性存储器之外,实施例中装置所在的任意具备数据处理能力的设备通常根据该任意具备数据处理能力的设备的实际功能,还可以包括其他硬件,对此不再赘述。
上述装置中各个单元的功能和作用的实现过程具体详见上述方法中对应步骤的实现过程,在此不再赘述。
对于装置实施例而言,由于其基本对应于方法实施例,所以相关之处参见方法实施例的部分说明即可。以上所描述的装置实施例仅仅是示意性的,其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的,作为单元显示的部件可以是或者也可以不是物理单元,即可以位于一个地方,或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本发明方案的目的。本领域普通技术人员在不付出创造性劳动的情况下,即可以理解并实施。
本发明实施例还提供一种计算机可读存储介质,其上存储有程序,该程序被处理器执行时,实现上述实施例中的基于经验模态分解和卡尔曼滤波的脑电信号去噪方法。
所述计算机可读存储介质可以是前述任一实施例所述的任意具备数据处理能力的设备的内部存储单元,例如硬盘或内存。所述计算机可读存储介质也可以是任意具备数据处理能力的设备的外部存储设备,例如所述设备上配备的插接式硬盘、智能存储卡(Smart Media Card,SMC)、SD卡、闪存卡(Flash Card)等。进一步的,所述计算机可读存储介质还可以既包括任意具备数据处理能力的设备的内部存储单元也包括外部存储设备。所述计算机可读存储介质用于存储所述计算机程序以及所述任意具备数据处理能力的设备所需的其他程序和数据,还可以用于暂时地存储已经输出或者将要输出的数据。
以上所述仅是本发明的优选实施方式,虽然本发明已以较佳实施例披露如上,然而并非用以限定本发明。任何熟悉本领域的技术人员,在不脱离本发明技术方案范围情况下,都可利用上述揭示的方法和技术内容对本发明技术方案做出许多可能的变动和修饰,或修改为等同变化的等效实施例。因此,凡是未脱离本发明技术方案的内容,依据本发明的技术实质对以上实施例所做的任何的简单修改、等同变化及修饰,均仍属于本发明技术方案保护的范围内。
Claims (9)
1.一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法,其特征在于,包括以下步骤:
(1)将电极在一段时间内采集的需去噪处理的脑电信号记作原始信号X(t),其中t=1,2,3,...,T;T为采样点数;
(2)初始化ri-1(t)=X(t),i=1;
(3)初始化Hj-1(t)=ri-1(t),j=1;
(4)对Hj-1(t)的极大极小值点进行三次样条插值,形成上下包络线,并计算出均值包络线Mj-1(t);
(5)令Hj(t)=Hj-1(t)-Mj-1(t);
(6)若Hj(t)是IMF分量,则求得本征模态函数IMFi=Hj(t),并进行步骤(7);否则,j=j+1,转到(4);
(7)令ri(t)=ri-1(t)-IMFi;
(8)如果ri(t)极值点仍多于2个,则i=i+1,转到(3);否则,分解结束,得到残余分量ri(t);令n=i,则原始信号X(t)用下式表示:
式中,X(t)表示原始信号,IMFi表示信号经过分解后的第i个IMF分量,且所有IMF分量频率按照从高到低的顺序进行排列,rn(t)表示残差分量;
(9)对于步骤(8)中的IMFi,通过计算与原信号相关系数进行筛选得到m个待去噪的信号IMFa,a=1,2,...,m;
(10)把一段长为t的信号视为一个t维系统,系统状态k表示为变量xk∈Rt,系统控制输入为uk,系统过程噪声为wk,则系统状态k由上一个状态k-1决定的随机差分方程为:
xk=Axk-1+Buk-1+wk-1
其中A表示在没有过程噪声下的状态转移矩阵,B是可选的控制输入uk的增益;wk-1为上一个状态的系统过程噪声;
定义观测变量z∈Rm,得到状态k下的观测方程为:
zk=Cxk+vk
其中C为观测矩阵,vk表示观测噪声,zk表示真实测量值;假设系统过程噪声和观测噪声是互相独立且正态分布的噪声,并且它们的协方差矩阵分别为Q和R;wk和vk的概率分布满足正态分布;
Gk是卡尔曼增益,由下式给出:
(12)对于步骤(9)中的m组IMF分量信号IMFa,令zka=IMFa,a=1,2,...,m,对于每个IMFa执行步骤(10)-步骤(12),最终得到向量X=[xk1 xk2 xk3...xkm]H;xkm表示第m个IMF分量滤波后的结果;对X按列求和得到Yk0,作为经过卡尔曼滤波系统处理后的脑电信号;
(13)将(9)中未筛选的信号和Yk0加和,得到最终脑电信号Yk。
2.根据权利要求1所述的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法,其特征在于,步骤(4)中,对步骤(3)中的信号Hj-1(t)的局部上、下极值点,采用三次样条插值函数分别拟合形成原始信号的上、下包络线Uj-1(t)和Vj-1(t);求解上下包络线的均值Mj-1(t)=(Uj-1(t)+Vj-1(t))/2,得到该信号的均值包络线;从步骤(4)中的信号Hj-1(t)中减去均值Mj-1(t),得到中间信号Hj(t)=Hj-1(t)-Mj-1(t)。
3.根据权利要求1所述的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法,其特征在于,对步骤(1)中的原始信号X(k)基于信号本身特性进行经验模态分解时,分解结果应当是满足以下两个条件的本征模态函数(即IMF分量),即每个固有模态函数需满足:
1)整个信号段内的过零点数和极值点数相同或最多相差1;
2)由信号局部极大值点构造的上包络线和极小值点构造的下包络线的平均值在任意时刻均等于零,即上下包络线关于时间轴局部对称。
5.根据权利要求1所述的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法,其特征在于,对于步骤(13)得到的最终处理后的脑电信号Yk通过计算和原信号X(k)的相关系数、计算信号Yk信噪比、计算信号Yk均方误差三种方法来衡量与原信号之间的去噪效果。
6.根据权利要求1所述的一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法,其特征在于,对于步骤(9)中的各项IMFi的筛选方法采用计算与原信号X(k)相关系数的准则实现,具体实现方法如下:
1)给出任意两个离散信号相关系数R(X,Yi)的计算公式如下:
2)去掉残差以及相关系数较小的分量;按照手动设置的区间(A,B)来划分IMFi,具体参考如下的准则:
a.R<A,该分量相关性差,舍去该分量;
b.R∈[A,B],该分量相关性一般,将该分量加入IMFa中,进行卡尔曼滤波过程;
c.R>B,该分量相关性强,直接保留,不进行滤波处理;
3)根据实际信号分解结果,A为0.1,B为0.6,若处理结果不达标,在计算资源允许的情况下继续调整A和B的大小;
4)根据实际信号分类结果,IMFi的非线性与震荡性与经验模态分解次数S有关,具体表现为:
7.一种基于经验模态分解和卡尔曼滤波的脑电信号去噪装置,包括存储器和一个或多个处理器,所述存储器中存储有可执行代码,其特征在于,所述处理器执行所述可执行代码时,用于实现如权利要求1-6中任一项所述的基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤。
8.一种计算机可读存储介质,其上存储有程序,其特征在于,该程序被处理器执行时,实现如权利要求1-6中任一项所述的基于经验模态分解和卡尔曼滤波的脑电信号去噪方法的步骤。
9.一种实现权利要求1-6任一项所述方法的基于经验模态分解和卡尔曼滤波的脑电信号去噪系统,其特征在于,该系统包括脑电信号采集模块、经验模态分解模块和卡尔曼滤波模块;
所述脑电信号采集模块用于采集脑电信号,将电极在一段时间内采集的需去噪处理的脑电信号记作原始信号X(t),其中t=1,2,3,...,T;T为采样点数;
所述经验模态分解模块用于基于脑电信号采集模块采集的原始信号,初始化ri-1(t)=X(t),i=1;初始化Hj-1(t)=ri-1(t),j=1;对Hj-1(t)的极大极小值点进行三次样条插值,形成上下包络线,并计算出均值包络线Mj-1(t);令Hj(t)=Hj-1(t)-Mj-1(t);并判断Hj(t)是否是IMF分量,若若Hj(t)不是IMF分量,令j=j+1,计算均值包络线,继续判断;若Hj(t)是IMF分量,则求得本征模态函数IMFi=Hj(t),令ri(t)=ri-1(t)-IMFi,判断ri(t)极值点是否仍多于2个,若是则i=i+1,继续经验模态分解;否则,分解结束,得到残余分量ri(t);令n=i,则原始信号X(t)用下式表示:
式中,X(t)表示原始信号,IMFi表示信号经过分解后的第i个IMF分量,且所有IMF分量频率按照从高到低的顺序进行排列,rn(t)表示残差分量;
所述卡尔曼滤波模块用于基于信号经过分解后的IMFi,通过计算与原信号相关系数进行筛选得到m个待去噪的信号IMFa,a=1,2,...,m;把一段长为t的信号视为一个t维系统,系统状态k表示为变量xk∈Rt,系统控制输入为uk,系统过程噪声为wk,则系统状态k由上一个状态k-1决定的随机差分方程为:
xk=Axk-1+Buk-1+wk-1
其中A表示在没有过程噪声下的状态转移矩阵,B是可选的控制输入uk的增益;wk-1为上一个状态的系统过程噪声;
定义观测变量z∈Rm,得到状态k下的观测方程为:
zk=Cxk+vk
其中C为观测矩阵,vk表示观测噪声,zk表示真实测量值;假设系统过程噪声和观测噪声是互相独立且正态分布的噪声,并且它们的协方差矩阵分别为Q和R;wk和vk的概率分布满足正态分布;
Gk是卡尔曼增益,由下式给出:
对于m组IMF分量信号IMFa,令zka=IMFa,a=1,2,...,m,对于每个IMFa,经过卡尔曼滤波后,最终得到向量X=[xk1 xk2 xk3...xkm]H;xkm表示第m个IMF分量滤波后的结果;对X按列求和得到Yk0,作为经过卡尔曼滤波系统处理后的脑电信号;将未筛选的信号和Yk0加和,得到最终脑电信号Yk。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211015305.0A CN115211869A (zh) | 2022-08-23 | 2022-08-23 | 一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211015305.0A CN115211869A (zh) | 2022-08-23 | 2022-08-23 | 一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN115211869A true CN115211869A (zh) | 2022-10-21 |
Family
ID=83615104
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211015305.0A Pending CN115211869A (zh) | 2022-08-23 | 2022-08-23 | 一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115211869A (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117093827A (zh) * | 2023-10-16 | 2023-11-21 | 欣灵电气股份有限公司 | 基于物联网的智慧消防给水数据处理系统 |
CN117310455A (zh) * | 2023-11-30 | 2023-12-29 | 上海采起电子技术服务有限公司 | 一种基于电信号的口腔扫描仪电路板故障检测方法 |
CN117349603A (zh) * | 2023-12-06 | 2024-01-05 | 小舟科技有限公司 | 脑电信号的自适应降噪方法及装置、设备、存储介质 |
CN117648536A (zh) * | 2024-01-26 | 2024-03-05 | 浙江威利坚科技股份有限公司 | 一种空调专用电弧信号的去噪方法 |
-
2022
- 2022-08-23 CN CN202211015305.0A patent/CN115211869A/zh active Pending
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117093827A (zh) * | 2023-10-16 | 2023-11-21 | 欣灵电气股份有限公司 | 基于物联网的智慧消防给水数据处理系统 |
CN117093827B (zh) * | 2023-10-16 | 2024-01-30 | 欣灵电气股份有限公司 | 基于物联网的智慧消防给水数据处理系统 |
CN117310455A (zh) * | 2023-11-30 | 2023-12-29 | 上海采起电子技术服务有限公司 | 一种基于电信号的口腔扫描仪电路板故障检测方法 |
CN117310455B (zh) * | 2023-11-30 | 2024-02-09 | 上海采起电子技术服务有限公司 | 一种基于电信号的口腔扫描仪电路板故障检测方法 |
CN117349603A (zh) * | 2023-12-06 | 2024-01-05 | 小舟科技有限公司 | 脑电信号的自适应降噪方法及装置、设备、存储介质 |
CN117349603B (zh) * | 2023-12-06 | 2024-03-12 | 小舟科技有限公司 | 脑电信号的自适应降噪方法及装置、设备、存储介质 |
CN117648536A (zh) * | 2024-01-26 | 2024-03-05 | 浙江威利坚科技股份有限公司 | 一种空调专用电弧信号的去噪方法 |
CN117648536B (zh) * | 2024-01-26 | 2024-04-19 | 浙江威利坚科技股份有限公司 | 一种空调专用电弧信号的去噪方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN115211869A (zh) | 一种基于经验模态分解和卡尔曼滤波的脑电信号去噪方法、装置及系统 | |
Singh et al. | Denoising of ECG signal by non-local estimation of approximation coefficients in DWT | |
CN109165556B (zh) | 一种基于grnn身份识别方法 | |
Tracey et al. | Nonlocal means denoising of ECG signals | |
Smital et al. | Adaptive wavelet wiener filtering of ECG signals | |
Luo et al. | A hierarchical method for removal of baseline drift from biomedical signals: application in ECG analysis | |
Mourad | ECG denoising algorithm based on group sparsity and singular spectrum analysis | |
CN114224360B (zh) | 一种基于改进emd-ica的eeg信号处理方法、设备及存储介质 | |
CN110680308A (zh) | 基于改进emd与阈值法融合的心电信号去噪方法 | |
WO2020228420A1 (zh) | 降噪自编码器的训练方法、心电信号的降噪方法及其装置 | |
Phadikar et al. | Automatic EEG eyeblink artefact identification and removal technique using independent component analysis in combination with support vector machines and denoising autoencoder | |
Madan et al. | Denoising of ECG signals using weighted stationary wavelet total variation | |
Sheoran et al. | A new method for automatic electrooculogram and eye blink artifacts correction of EEG signals using CCA and NAPCT | |
Mir et al. | ECG denoising and feature extraction techniques–a review | |
Vargas et al. | Electrocardiogram signal denoising by clustering and soft thresholding | |
Fatemi et al. | An online subspace denoising algorithm for maternal ECG removal from fetal ECG signals | |
CN114886378A (zh) | 一种基于改进的互补集合模态分解联合去噪方法及系统 | |
Vargas et al. | Empirical mode decomposition, viterbi and wavelets applied to electrocardiogram noise removal | |
Prateek et al. | Sparsity-assisted signal denoising and pattern recognition in time-series data | |
Sardouie et al. | Interictal EEG noise cancellation: GEVD and DSS based approaches versus ICA and DCCA based methods | |
Walters-Williams et al. | BMICA-independent component analysis based on B-spline mutual information estimation for EEG signals | |
CN107895387B (zh) | Mri图像重建方法及装置 | |
Zhong et al. | Non-invasive fetal electrocardiography denoising using deep convolutional encoder-decoder networks | |
Chang | A comparative analysis of various respiratory sound denoising methods | |
Mourad | ECG denoising based on 1-D double-density complex DWT and SBWT |
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 |