CN112294340A - 快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备 - Google Patents
快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备 Download PDFInfo
- Publication number
- CN112294340A CN112294340A CN202011164260.4A CN202011164260A CN112294340A CN 112294340 A CN112294340 A CN 112294340A CN 202011164260 A CN202011164260 A CN 202011164260A CN 112294340 A CN112294340 A CN 112294340A
- Authority
- CN
- China
- Prior art keywords
- electroencephalogram
- mimf
- artifacts
- imf
- artifact
- 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
- 238000000034 method Methods 0.000 title claims abstract description 224
- 230000003183 myoelectrical effect Effects 0.000 title claims abstract description 79
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 40
- 241000932075 Priacanthus hamrur Species 0.000 claims abstract description 19
- 230000003044 adaptive effect Effects 0.000 claims abstract description 15
- 238000004364 calculation method Methods 0.000 claims abstract description 13
- 230000006870 function Effects 0.000 claims description 38
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 238000004590 computer program Methods 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 7
- 239000013598 vector Substances 0.000 claims description 7
- 239000002131 composite material Substances 0.000 claims description 6
- 238000012935 Averaging Methods 0.000 claims description 4
- 238000005259 measurement Methods 0.000 claims description 3
- 230000001105 regulatory effect Effects 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 9
- 238000000926 separation method Methods 0.000 abstract description 7
- 238000001228 spectrum Methods 0.000 abstract description 6
- 238000004088 simulation Methods 0.000 description 31
- 238000000537 electroencephalography Methods 0.000 description 18
- 238000012545 processing Methods 0.000 description 15
- 210000004556 brain Anatomy 0.000 description 13
- 238000002474 experimental method Methods 0.000 description 13
- 238000001914 filtration Methods 0.000 description 12
- 230000003595 spectral effect Effects 0.000 description 11
- 210000003205 muscle Anatomy 0.000 description 7
- 238000004458 analytical method Methods 0.000 description 6
- 238000011156 evaluation Methods 0.000 description 5
- 230000005611 electricity Effects 0.000 description 4
- 238000012880 independent component analysis Methods 0.000 description 4
- 230000014759 maintenance of location Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 210000004761 scalp Anatomy 0.000 description 3
- 241000277331 Salmonidae Species 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000007177 brain activity Effects 0.000 description 2
- 230000008602 contraction Effects 0.000 description 2
- 238000010219 correlation analysis Methods 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 210000004709 eyebrow Anatomy 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 230000000284 resting effect Effects 0.000 description 2
- 230000033764 rhythmic process Effects 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 229910021607 Silver chloride Inorganic materials 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000002001 electrophysiology Methods 0.000 description 1
- 230000007831 electrophysiology Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000005764 inhibitory process Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- HKZLPVFGJNLROG-UHFFFAOYSA-M silver monochloride Chemical compound [Cl-].[Ag+] HKZLPVFGJNLROG-UHFFFAOYSA-M 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
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/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Biophysics (AREA)
- Heart & Thoracic Surgery (AREA)
- Physiology (AREA)
- Psychiatry (AREA)
- Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Power Engineering (AREA)
- Pathology (AREA)
- Biomedical Technology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了一种快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备，本发明的方法首先利用最小二乘多维经验模态分解方法对脑电信号进行时频域自适应分解，得到多维一致本征模态函数；然后针对每个本征模态函数中携带的肌电伪迹强度不同，利用SCAD区间阈值化方法自适应计算每个本征模态函数去除肌电伪迹的阈值，再将本征模态函数和阈值同时输入至SCAD区间阈值函数中得到无伪迹本征模态函数；最后重构得到无肌电伪迹的脑电。本发明避免了源分离方法中的计算误差，同时减小了部分重构方法误去除脑电导致的信号失真问题；解决了目前方法自适应性较差，导联数量过少去除效果不佳以及脑电与肌电伪迹频谱严重重叠时脑电误去除的问题。
Description
技术领域
本发明涉及电生理学技术领域，特别涉及一种快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备。
背景技术
脑电图(Electroencephalography,EEG)是一种电生理学技术，通过在头皮上放置电极来记录大脑的电活动。脑电图因其具有时间分辨率高、无侵入性、成本低、适用于长期监测等特点，在医学和科学领域有着广泛的应用，包括诊断、分析疾病和脑机接口等。近年来，可穿戴便携式的脑电设备广泛应用，用于随时随地的健康监测。这种移动脑电设备通常导联数量比较少，方便实时使用。但脑电信号经常受到大脑活动以外的各种生理信号影响，其中，肌电伪迹(Electromygraphic，EMG)取决于肌肉的收缩程度和收缩类型，具有幅值大、频域分布广等特征，导致肌电伪迹难以去除。
目前去除脑电中肌电伪迹的方法主要有滤波、盲源分离(Blind sourcesseparation,BSS)和信号分解。在滤波方法中，经典带通、低通等滤波方法在应用于脑电等非平稳信号时，分解不准确，为了克服这一缺陷出现了自适应滤波方法，自适应滤波方法假设脑电信号和肌电伪迹是不相关的，利用参考信号产生与肌电伪迹相关的估计信号，然后从脑电信号中减去估计信号。在去除肌电伪迹时，由于肌电伪迹复杂的特征，自适应滤波的实现需要大量的肌电伪迹参考信号来准确估计肌电伪迹，这增加了硬件设备的复杂性。
BSS方法将脑电信号与肌电伪迹分离成多个分量，然后去除与肌电伪迹相关的分量。许多BSS方法已经应用到去除肌电伪迹的研究上，例如独立成分分析(Independentcomponent analysis,ICA)、典型相关分析(Canonical correlation analysis,CCA)和独立向量分析(Independent vector analysis,IVA)等。BSS方法使用时，通常假设隐藏源的数量小于或等于导联数量，这使得BSS方法只有在足够数量导联时才有效，所以BBS方法受到导联数量限制，在少导联场景下使用效果不佳。
为了实现有效的在少导联场景下去除脑电中肌电伪迹，一些学者提出先利用信号分解方法将脑电分解成多个分量，再对分量进行细致的处理的策略。信号分解方法主要有小波分解(Wavelet decomposition,WD)[11]和经验模态分解(Empirical modedecomposition，EMD)。WD方法将信号分解为小波系数，将与肌电伪迹相关的小波系数置零，最后利用逆小波重构得到去除肌电伪迹信号。但脑电和肌电伪迹准确分离依赖于小波基的选择和分解层数，难以在实际中使用。EMD类方法是数据驱动的信号分解方法，它将信号自适应地分解为若干AM/FM零均值信号，称为本征模态函数(Intrinsic mode function,IMF)。基于EMD类去除少导联脑电中肌电伪迹的方法，根据各方法对IMF处理方式不同可分为三类：(1)IMF部分重构方法，Teng等[1]使用多维经验模态分解(Multivariate EMD，MEMD)对6导联脑电进行分解得到MIMF，利用样本熵来判定肌电IMF并去除，但这种直接将IMF去除的方法会导致重构信号丢失脑电信息。(2)IMF源分离方法，即联合EMD类方法和BSS方法对IMF信息再挖掘，将所有IMF由ICA、CCA等BSS方法处理计算隐藏源(包括脑电和肌电等隐藏源)，将自相关性低的隐藏源设置为肌电伪迹并去除，例如，集成经验模态分解典型相关分析(Ensemble EMD-CCA,EEMD-CCA)[2]，多变量EEMD-CCA[3]，MEMD-CCA[4]和最小二乘多维经验模态分解典型相关分析(Least square MEMD-CCA,LSMEMD-CCA)[5]等是常见有效的少导联脑电中肌电伪迹去除方法。但该类方法在利用自相关系数判定肌电伪迹时缺乏自适应性，阈值大小对伪迹去除效果有很大影响；另外，这类方法在利用逆BSS计算估计信号时，会产生不可避免的计算误差。(3)IMF自适应阈值滤波方法，基于EMD的区间阈值化(Intervalthresholding，IT)[6-7]方法通过对分解出的IMF的过零区间进行自适应阈值处理去除肌电伪迹，具有很好的自适应性。Safieddine等[8]通过实验证明EMD-IT方法性能优于BSS方法，但该方法在处理多维信号时，使用一维信号分解方法对导联逐一处理，未考虑脑电信号整体性和导联间的关系，无法得到理想的效果。
综上所述，基于EMD类去除少导联脑电中肌电伪迹是行之有效的方法，但如上文所述，依然存在缺陷。所以现在急需一种有效、具有很好自适应性的少导联脑电中肌电伪迹去除方法。
引用文献：
1 Teng C,Zhang Y,Wang G.The removal of EMG artifact from EEG signalsby the multivariate empirical mode decomposition.In:2014IEEE InternationalConference on Signal Processing,Communications and Computing.Guilin,China:IEEE,2014.873-876.
2 Sweeney K T,Mcloone S,Ward T E.The Use of Ensemble Empirical ModeDecomposition With Canonical Correlation Analysis as a Novel Artifact RemovalTechnique.IEEE Transactions on Biomedical Engineering,2013,60(1):97-105.
3 Chen X,Chen Q,Zhang Y,Wang Z J.A Novel EEMD-CCA Approach toRemoving Muscle Artifacts for Pervasive EEG.IEEE Sensors Journal,2019,19(19):8420-8431.
4 Chen X,Xu X,Liu A,McKeown M J,Wang Z J.The use of multivariate EMDand CCA for denoising muscle artifacts from few-channel EEG recordings.IEEETransactions on Instrumentation and Measurement,2017,67(2):359-370.
5 Liu Y,Zhou Y,Lang X,Lang X,Liu Y,Zheng Q,et al.An efficient androbust muscle artifact removal method for few-channel EEG.IEEE Access,2019,7:176036-176050.
6 Kopsinis Y,Mclaughlin S.Development of EMD-based denoising methodsinspired by wavelet thresholding.IEEE Transactions on Signal Processing,2009,57(4):1351-1362.
7 Kopsinis Y,Mclaughlin S.Empirical mode decomposition based soft-thresholding.In:European Signal Processing Conference.Lausanne,Switzerland:IEEE,2008.1-5.
8 Safieddine D,Kachenoura A,Albera L,Birot G,Karfoul A,Pasnicu A,etal.Removal of muscle artifact from EEG data:comparison between stochastic(ICAand CCA)and deterministic(EMD and wavelet-based)approaches.EURASIP Journal onAdvances in Signal Processing,2012,2012(1):127.
9 Welch P D.The use of fast Fourier transform for the estimation ofpower spectra:A method based on time averaging over short,modifiedperiodograms.IEEE Transactions on Audio and Electroacoustics,1967,15(2):70-73.
10 Mandic D P.EEG data[Online],available:http://www.commsp.ee.ic.ac.uk/～mandic/re-search/BSS_Stuff.htm,Sept 26,2019.
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备。本发明提出了一种快速自动去除肌电伪迹的方法，首先利用最小二乘多维经验模态分解方法对脑电信号进行时频域自适应分解，得到多维一致本征模态函数；然后针对每个本征模态函数中携带的肌电伪迹强度不同，利用SCAD区间阈值化方法自适应计算每个本征模态函数去除肌电伪迹的阈值，再将本征模态函数和阈值同时输入至SCAD区间阈值函数中得到无伪迹本征模态函数；最后将无肌电伪迹本征模态函数重构，得到无肌电伪迹的脑电。
为实现上述目的，本发明采用的技术方案是：提供一种快速自动去除肌电伪迹的方法，包括以下步骤：
S1、MIMF求解：将N通道含有肌电伪迹的真实脑电和M通道的辅助噪声构成(N+M)维复合信号，使用最小二乘多维经验模态分解方法将该(N+M)维复合信号分解为(N+M)个MIMF集合，去除辅助噪声分解得到的MIMF，只保留脑电信号分解出的N个MIMF集合；
S2、肌电伪迹去除：针对基于所述步骤S1所得的N个MIMF集合，使用SCAD区间阈值化方法，根据IMF能量和特性自适应计算N个MIMF集合中的每个IMF的阈值，联合阈值函数对MIMF中每个过零区间进行肌电伪迹阈值去除，得到无肌电伪迹的MIMF；其中，IMF表示本征模态函数，MIMF表示多维本征模态函数，即1个MIMF集合中包括多个IMF；
S3、信号重构：对每一个通道的无肌电伪迹MIMF进行重构，得到N维不含肌电伪迹的脑电。
优选的是，所述步骤S1中采用最小二乘多维经验模态分解方法对脑电信号进行自适应的时频分解，其具体步骤包括：
1-1)假设X(t)是一个N通道脑电信号，t＝1,…,T，Y(t)是一个M通道辅助噪声信号，M是任意的正整数，输入信号Z(t)＝[X(t),Y(t)]；
1-2)设i＝1，并产生L组用于均匀投影的超球面向量点集；
1-3)计算输入信号Z(t)的第l个投影函数fl(t)，计算方式为：
fl(t)＝Z(t)·(Vl(t))T (2)；
1-4)对于所有l，利用EMD算法提取相应投影函数fl(t)的第一层单变量IMFul(t)；
1-5)联合所有IMF函数
以及它们的投影向量
则相应的(N+M)维MIMF函数Di(t)可以通过解下列超定线性方程组获得
1-6)计算残差函数：
S(t)＝Z(t)-Ui(t) (4)；
1-7)如果S(t)没有足够的极值组成包络，停止当前算法并获取最终残差项R(t)＝S(t)；反之则更新当前输入为Z(t)＝S(t)及i＝i+1，然后转至步骤1-3)；
1-8)在停止上述筛选过程后，获得Q组(N+M)维
和残差R(t)，信号Z(t)表示为：
优选的是，所述步骤S2中，采用基于SCAD区间阈值化的方法进行MIMF肌电去除，具体包括以下步骤：
2-1)对于每一个MIMF集合，其中包括Q个IMF，对其中的第2～Q个IMF进行重组，得到
2-2)随机改变第1个IMF采样点的位置，得到
2-3)构建一个不同于原始信号的新的含肌电伪迹信号，得到
2-4)对新信号Xa(t)做最小二乘多维经验模态分解，得到
2-5)对步骤2-4)分解出的MIMF使用区间阈值化方法去除肌电伪迹，重构得到信号
2-6)将步骤2-2)至2-5)迭代K-1次，得到K个去除肌电伪迹信号：
2-7)将步骤2-6)得到的结果平均，获得最后的无肌电伪迹的脑电信号，
优选的是，所述步骤2-5)中使用区间阈值化方法去除肌电伪迹的步骤包括：
2-5-1)阈值计算：
对于每一个MIMF集合中的任意一个IMF计算其自适应阈值Ti：
其中C是一个调节系数，Ei为第i个IMF能量，T是采样长度；
先利用下式(7)所示的估计器计算第一个IMF的能量E1，然后利用式(8)计算其他IMF的能量Ei，u1(t)为第一个IMF：
其中，式(8)中β和ρ为设定的参数；
2-5-2)基于SCAD区间阈值函数的肌电伪迹去除：
SCAD区间阈值函数为：
其中，第i个IMF中第j和j+1个过零点构成的过零区间，
为第i个IMF中区间内所有测量值，
为第i个IMF中的极值，Ti为第i个IMF的阈值，
为去除肌电伪迹后的估计值；其中α为调节因子。
优选的是，所述步骤2-6)中，K＝9。
优选的是，步骤2-5-1)中，β为0.719，ρ为2.01。
优选的是，步骤2-5-2)中，α为3.7。
本发明还提供一种快速自动去除肌电伪迹的系统，该系统采用如上所述的方法去除脑电中的肌电伪迹。
本发明还提供一种存储介质，其上存储有计算机程序，该程序被执行时用于实现如上所述的方法。
本发明还提供一种计算机设备，包括存储器、处理器以及存储在所述存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上所述的方法。
本发明的有益效果是：
本发明提供了一种自适应、有效的少导联脑电肌电伪迹去除方法，本发明利用LSMEMD对少导联脑电进行自适应时频同步分解，符合导联间的同步和一致性；本发明采用的SCAD区间阈值化是一种符合IMF特性的自适应阈值计算方法，对LSMEMD分解产物MIMF进行更加细致的处理，避免了源分离方法中的计算误差，同时减小了部分重构方法误去除脑电导致的信号失真问题；本发明的方法基于LSMEMD和区间阈值进行改进，解决了目前方法自适应性较差，导联数量过少去除效果不佳以及脑电与肌电伪迹频谱严重重叠时脑电误去除的问题；通过仿真和真实数据实验，与传统的BSS方法CCA，信号分解部分重构方法LSMEMD-PR和目前较有效的少导联去除肌电伪迹方法LSMEMD-CCA进行比较，进一步证明本发明提供的LSMEMD-IT方法优于现有方法的性能。
附图说明
图1为本发明的实施例1中的快速自动去除肌电伪迹的方法的原理框图；
图2为本发明的实施例1中的区间阈值化算法的流程图；
图3为本发明的实施例2中的仿真信号的示意图；
图4为本发明的实施例2中的受肌电伪迹影响的真实数据；
图5a为本发明的实施例2中的RRMSE仿真性能比较结果；
图5b为本发明的实施例2中的PSNR仿真性能比较结果；
图5c为本发明的实施例2中的去噪后信噪比仿真性能比较结果；
图5d为本发明的实施例2中的ACC仿真性能比较结果；
图6为本发明的实施例2中的仿真数据隐藏源和LSMEMD-CCA处理后的隐藏源；
图7a为本发明的实施例2中的3个导联SNR＝0.5时不同方法去除肌电伪迹结果；
图7b为本发明的实施例2中的图7a中的3个导联中的第2个导联使用不同方法去除肌电伪迹后的IMF；
图7c为本发明的实施例2中的图7a中的3个导联中的第2个导联使用不同方法去除肌电伪迹后IMF的功率谱密度；
图7d为本发明的实施例2中的3导联SNR＝0.5时不同方法去除肌电伪迹后功率谱密度；
图8a为本发明的实施例2中的6个导联SNR＝0.5时不同方法去除肌电伪迹结果；
图8b为本发明的实施例2中的图8a中的6个导联中的第2个导联使用不同方法去除肌电伪迹后的IMF；
图8c为本发明的实施例2中的图8a中的6个导联中的第2个导联使用不同方法去除肌电伪迹后IMF的功率谱密度；
图8d为本发明的实施例2中的6导联SNR＝0.5时不同方法去除肌电伪迹后功率谱密度；
图9为本发明的实施例2中的真实数据源采用不同方法的自相关系数；
图10a为本发明的实施例2中的真实数据源1采用不同方法的功率谱密度；
图10b为本发明的实施例2中的真实数据源1采用LSMEMD-IT方法的伪迹去除结果；
图11a为本发明的实施例2中的真实数据源2采用不同方法的功率谱密度；
图11b为本发明的实施例2中的真实数据源2采用LSMEMD-IT方法的伪迹去除结果；
图12a为本发明的实施例2中的真实数据源3采用不同方法的功率谱密度；
图12b为本发明的实施例2中的真实数据源3采用LSMEMD-IT方法的伪迹去除结果。
具体实施方式
下面结合实施例对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不排除一个或多个其它元件或其组合的存在或添加。
实施例1
参照图1，本实施例的一种快速自动去除肌电伪迹的方法，包括以下步骤：
S1、MIMF求解：将N通道含有肌电伪迹的真实脑电和M通道的辅助噪声构成(N+M)维复合信号，使用最小二乘多维经验模态分解方法将该(N+M)维复合信号分解为(N+M)个MIMF集合，去除辅助噪声分解得到的MIMF，只保留脑电信号分解出的N个MIMF集合；
S2、肌电伪迹去除：针对基于所述步骤S1所得的N个MIMF集合，使用SCAD区间阈值化方法，根据IMF能量和特性自适应计算N个MIMF集合中的每个IMF的阈值，联合阈值函数对MIMF中每个过零区间进行肌电伪迹阈值去除，得到无肌电伪迹的MIMF；其中，IMF表示本征模态函数，MIMF表示多维本征模态函数，即1个MIMF集合中包括多个IMF；
S3、信号重构：对每一个通道的无肌电伪迹MIMF进行重构，得到N维不含肌电伪迹的脑电。
其中，所述步骤S1中采用LSMEMD方法(最小二乘多维经验模态分解方法)对脑电信号进行自适应的时频分解，其具体步骤包括：
1-1)假设X(t)是一个N通道脑电信号，t＝1,…,T，Y(t)是一个M通道辅助噪声(高斯白噪声)信号，M是任意的正整数，输入信号Z(t)＝[X(t),Y(t)]；通过平衡分解性能和计算量，在一种优选的实施例中，设置5个辅助噪声通道；
1-2)设i＝1，并产生L组用于均匀投影的超球面向量点集；
1-3)计算输入信号Z(t)的第l个投影函数fl(t)，计算方式为：
fl(t)＝Z(t)·(Vl(t))T (2)；
1-4)对于所有l，利用EMD算法提取相应投影函数fl(t)的第一层单变量IMFul(t)；
1-5)联合所有IMF函数
以及它们的投影向量
则相应的(N+M)维MIMF函数Di(t)可以通过解下列超定线性方程组获得
1-6)计算残差函数：
S(t)＝Z(t)-Ui(t) (4)；
1-7)如果S(t)没有足够的极值组成包络，停止当前算法并获取最终残差项R(t)＝S(t)；反之则更新当前输入为Z(t)＝S(t)及i＝i+1，然后转至步骤1-3)；
1-8)在停止上述筛选过程后，获得Q组(N+M)维
和残差R(t)，信号Z(t)表示为：
其中，所述步骤S2中，采用基于SCAD区间阈值化的方法进行MIMF肌电去除，参照图2，具体包括以下步骤：
2-1)对于每一个MIMF集合，其中包括Q个IMF，对其中的第2～Q个IMF进行重组，得到
2-2)随机改变第1个IMF采样点的位置，得到
2-3)构建一个不同于原始信号的新的含肌电伪迹信号，得到
2-4)对新信号Xa(t)做最小二乘多维经验模态分解，得到
2-5)对步骤2-4)分解出的MIMF使用区间阈值化方法去除肌电伪迹，重构得到信号
2-6)将步骤2-2)至2-5)迭代K-1次，得到K个去除肌电伪迹信号：
考虑到脑电信号的非平稳，不规则性，在一种优选的实施例中，迭代次数设置为8次；
2-7)将步骤2-6)得到的结果平均，获得最后的无肌电伪迹的脑电信号，
其中，所述步骤2-5)中使用区间阈值化方法去除肌电伪迹的步骤包括：
2-5-1)阈值计算：
LSMEMD将含有肌电伪迹的信号分解成多个MIMF.每个IMF中包含肌电伪迹的强度各不相同，因此需要对每一个IMF自适应计算肌电伪迹的水平，作为区分肌电伪迹和脑电的阈值。
对于每一个MIMF集合中的任意一个IMF计算其自适应阈值Ti：
其中C是一个调节系数，Ei为第i个IMF能量，T是采样长度；
除第一个IMF外，其他IMF功率谱表现出自相似特征，且IMF的能量呈半对数线性下降趋势。本实施例中，先利用下式(7)所示的估计器计算第一个IMF的能量E1，然后利用式(8)计算其他IMF的能量Ei，u1(t)为第一个IMF：
其中，式(8)中β和ρ为设定的参数；在一种优选的实施例中，β为0.719，ρ为2.01。
2-5-2)基于SCAD区间阈值函数的肌电伪迹去除：
每一个被LSMEMD分解出的MIMF都有可能包含肌电伪迹和脑电信号，为了区分并去除肌电伪迹，本实施例中对MIMF进行了更加细致的处理。将MIMF将中孤立的样本点看做肌电或脑电是没有意义的，且会导致重构信号不连续，如前文所述，脑电和肌电伪迹振荡存在差异，因此将MIMF中两个相邻零点(
和
)组成的过零区间
作为研究对象较为合理，利用过零区间内的极值
与步骤2-5-1)中计算的阈值Ti比较，判断该过零区间内主要是脑电还是肌电伪迹，然后利用相应的阈值函数进行处理。
传统的软阈值函数在处理时会使估计结果偏移一个常数，导致去噪后的信号有失真现象；因为SCAD惩罚估计器同时满足无偏、稀疏和连续等特性，使去噪后的信号具有连续性，同时减小信号的失真。本实施例中采用的SCAD区间阈值函数为：
其中，第i个IMF中第j和j+1个过零点构成的过零区间，
为第i个IMF中区间内所有测量值，
为第i个IMF中的极值，Ti为第i个IMF的阈值，
为去除肌电伪迹后的估计值；其中α为调节因子，在一种优选的实施例中，α为3.7。
本发明还提供一种快速自动去除肌电伪迹的系统，该系统采用如实施例1所述的方法去除脑电中的肌电伪迹。
本发明还提供一种存储介质，其上存储有计算机程序，该程序被执行时用于实现如实施例1所述的方法。
本发明还提供一种计算机设备，包括存储器、处理器以及存储在所述存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现实施例1所述的方法。
实施例2对本发明的方法进行性能检测
1、性能评价指标
为了充分验证方法性能，本实施例选择五种常用的性能评价指标：相对均方根误差(Relative root mean squared error,RRMSE)，平均相关系数(Average correlationcoefficient,ACC)，信噪比(Signal-noise ratio,SNR)，峰值信噪比(Peak SNR,PSNR)和功率谱密度(Power spectral density,PSD)。
1-1、相对均方根误差
相对均方根误差用来评价肌电伪迹去除能力，RRMSE越小，表示去除肌电伪迹后的脑电与真实纯净脑电信号相似度越高，去除肌电伪迹能力越强。
其中，N和T分别代表导联数量和采样长度，XEEG(t)和
分别代表真实纯净脑电信号和去除肌电伪迹后的脑电信号。
1-2、平均相关系数
相关系数是用来反映两个变量之间关系密切程度的指标，在本实施例中用来评估对真实脑电信号的保留性能。计算每个导联真实纯净脑电和去除肌电伪迹后脑电之间的相关系数，用ρ来表示：
1-3、信噪比
信噪比越大，说明信号里的噪声越小，信号的质量越高，否则相反。在实施例中SNR计算公式如下：
1-4、峰值信噪比
峰值信噪比是用来衡量重构信号的质量峰值信噪比较大时，说明真实纯净脑电信号与去除肌电伪迹后的信号有较高的相似度，定义由公式(14)给出：
1-5、功率谱密度
PSD可以很好地表示信号在不同频率的分布情况，通常静息态下脑电多集中于35Hz以下，因此，通过评估去伪迹前后信号的PSD，如果去除肌电伪迹后的信号在35Hz以上频段减少及在1-35Hz频段接近真实信号，可间接评估方法在去除肌电伪迹的同时保留有用脑电信号的能力.
Welch提出了一种改进估计的PSD估计器[9]，它将时间信号分成连续的时间窗，通过将每个时间窗平均得到周期图谱.对于第n个导联xn(t)，在公式(16)给出了功率谱密度的Welch估计，为了简单起见，这里xn(t)简写为x：
其中，xm是x的第m个时间窗，M和K分别代表每个窗口的样本数量和可用窗的数量，
是第m个时间窗的周期图谱；因为在真实数据中无法获取真实纯净的脑电，所以无法将真实数据去除肌电伪迹后的信号与真实纯净的脑电作比较，所以真实数据只能用功率谱密度来评价。
2、数据获取
2-1、仿真数据获取
因为本实施例主要研究脑电中肌电伪迹的去除，因此仿真数据X(t)由真实纯净脑电XEEG(t)和肌电XEMG(t)构成如式(17)所示，不含其它伪迹，使用陷波滤波器去除50Hz工频，低通滤波至32Hz，仿真数据如图3所示。
X(t)＝XEEG(t)+λ·XEMG(t) (17)
XEMG(t)＝ASEMG(t)
其中A是混合矩阵用来控制肌电源SEMG(t)的分布，为了引入足够的空间结构，A的每一列包含5-8个非零项，使每一个独立的肌电源同时存在于5-8个导联中；每一列非零项的数量、位置和具体的值都是根据正态分布随机确定。
在仿真实验中，通过改变(公式21)中参数λ调整SNR，文中SNR值设置0.5-3，步长0.5；
1)导联选择
本实施例研究的重点是少导联脑电场景。导联数量N通常被设置为2-7，但电极的位置是根据不同应用设计的，因为在广泛使用10-20系统中，19导联信号具有更加真实的空间结构。为了模拟真实情况，本实施例在19导联数据中随机选取N导联。
2)真实纯净脑电XEEG(t)
XEEG(t)采集于上海复旦大学附属儿童医院，使用21导联AgCl电极(Fp1，Fp2，F7，F3，Fz，F4，F8，T3，Cz，C4，T4，T5，P3，Pz，P4，T6，O1，O2，A1，A2，其中A1和A2为参考电极)Nicolet One系统采集，按10-20系统安置于头皮.采样频率为500Hz，测量单位uV，使用Nicolet V32放大器采集.为了保证实验的严谨，有经验的临床医生从所有数据中选出50段10s几乎无伪迹的脑电且考虑正常安静时采集到的脑电数据多集中于32Hz以下，对数据进行32Hz低通滤波，所得结果作为真实纯净脑电XEEG(t)。
3)仿真肌电源SEMG(t)
因为肌电伪迹的产生时刻和持续时间是不可预测的，为了充分模拟不同的肌电源，利用随机噪声生成生成连续和瞬时的肌电，采用任意起始时刻和持续时间，振幅也随机给出，并带通滤波至20-100Hz；本实施例生成19导联肌电源SEMG(t)，然后再随机选取N导联肌电源。
2-2、真实数据
本实施例中，真实数据来源于[10]，根据国际10-20系统放置的6通道头皮记录中收集，电极为Fp1，Fp2，C3，C4，O1，O2，参考电极置于Cz位置。采样率为256Hz，持续时长30s，大部分的肌电是由抬眉毛产生的，图4即为整段真实数据，图中可以明显看到在10s、17s和25s左右可以观察到肌电伪迹，在Fp1和Fp2最明显，O1和O2电极因为离眉毛最远所以受影响最轻；同时，在数据的采集中，被试者属于静息状态，几乎不含有γ波，所以纯净脑电主要分布在35Hz甚至更低的频率以下，由于本实施例提出的方法是计划实时使用的，截取3段持续时间为8秒的受肌电伪迹影响的脑电源，起始时间分别是6s、14s和21s。
3、实验结果与分析
为了全面评价本实施例算法的性能，进行了两类实验：首先，由于在现实中无法获得未受伪迹影响的纯净脑电，使用仿真数据和RRMSE、SNR、PSNR、ACC、PSD评价指标对本实施例算法进行验证；其次，使用采集到的被肌电伪迹影响真实数据和PSD验证本实施例方法的真实性能。
本实施例选择CCA、LSMEMD-PR、LSMEMD-CCA为对比方法，其原因如下：如前文所述，整体上肌电伪迹去除方法包括自适应滤波、BSS和信号分解等方法，其中自适应滤波方法需要大量参考电极，本实施例数据无法满足，不作为本实施例研究的对比方法；BSS类方法中，考虑CCA被证明去除肌电伪迹性能优于ICA等其他BSS方法，所以选择CCA作为本实施例对比方法之一；信号分解类方法，因WD类方法选择小波基等因素对实验结果影响较大，不列入本实施例对比方法；对于EMD类方法，MEMD类方法是目前被证明优于EEMD和和EMD等方法的方法，所以选择基于LSMEMD的方法进行比较，主要包括LSMEMD-PR和LSMEMD-CCA，即分别为基于LSMEMD分解的IMF部分重构方法和IMF源分离方法。
3-1、仿真结果和分析
3-1-1、仿真实验及结果
首先，为了验证所提出方法去除肌电伪迹的性能，将实施例1的快速自动去除肌电伪迹的方法(以下简称为：LSMEMD-IT)与CCA、LSMEMD-PR、LSMEMD-CCA 3种方法进行比较，正如前文中提到的，将生成的19导联被肌电影响的脑电随机生成6组2-7导联仿真数据，每组仿真数据均采用不同信噪比进行测试，每组实验重复50次。仿真性能评价指标RRMSE、PSNR、去除肌电伪迹后SNR和ACC如图5a-5d所示，CCA、LSMEMD-PR、LSMEMD-CCA和LSMEMD-IT对应的结果分别用不同曲线表示。
其次，因为LSMEMD-PR、LSMEMD-CCA和LSMEMD-IT三种方法性能接近，同时三种方法都是对IMF进行操作，只是处理方式不同，因此设计实验，进一步研究三种方法去除肌电伪迹后的IMF对重构信号的影响；通常采集脑电时，肌电伪迹较大，实验选择SNR＝0.5，导联选择3个和6个，实验结果如图7a-7d和8a-8d所示。
最后，关于算法参数的设置，在CCA方法中，归一化参数设置为0.001；因为真实纯净脑电低通滤波至32Hz.考虑到滤波器的衰减，考虑脑电典型的节律分布包括δ(1－3Hz)、θ(4－7Hz)、α(8－13Hz)、β(14－32Hz)等，本实施例中将信号分解为4个imf和1个余差，其中imf的频带分布与上述脑电典型节律分布一致。同理，仿真数据低通滤波至70Hz，将仿真数据分解为6个IMF，不失一般性，辅助噪声通道设置为5，对于投影方向数量，在优选的实施例中，投影数量通常大于6×N(N为导联数量)，本实施例设置为64。对于CCA、LSMEMD-PR和LSMEMD-CCA方法阈值的选择，在仿真实验中，因为真实脑电是已知的，为了公平起见，本实施例为这3种方法选择最优阈值，最优阈值是指通过移除某些成分得到RRMSE最小，本实施例中，LSMEMD-IT方法中通用阈值常数C设置为0.8。
3-1-2、LSMEMD-IT对比CCA分析
在图5中可以明显看出，CCA方法相比本实施例方法有着较大差距，CCA方法得到了最大的RRMSE和最小的ACC。这一结果表明，CCA方法不能有效地去除肌电伪迹同时保留脑电信号，更具体地说，如果脑电信号导联的数量远小于隐藏源的数量，那么数量非常有限的源就会被CCA解混；由于这些源很可能来自大脑和肌肉活动，所以对任何源的直接去除都有可能导致大脑活动的潜在损失，从而导致大的RRMSE和低的ACC，所以CCA方法在少导联时并不能很好地工作。
3-1-3、LSMEMD-IT对比LSMEMD-PR和LSMEMD-CCA分析
这三种方法不同点在于对IMF的处理方式：LSMEMD-PR将原始信号自适应地分解成多个MIMF，计算每个IMF的自相关性，肌电通常具有较低的自相关性，因此将自相关性较低的IMF整个去除，其余IMF重构生成处理后的信号；LSMEMD-CCA方法使用CCA对IMF进行处理，计算隐藏源，自相关性较低的源被视为肌电伪迹并去除，利用逆CCA得到干净的IMF，最后将IMF重构生成处理后的信号；LSMEMD-IT方法依据IMF能量和特性自动计算阈值，利用阈值和阈值函数对每个IMF中相应的过零区间进行阈值操作，从而去除肌电伪迹得到无肌电伪迹IMF，重构生成无肌电伪迹信号。
通过图5可以发现，LSMEMD-PR、LSMEMD-CCA、LSMEMD-IT三种方法在导联数量较少时性能差距相对较大。合理的解释是，对于LSMEMD-PR方法，因为肌电伪迹和脑电并不会严格的集中在某个IMF上，直接将整个IMF去除，可能会导致去除脑电信号，使重构信号出现偏差，当导联数量越少时，可用信息越少，对任何信息的错误处理，导致重构信号的偏差会加大；对于LSMEMD-CCA方法，虽然其利用CCA计算隐藏源，但当导联数量过少时，该方法计算出的源的数量也会随之减少，对少量的源利用阈值去除，依然会导致误去除脑电问题，如图6所示，在两个导联，SNR＝0.5时，利用LSMEMD-CCA方法分解仿真数据得到的隐藏源(图6中的original sources)和去除肌电伪迹后的隐藏源(图6中的calculatied sources)，分别用标号1(波浪线)和2(近似直线)表示；在图中可以看到，第6个源明显包含脑电成分，然而LSMEMD-CCA方法却将其判定为肌电伪迹并大量去除，使重构信号失真，从而导致具有较大的RRMSE和较小的ACC。综上所述，LSMEMD-PR和LSMEMD-CCA方法在导联数量较少时，不能充分发挥作用。
上述三种方法主要不同点在于IMF的处理方式，下面本实施例将用两组实验，来说明不同方法去除肌电伪迹后IMF对重构信号的影响。为了直观的表示三种方法的肌电伪迹去除性能，从仿真数据中随机选取3导联和6导联以及SNR＝0.5时进行去除肌电伪迹实验，实验结果在图7(a)和8(a)中展示，其中，真实脑电(图7a中的Ground truth EEG，标号1)和含有肌电伪迹的仿真数据(图7a中的Contaminated EEG，标号2)以及LSMEMD-PR(标号3)、LSMEMD-CCA(标号4)和LSMEMD-IT(标号5)方法去除肌电伪迹结果分别如图所示，需要理解的是，其中线条相互重叠，故再通过标号进行区分；可以看到三种方法的性能总体差别不大，但通过随机截取的部分信号(1500-2000采样点)进行放大展示可以看到，LSMEMD-PR和LSMEMD-CCA方法在3导联时过度去除脑电，出现了失真现象，而本实施例方法在3导联和6导联时都能够充分抑制肌电伪迹同时保留脑电信号。
下面详细讨论IMF对重构结果的影响：首先，三种方法去除肌电伪迹后IMF的差异如图7(b)和8(b)所示，本实施例随机挑选了第2导联进行示意，真实脑电(图7b中的Groundtruth EEG)、LSMEMD-IT、LSMEMD-CCA和LSMEMD-PR结果分别用标号1、2、3、4表示(需要理解的是，其中线条相互重叠，故再通过标号进行区分)，其中真实脑电分解出的IMF用G-IMF表示，因为随IMF阶次增大，频率降低，所以频率较大的肌电通常在阶次较小的IMF中；图中可以看到，三种方法对IMF1、IMF2和IMF3去除肌电伪迹比较明显，说明肌电伪迹主要在前3个IMF中.LSMEMD-PR方法将前3个IMF认定为肌电伪迹，全部去除；LSMEMD-CCA和LSMEMD-IT方法在去除肌电伪迹后IMF1和IMF2保留了极少的信号，能量非常小，几乎可以忽略不计，可以认为IMF1和IMF2是纯肌电伪迹，但在IMF3上都有明显的信号，并且存在较大差异。
为了进一步研究上述差异对去除肌电伪迹性能的影响，利用PSD和中心频率对三种方法去除肌电伪迹后IMF和真实脑电IMF进行对比，结果如图7(c)和8(c)所示；真实脑电(图7c中的Ground truth EEG)、LSMEMD-IT、LSMEMD-CCA和LSMEMD-PR结果分别用标号1、2、3和4表示，中心频率用5表示(需要理解的是，其中线条相互重叠，故再通过标号进行区分)。由上文所述，只对三种方法处理后的IMF3…IMF5…R进行展示，从图中可以看到真实脑电G-IMF1的中心频率在25Hz左右，这时肌电伪迹与脑电频率重叠最显著；LSMEMD-PR方法直接将仿真数据分解出的IMF3去除，IMF3同时包含肌电伪迹和脑电，直接去除会导致丢失大量脑电；而LSMEMD-CCA和LSMEMD-IT方法都对IMF3做了更加细致的处理，但LSMEMD-CCA方法处理之后，IMF3与G-IMF1差异非常大，出现了明显的下降趋势，并且中心频率偏移至10Hz左右，去除了大量脑电信号；不仅如此，在IMF4中15Hz以上也出现了脑电误去除问题，而LSMEMD-IT方法对IMF去除肌电伪迹后，中心频率几乎没有改变，PSD曲线与真实脑电更接近，说明LSMEMD-IT方法在肌电与脑电频谱重合较为严重时，依然能够去除肌电伪迹同时保留脑电。
最后，不失一般性，在图7(d)和8(d)给出了所有导联三种方法去除伪迹前后的PSD；图7(d)中真实脑电(Ground truth EEG)用标号1标记，含有肌电伪迹的仿真数据(Contaminated EEG)用标号2标记。需要理解的是，图中线条相互重叠，故再通过标号进行区分；图8(a)中，含有肌电伪迹的仿真数据(Contaminated EEG，标号1)、真实脑电(Groundtruth EEG，标号2)和以及LSMEMD-PR(标号3)、LSMEMD-CCA(标号4)和LSMEMD-IT(标号5)如图所示；图8(b)中，真实脑电(Ground truth EEG，标号1)LSMEMD-IT(标号2)、LSMEMD-CCA(标号3)和LSMEMD-PR(标号4)；图8(c)中，真实脑电(Ground truth EEG，标号1)、LSMEMD-IT(标号2)、LSMEMD-CCA(标号3)和LSMEMD-PR(标号4)；图8(d)中，真实脑电(Ground truthEEG)用标号1标记，含有肌电伪迹的仿真数据(Contaminated EEG)用标号2标记。
图中清楚地看出三种方法在15Hz以下时，脑电保留性能良好，处理后的曲线与真实脑电曲线几乎重合，但超过15Hz以后LSMEMD-PR和LSMEMD-CCA方法开始呈现下降趋势，并且下降较大，将脑电误判为肌电伪迹去除。因为LSMEMD-PR方法将IMF3直接去除，去除的脑电最多，所以曲线最低；而LSMEMD-IT方法在超过15Hz以后与其他方法相比整体趋势更接近真实脑电，说明LSMEMD-IT方法能够最大限度去除肌电伪迹同时保留脑电信号。
综上所述，在仿真数据中与对比方法相比，本发明提供的LSMEMD-IT方法在少导联场景下能够充分抑制肌电伪迹并保留真实脑电信号；当脑电与肌电伪迹频谱严重重合时，这种优势更加明显。
3-2、真实结果和讨论
3-2-1、真实实验及结果
为了充分说明本实施例方法在真实数据中的性能，选择3段具有明显肌电伪迹的数据，如图10(b)、11(b)和12(b)所示，6个导联，持续时长均为8s；与仿真实验参数相同，IMF数量为6个和1个余差，通过真实实验得到三个数据源不同方法的PSD结果图10(a)、11(a)和12(a)。需要理解的是，图中线条相互重叠，故再通过标号进行区分：图10(a)、11(a)和12(a)中真实数据(real)、CCA、LSMEMD-PR、LSMEMD-CCA-0.8、LSMEMD-CCA-0.95分别用标号1、2、3、4、5、6标注(需要理解的是，图中线条相互重叠，故再通过标号进行区分)；图10(b)、图11(b)中、图12(b)，真实数据(real)、LSMEMD-IT分别用标号1、2标注。
由于真实脑电不可用，所以CCA、LSMEMD-PR和LSMEMD-CCA方法无法使用最优阈值，必须在使用去除方法之前确定阈值，肌电伪迹与脑电信号相比，有较低的自相关性，因此需要人为设置阈值来移除自相关性低于阈值的肌电伪迹。
表1给出LSMEMD-PR方法分解出IMF的自相关系数，可以看到IMF2和IMF3相关系数具有相对较大的差值且IMF2数值较低，因此在本实施例中将每个导联分解出的前两个IMF设置为与肌电伪迹相关的IMF并去除。CCA和LSMEMD-CCA方法源的自相关系数在图9中给出，CCA方法中，最多可以产生6个隐藏源，因此CCA方法自相关系数相对较低，设置0.8为CCA方法的阈值，对于LSMEMD-CCA方法，在第24和30个源的位置均有较为明显的下降，认为这些源以下是肌电成分.表2给出在下降处附近源的自相关系数，依据表2给出的结果，设置0.8和0.95为LSMEMD-CCA方法的阈值；LSMEMD-IT方法中设置常数C为2。
表1真实数据三个源分解后IMF的自相关系数
表2 LSMEMD-CCA方法中自相关系数下降处附近的值
3-2-2、利用PSD对比讨论真实结果
对于CCA方法在真实数据下的去除肌电伪迹性能，通过图10(a)、11(a)和12(a)可以发现，CCA方法去除肌电伪迹后的PSD在所有频段都有下降，但在肌电伪迹较多的高频段下降幅度并不大。正如之前所述，真实数据只有6个导联，对于CCA方法，计算源的数量也在6个以内，少量的源就表示每个源中可能同时存在肌电和脑电，对任何一个源的去除就等于同时去除了肌电伪迹和脑电信号。总之，对于少导联脑电，CCA不能最大限度的抑制肌电伪迹同时保留脑电。
对于LSMEMD-PR和LSMEMD-CCA方法：当LSMEMD-CCA阈值等于0.95时，在15Hz甚至更低频率时便开始下降，说明去除了部分脑电信号，阈值选取过大；当LSMEMD-CCA阈值等于0.8时，虽然在低频段保留了脑电信号，但在高频段(50Hz以上)对纯肌电的抑制不足，说明阈值选取过小；而对于LSMEMD-PR方法，它的性能和LSMEMD-CCA阈值等于0.8时几乎一样；从自相关系数可以看到，LSMEMD-CCA方法在系数较大时非常密集，虽然阈值只相差0.15，但包含6个源，且去除肌电伪迹能力差异较大.这也说明LSMEMD-CCA方法不仅自适应性较差，并且对阈值的准确性要求较高。
对于LSMEMD-IT方法：在低频段具有很好的脑电保留能力，在20Hz之前与真实数据PSD几乎重合；在20Hz之后曲线开始下降，并在50Hz以上达到较低的PSD值，充分抑制高频肌电伪迹.不仅如此，本实施例方法在中间频段依然有着优于对比方法的性能，虽然LSMEMD-CCA阈值等于0.8时和LSMEMD-IT方法都在20Hz左右开始下降，但可以看到本实施例方法下降速度快于LSMEMD-CCA方法，在同一频率下，有更低的PSD.正如仿真实验中所述，中间频段肌电和脑电频率重叠最严重，本实施例方法在中间频段使肌电伪迹得到充分抑制，在充分去除肌电伪迹的同时最大限度保留脑电信号。
最后，伪迹去除结果如图10(b)、11(b)和12(b)所示，LSMEMD-IT充分去除了肌电伪迹，保留了脑电信号。综上所述，本实施例的方法在真实数据中与对比方法相比，自适应性更强，使其成为一种可靠的少导联脑电肌电伪迹去除方法。
4、结论
本实施例提出的一种自适应、有效的少导联脑电肌电伪迹去除方法LSMEMD-IT，该方法基于LSMEMD和区间阈值进行改进，解决了目前方法自适应性较差，导联数量过少去除效果不佳以及脑电与肌电伪迹频谱严重重叠时脑电误去除问题；通过仿真和真实数据实验，与传统的BSS方法CCA，信号分解部分重构方法LSMEMD-PR和目前较有效的少导联去除肌电伪迹方法LSMEMD-CCA进行比较，进一步证明本发明提供的LSMEMD-IT方法优于现有方法的性能。
本发明基于LSMEMD和IT提出一种自适应去除肌电伪迹的方法：LSMEMD-IT，利用LSMEMD对少导联脑电进行自适应时频同步分解，符合导联间的同步和一致性；SCAD区间阈值化是一种符合IMF特性的自适应阈值计算方法，对LSMEMD分解产物MIMF进行更加细致的处理，避免了源分离方法中的计算误差，同时减小了部分重构方法误去除脑电导致的信号失真问题。本发明使用仿真和真实数据验证LSMEMD-IT方法，与CCA、最小二乘多维经验模态分解部分重构(LSMEMD-partial reconstruction,LSMEMD-PR)和LSMEMD-CCA进行比较，证明了LSMEMD-IT方法的有效性。
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节。
Claims (10)
1.一种快速自动去除肌电伪迹的方法，其特征在于，包括以下步骤：
S1、MIMF求解：将N通道含有肌电伪迹的真实脑电和M通道的辅助噪声构成(N+M)维复合信号，使用最小二乘多维经验模态分解方法将该(N+M)维复合信号分解为(N+M)个MIMF集合，去除辅助噪声分解得到的MIMF，只保留脑电信号分解出的N个MIMF集合；
S2、肌电伪迹去除：针对基于所述步骤S1所得的N个MIMF集合，使用SCAD区间阈值化方法，根据IMF能量和特性自适应计算N个MIMF集合中的每个IMF的阈值，联合阈值函数对MIMF中每个过零区间进行肌电伪迹阈值去除，得到无肌电伪迹的MIMF；其中，IMF表示本征模态函数，MIMF表示多维本征模态函数，即1个MIMF集合中包括多个IMF；
S3、信号重构：对每一个通道的无肌电伪迹MIMF进行重构，得到N维不含肌电伪迹的脑电。
2.根据权利要求1所述的快速自动去除肌电伪迹的方法，其特征在于，所述步骤S1中采用最小二乘多维经验模态分解方法对脑电信号进行自适应的时频分解，其具体步骤包括：
1-1)假设X(t)是一个N通道脑电信号，t＝1,…,T，Y(t)是一个M通道辅助噪声信号，M是任意的正整数，输入信号Z(t)＝[X(t),Y(t)]；
1-2)设i＝1，并产生L组用于均匀投影的超球面向量点集；
1-3)计算输入信号Z(t)的第l个投影函数fl(t)，计算方式为：
fl(t)＝Z(t)·(Vl(t))T (2)；
1-4)对于所有l，利用EMD算法提取相应投影函数fl(t)的第一层单变量IMFul(t)；
1-5)联合所有IMF函数
以及它们的投影向量
则相应的(N+M)维MIMF函数Di(t)可以通过解下列超定线性方程组获得
1-6)计算残差函数：
S(t)＝Z(t)-Ui(t) (4)；
1-7)如果S(t)没有足够的极值组成包络，停止当前算法并获取最终残差项R(t)＝S(t)；反之则更新当前输入为Z(t)＝S(t)及i＝i+1，然后转至步骤1-3)；
1-8)在停止上述筛选过程后，获得Q组(N+M)维MIMF和残差R(t)，信号Z(t)表示为：
3.根据权利要求2所述的快速自动去除肌电伪迹的方法，其特征在于，所述步骤S2中，采用基于SCAD区间阈值化的方法进行MIMF肌电去除，具体包括以下步骤：
2-1)对于每一个MIMF集合，其中包括Q个IMF，对其中的第2～Q个IMF进行重组，得到
2-2)随机改变第1个IMF采样点的位置，得到
2-3)构建一个不同于原始信号的新的含肌电伪迹信号，得到
2-4)对新信号Xa(t)做最小二乘多维经验模态分解，得到
2-5)对步骤2-4)分解出的MIMF使用区间阈值化方法去除肌电伪迹，重构得到信号
2-6)将步骤2-2)至2-5)迭代K-1次，得到K个去除肌电伪迹信号：
2-7)将步骤2-6)得到的结果平均，获得最后的无肌电伪迹的脑电信号，
4.根据权利要求3所述的快速自动去除肌电伪迹的方法，其特征在于，所述步骤2-5)中使用区间阈值化方法去除肌电伪迹的步骤包括：
2-5-1)阈值计算：
对于每一个MIMF集合中的任意一个IMF计算其自适应阈值Ti：
其中C是一个调节系数，Ei为第i个IMF能量，T是采样长度；
先利用下式(7)所示的估计器计算第一个IMF的能量E1，然后利用式(8)计算其他IMF的能量Ei，u1(t)为第一个IMF：
其中，式(8)中β和ρ为设定的参数；
2-5-2)基于SCAD区间阈值函数的肌电伪迹去除：
SCAD区间阈值函数为：
其中，第i个IMF中第j和j+1个过零点构成的过零区间，
为第i个IMF中区间内所有测量值，
为第i个IMF中的极值，Ti为第i个IMF的阈值，
为去除肌电伪迹后的估计值；其中α为调节因子。
5.根据权利要求4所述的快速自动去除肌电伪迹的方法，其特征在于，所述步骤2-6)中，K＝9。
6.根据权利要求5所述的快速自动去除肌电伪迹的方法，其特征在于，步骤2-5-1)中，β为0.719，ρ为2.01。
7.根据权利要求6所述的快速自动去除肌电伪迹的方法，其特征在于，步骤2-5-2)中，α为3.7。
8.一种快速自动去除肌电伪迹的系统，其特征在于，该系统采用如权利要求1-7中任意一项所述的方法去除脑电中的肌电伪迹。
9.一种存储介质，其上存储有计算机程序，其特征在于，该程序被执行时用于实现如权利要求1-7中任意一项所述的方法。
10.一种计算机设备，包括存储器、处理器以及存储在所述存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-7中任意一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011164260.4A CN112294340A (zh) | 2020-10-27 | 2020-10-27 | 快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011164260.4A CN112294340A (zh) | 2020-10-27 | 2020-10-27 | 快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112294340A true CN112294340A (zh) | 2021-02-02 |
Family
ID=74331042
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011164260.4A Pending CN112294340A (zh) | 2020-10-27 | 2020-10-27 | 快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112294340A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112957055A (zh) * | 2021-02-05 | 2021-06-15 | 中国科学院深圳先进技术研究院 | 基于eemd-pca去除eeg信号中运动伪迹的方法及装置 |
| CN112975982A (zh) * | 2021-03-16 | 2021-06-18 | 北京理工大学 | 一种基于脑机融合的空地协同多机器人系统 |
| CN114288520A (zh) * | 2021-12-31 | 2022-04-08 | 广州酷狗计算机科技有限公司 | 一种基于脑电波的辅助睡眠方法、装置、设备及存储介质 |
| CN114886388A (zh) * | 2022-07-12 | 2022-08-12 | 浙江普可医疗科技有限公司 | 一种麻醉深度监测过程中脑电信号质量的评估方法及装置 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101869477A (zh) * | 2010-05-14 | 2010-10-27 | 北京工业大学 | 一种自适应脑电信号中眼电伪迹的自动去除方法 |
| CN105342605A (zh) * | 2015-12-09 | 2016-02-24 | 西安交通大学 | 一种去除脑电信号中肌电伪迹的方法 |
| CN106805945A (zh) * | 2017-01-22 | 2017-06-09 | 合肥工业大学 | 一种少数通道的脑电信号中肌电伪迹的消除方法 |
| CN108309290A (zh) * | 2018-02-24 | 2018-07-24 | 华南理工大学 | 单通道脑电信号中肌电伪迹的自动去除方法 |
-
2020
- 2020-10-27 CN CN202011164260.4A patent/CN112294340A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101869477A (zh) * | 2010-05-14 | 2010-10-27 | 北京工业大学 | 一种自适应脑电信号中眼电伪迹的自动去除方法 |
| CN105342605A (zh) * | 2015-12-09 | 2016-02-24 | 西安交通大学 | 一种去除脑电信号中肌电伪迹的方法 |
| CN106805945A (zh) * | 2017-01-22 | 2017-06-09 | 合肥工业大学 | 一种少数通道的脑电信号中肌电伪迹的消除方法 |
| CN108309290A (zh) * | 2018-02-24 | 2018-07-24 | 华南理工大学 | 单通道脑电信号中肌电伪迹的自动去除方法 |
Non-Patent Citations (3)
| Title |
|---|
| YAN LIU 等: "An Efficient and Robust Muscle Artifact Removal Method for Few-Channel EEG", 《IEEE ACCESS 2019》 * |
| YANNIS KOPSINIS 等: "Development of EMD-Based Denoising Methods Inspired by Wavelet Thresholding", 《IEEE TRANSACTIONS ON SIGNAL PROCESSING》 * |
| YANNIS KOPSINIS 等: "EMPIRICAL MODE DECOMPOSITION BASED SOFT-THRESHOLDING", 《2008 16TH EUROPEAN SIGNAL PROCESSING CONFERENCE》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112957055A (zh) * | 2021-02-05 | 2021-06-15 | 中国科学院深圳先进技术研究院 | 基于eemd-pca去除eeg信号中运动伪迹的方法及装置 |
| CN112957055B (zh) * | 2021-02-05 | 2023-12-29 | 中国科学院深圳先进技术研究院 | 基于eemd-pca去除eeg信号中运动伪迹的方法及装置 |
| CN112975982A (zh) * | 2021-03-16 | 2021-06-18 | 北京理工大学 | 一种基于脑机融合的空地协同多机器人系统 |
| CN112975982B (zh) * | 2021-03-16 | 2021-11-09 | 北京理工大学 | 一种基于脑机融合的空地协同多机器人系统 |
| CN114288520A (zh) * | 2021-12-31 | 2022-04-08 | 广州酷狗计算机科技有限公司 | 一种基于脑电波的辅助睡眠方法、装置、设备及存储介质 |
| CN114886388A (zh) * | 2022-07-12 | 2022-08-12 | 浙江普可医疗科技有限公司 | 一种麻醉深度监测过程中脑电信号质量的评估方法及装置 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chen et al. | A novel EEMD-CCA approach to removing muscle artifacts for pervasive EEG | |
| CN112294340A (zh) | 快速自动去除肌电伪迹的方法、系统、存储介质及计算机设备 | |
| US7809433B2 (en) | Method and system for limiting interference in electroencephalographic signals | |
| CN110338786B (zh) | 一种癫痫样放电的识别与分类方法、系统、装置和介质 | |
| CN108577834B (zh) | 一种用于癫痫间期棘波自动检测的方法 | |
| EP1788937A2 (en) | Method for adaptive complex wavelet based filtering of eeg signals | |
| Hossain et al. | A robust ECG denoising technique using variable frequency complex demodulation | |
| Dai et al. | Removal of ECG artifacts from EEG using an effective recursive least square notch filter | |
| CN105342605A (zh) | 一种去除脑电信号中肌电伪迹的方法 | |
| Liu et al. | An efficient and robust muscle artifact removal method for few-channel EEG | |
| CN106236080A (zh) | 基于多通道的脑电信号中肌电噪声的消除方法 | |
| CN109480832A (zh) | 一种单通道的脑电信号中肌电伪迹的消除方法 | |
| Al-kadi et al. | Compatibility of mother wavelet functions with the electroencephalographic signal | |
| Uthayakumar et al. | Multifractal-wavelet based denoising in the classification of healthy and epileptic EEG signals | |
| CN110292374B (zh) | 基于奇异谱分析和变分模态分解的心电信号去基线漂移方法 | |
| Ren et al. | Noise reduction based on ICA decomposition and wavelet transform for the extraction of motor unit action potentials | |
| CN115054269A (zh) | 一种基于apso-vmd算法的ecg信号去噪方法及系统 | |
| CN117349598B (zh) | 脑电信号处理方法及装置、设备、存储介质 | |
| Elouaham et al. | Combination time-frequency and empirical wavelet transform methods for removal of composite noise in EMG signals | |
| El Khadiri et al. | A Comparison of the Denoising Performance Using Capon Time-Frequency and Empirical Wavelet Transform Applied on Biomedical Signal. | |
| Shen et al. | Method for extracting time-varying rhythms of electroencephalography via wavelet packet analysis | |
| Huang et al. | A novel application of the S-transform in removing powerline interference from biomedical signals | |
| Kaushal et al. | Better approach for denoising EEG signals | |
| Hu et al. | ECG cancellation for surface electromyography measurement using independent component analysis | |
| Malan et al. | Removal of Ocular Atrifacts from Single Channel EEG Signal Using DTCWT with Quantum Inspired Adaptive Threshold |
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: 20210202 |
|
| RJ01 | Rejection of invention patent application after publication |