CN116584960A - 一种膈肌肌电信号降噪方法 - Google Patents
一种膈肌肌电信号降噪方法 Download PDFInfo
- Publication number
- CN116584960A CN116584960A CN202310353452.7A CN202310353452A CN116584960A CN 116584960 A CN116584960 A CN 116584960A CN 202310353452 A CN202310353452 A CN 202310353452A CN 116584960 A CN116584960 A CN 116584960A
- Authority
- CN
- China
- Prior art keywords
- wavelet
- noise reduction
- noise
- signals
- sequence
- 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
- 230000009467 reduction Effects 0.000 title claims abstract description 99
- 238000000034 method Methods 0.000 title claims abstract description 62
- 238000012880 independent component analysis Methods 0.000 claims abstract description 52
- 230000009466 transformation Effects 0.000 claims abstract description 36
- 238000012545 processing Methods 0.000 claims abstract description 34
- 230000003183 myoelectrical effect Effects 0.000 claims abstract description 33
- 239000011159 matrix material Substances 0.000 claims description 19
- 238000004364 calculation method Methods 0.000 claims description 11
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 238000004590 computer program Methods 0.000 claims description 10
- 108010076504 Protein Sorting Signals Proteins 0.000 claims description 6
- 238000002156 mixing Methods 0.000 claims description 5
- 210000003205 muscle Anatomy 0.000 claims description 5
- 238000005516 engineering process Methods 0.000 abstract description 2
- 230000009977 dual effect Effects 0.000 abstract 1
- 230000000694 effects Effects 0.000 description 21
- 238000010586 diagram Methods 0.000 description 6
- 238000001228 spectrum Methods 0.000 description 6
- 230000004048 modification Effects 0.000 description 5
- 238000012986 modification Methods 0.000 description 5
- 238000011156 evaluation Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000011946 reduction process Methods 0.000 description 2
- 230000000241 respiratory effect Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 208000006545 Chronic Obstructive Pulmonary Disease Diseases 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000004556 brain Anatomy 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000001035 drying Methods 0.000 description 1
- 238000004070 electrodeposition Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 210000003238 esophagus Anatomy 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 208000001797 obstructive sleep apnea Diseases 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 210000003019 respiratory muscle Anatomy 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
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/389—Electromyography [EMG]
-
- 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/25—Bioelectric electrodes therefor
- A61B5/279—Bioelectric electrodes therefor specially adapted for particular uses
- A61B5/28—Bioelectric electrodes therefor specially adapted for particular uses for electrocardiography [ECG]
-
- 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/25—Bioelectric electrodes therefor
- A61B5/279—Bioelectric electrodes therefor specially adapted for particular uses
- A61B5/296—Bioelectric electrodes therefor specially adapted for particular uses for electromyography [EMG]
-
- 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/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
-
- 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/389—Electromyography [EMG]
- A61B5/397—Analysis of electromyograms
-
- 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/7253—Details of waveform analysis characterised by using transforms
- A61B5/726—Details of waveform analysis characterised by using transforms using Wavelet transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/10—Pre-processing; Data cleansing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Veterinary Medicine (AREA)
- Pathology (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Biophysics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Artificial Intelligence (AREA)
- Signal Processing (AREA)
- Cardiology (AREA)
- Physiology (AREA)
- Psychiatry (AREA)
- Theoretical Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明涉及一种膈肌肌电信号降噪方法，包括以下步骤，S1，借助于采集装置获取膈肌肌电信号和参考心电噪声信号；S2，对膈肌肌电信号和参考心电噪声信号进行小波变换，得到第一小波系数序列；S3，对第一小波变换系数序列利用独立成分分析法进行第一降噪处理，得到第二小波系数序列；S4，对第二小波系数序列利用小波熵法进行第二降噪处理，得到第三小波系数序列；S5，对第三小波系数序列进行小波逆变换，得到降噪的膈肌肌电信号。本发明采用独立成分分析一次降噪叠加小波熵二次降噪的二重降噪技术，能够有效抑制和降低心电噪声对膈肌肌电信号采集的影响，提升膈肌肌电信号采集的有效性和准确性。
Description
技术领域
本发明涉及一种膈肌肌电信号降噪方法。
背景技术
膈肌是人体主要的呼吸肌，占呼吸过程的80％，膈肌肌电信号是由膈肌运动而产生的一种非平稳的生物电信号，膈肌肌电可以用来评估慢性阻塞性肺疾病(COPD)患者的膈肌疲劳状态，与其他呼吸道信号结合可以用来计算阻塞性睡眠呼吸暂停患者的气道阻力等生理特征，以作为临床的诊断和评估提供科学的依据。
膈肌肌电的采集方式有三种，分别是：针电极、食道电极和体表电极，其中针电极和食道电极的采集方式噪声干扰小，有创，病人耐受度低，在临床推广度有限。而体表电极虽噪声干扰较大但因为其无创、采集方便的特点而在临床上具有广泛的应用前景。
在利用体表电极采集的膈肌肌电信号中，心电噪声的干扰是不可忽略的。心电噪声的幅度较膈肌肌电高出数倍，且心电噪声的频谱和膈肌肌电信号的频谱有较大的重叠。而膈肌肌电的频谱特征也是临床评估膈肌功能的重要指标，因此心电噪声是必须要被去除的一种干扰。
现有技术中，大多采用独立成分分析ICA(Independent Component Analysis)的方法进行多通道膈肌肌电信号的降噪处理，但是由于膈肌肌电信号中心电噪声幅度较强且与膈肌肌电频谱存在混叠，并且在实际应用中，采集数据的各个通道情况不同、电极位置距离心脏远近等的差异，导致各通道收集的心电噪声在幅度和波形上出现变化，单独使用独立成分分析法对心电噪声进行去除会出现多个通道之间降噪效果不均匀且整体降噪效果不佳的问题。
针对这样的情况，本发明在采用独立成分分析法进行一次降噪的基础上，利用小波熵的方法对一次降噪后效果不佳的通道进行二次降噪，能够有效解决独立成分分析法对心电噪声在各通道去除效果不统一、部分通道降噪效果不佳的问题。
同时，为了增强独立成分分析方法处理过程中对噪声的分离效果，本发明在独立成分分析方法的输入中加入参考心电噪声信号，从而在独立成分分析中能够更有效地提取噪声信息，增强降噪效果。其中，参考心电噪声信号的采集可利用采集装置在人体体表左侧靠近心脏位置的体表电极实现。
发明内容
(一)要解决的技术问题
鉴于现有技术的上述缺点、不足，本发明提供了一种膈肌肌电信号降噪方法。
(二)技术方案
为了达到上述目的，本发明采用的主要技术方案包括：
第一方面，本发明实施例提供一种膈肌肌电信号降噪方法，包括以下步骤：
S1，借助于采集装置获取膈肌肌电信号和参考心电噪声信号；
S2，对所述膈肌肌电信号和所述参考心电噪声信号进行小波变换，得到第一小波系数序列；
S3，对所述第一小波变换系数序列利用独立成分分析法进行第一降噪处理，得到第二小波系数序列；
S4，对所述第二小波系数序列利用小波熵法进行第二降噪处理，得到第三小波系数序列；
S5，对所述第三小波系数序列进行小波逆变换，得到降噪的膈肌肌电信号。
可选地，所述S2包括：
对所述膈肌肌电信号和所述参考心电噪声信号分别进行5层Harr小波变换，将变换后的各信号的小波系数合并为所述第一小波系数序列。
可选地，所述S3包括：
S31，将所述第一小波系数序列Wx＝[Wx1；Wx2；…；Wxn]作为独立成分分析算法的输入，经过独立成分分析得到源信号矩阵Y＝[y1；y2；…；yn]，以及解混矩阵W和混合矩阵A，
其中，Wx和Y为n×N的矩阵，W与A都是n×n的矩阵，N为小波系数向量长度，n为所有信号数量，包括若干个通道的膈肌肌电信号和参考心电噪声信号；
S32，将得到的各个源信号序列yi与参考心电噪声对应的小波变换系数序列Wx噪声做相关运算得到相关系数的绝对值θi，计算公式如下：
S33，根据预设相关度判定阈值T1，按照以下公式处理各个源信号序列yi：
S34，将处理过的yi放入对应Y中的第i行，得到处理后的源数据矩阵Y1；
S35，将Y1用混合矩阵A进行独立成分分析反变换，得到第一降噪后的各通道信号的第二小波变换系数序列Wx1，计算公式如下：
Wx1＝A*Y1。
可选地，对所述第一小波变换系数序列利用FASTICA算法执行独立成分分析和第一降噪处理。
可选地，所述S4包括：
S41，对第二小波变换系数序列根据预设噪声指标，选取超过预设噪声指标的通道组成第二降噪处理序列Wx2，其他通道组成Wx11；
S42，将第二降噪处理序列Wx2利用小波熵法进行第二降噪，得到Wx22；
S43，将第二降噪后的Wx22与Wx11合并为第三小波变换系数。
可选地，所述S42包括：
对所述第二降噪处理序列Wx2每一通道每一尺度小波变换的小波系数序列Cxni进行如下处理：
A1，根据其序列长度L将其平均分为M个区间，每一区间序列的长度为β＝floor(L/M)，计算每一区间小波系数的熵值H，
若该小区间标记为低小波系数区间，否，则标记为高小波系数区间；
A2，对Cxni小波系数区间中所有高小波系数区间中的点按照以下公式进行降噪处理：
其中，Hxni(k，w)为第k个高小波系数区间中的第w个小波系数值，q是预设小波熵高调整阈值，H1为Cxni中所有高小波系数区间的小波系数绝对值和的平均数；
A3，对Cxni小波系数区间中所有低小波系数区间中的点按照以下公式进行降噪处理：
其中，Lxni(k，w)为Cxni小波系数区间中第k个低小波系数区间中的第w个小波系数值，p是预设小波熵低调整阈值，L1为Cxni中所有低小波系数区间的小波系数绝对值和的平均数；
A5，将处理过的Cxni小波系数代替原来的Cxni小波系数。
第二方面，本发明提供一种计算机装置，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任一项所述膈肌肌电信号降噪方法的步骤。
第三方面，本发明提供一种膈肌肌电信号采集装置，包括采集器、处理器和输出单元，所述采集器用于获取所述膈肌肌电信号和所述参考心电噪声信号，所述处理器用于执行上述任一项所述膈肌肌电信号降噪方法；所述输出单元用于输出经过降噪处理的无心电噪声膈肌肌电信号。
可选地，所述采集器采用体表电极获取膈肌肌电信号。
可选地，所述采集器包括4组体表电极，对称放置于人体体表左右两侧靠近膈肌的位置，用于膈肌肌电信号采集，以及1组体表电极，放置于人体体表左侧靠近心脏的位置，用于参考心电噪声信号采集。
(三)有益效果
与现有技术相比，本发明采用独立成分分析一次降噪叠加小波熵二次降噪的二重降噪技术，克服了现有单独使用独立成分分析法降噪时多个通道降噪效果不均匀且整体降噪效果不佳的问题，能够有效抑制和降低心电噪声对膈肌肌电信号采集的影响，提升膈肌肌电信号采集的有效性和准确性。
另外，本发明在独立成分分析方法中的输入中加入参考心电噪声信号，从而在独立成分分析中能够更有效提取噪声信息，增强降噪效果。
附图说明
图1为本发明降噪方法流程图；
图2为本发明一实施例采集的各道数据；
图3为本发明一实施例各道数据小波变换后的小波系数序列；
图4为本发明一实施例经过独立成分分析分解的源信号；
图5为本发明一实施例经过独立成分分析降噪后的多通道小波系数；
图6为本发明一实施例经过独立成分分析降噪后的多通道数据；
图7为本发明一实施例采用小波熵去噪前后的对比图；
图8为本发明一实施例多通道膈肌肌电去除心电噪声后的数据。
具体实施方式
为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。
实施例一
如图1所示，本实施例提供一种膈肌肌电信号降噪方法，包括如下步骤：
S1，借助于采集装置获取膈肌肌电信号和参考心电噪声信号；
S2，对所述膈肌肌电信号和所述参考心电噪声信号进行小波变换，得到第一小波系数序列；
S3，对所述第一小波变换系数序列利用独立成分分析法进行第一降噪处理，得到第二小波系数序列；
S4，对所述第二小波系数序列利用小波熵法进行第二降噪处理，得到第三小波系数序列；
S5，对所述第三小波系数序列进行小波逆变换，得到降噪的膈肌肌电信号。
为了更好的说明步骤S2，下面对S2中每道信号(包括多个通道的膈肌肌电信号以及1个通道的参考心电噪声信号)进行小波变换的具体步骤进行解释：
S21，将这多道数据X(t)＝[x1(t)；x2(t)；…；xn(t)]中的每道进行5层Harr小波变换，得到每道数据的小波变换系数Cxi(j，k)。
其中，Cxi(j，k)是第i道数据的第j尺度的第k个系数，j＝1，2…，5，k是平移因子，x的取值是a或d，a表示低频成分，d表示高频成分。
计算公式如下：
S22，将分解得到的每道数据的5层离散小波系数进行拼接得到每道数据的小波系数序列，其拼接形式为Wxi＝[Ca5i，Cd5i，Cd4i，Cd3i，Cd2i，Cd1i]，其中Wxi表示第i通道的小波系数序列，是1×N的行向量。Ca5i是第i通道第5尺度的低频小波系数，Cdni是第i通道第n尺度的高频小波系数，n＝1，2，…，5。
为了详细说明S3的工作步骤，下面将S3分解为S31至S35的子步骤，包括：
S31，将所述第一小波系数序列Wx＝[Wx1；Wx2；…；Wxn]作为独立成分分析算法的输入，经过独立成分分析得到源信号矩阵Y＝[y1；y2；…；yn]，以及解混矩阵W和混合矩阵A，
其中，Wx和Y为n×N的矩阵，W与A都是n×n的矩阵，N为小波系数向量长度，n为所有信号数量，包括若干个通道的膈肌肌电信号和参考心电噪声信号；
S32，将得到的各个源信号序列yi与参考心电噪声对应的小波变换系数序列Wx噪声做相关运算得到相关系数的绝对值θi，计算公式如下：
θi的取值范围是[0，1]，作为相关性的判据，其物理意义是其绝对值越大说明与噪声源信号越相似。
S33，根据预设相关度判定阈值T1，按照以下公式处理各个源信号序列yi：
S34，将处理过的yi放入对应Y中的第i行，得到处理后的源数据矩阵Y1；
S35，将Y1用混合矩阵A进行独立成分分析反变换，得到第一降噪后的各通道信号的第二小波变换系数序列Wx1，计算公式如下：
Wx1＝A*Y1。
其中，独立成分分析(ICA)有多种实现方法，本发明采用基于负熵判据的FASTICA算法，该算法采用分布式和并行的计算方法，具有运行速度快、内存占用量小的特点。
为了更好的说明S4，将S4的工作方法分解为子步骤S41至S43，包括以下内容：
S41，对第二小波变换系数序列根据预设噪声指标，选取超过预设噪声指标的通道组成第二降噪处理序列Wx2，其他通道组成Wx11；
S42，将第二降噪处理序列Wx2利用小波熵法进行第二降噪，得到Wx22；
S43，将第二降噪后的Wx22与Wx11合并为第三小波变换系数。
其中，S42包括：
对所述第二降噪处理序列Wx2每一通道每一尺度小波变换的小波系数序列Cxni进行如下处理：
A1，根据其序列长度L将其平均分为M个区间，每一区间序列的长度为β＝floor(L/M)，计算每一区间小波系数的熵值H，
需要说明的是，小波系数的熵值H的计算可以采用如下步骤：
AA1，算出Cxni系数中每一小等份中小波系数绝对值的和，其中第j小等份的小波系数绝对值和的计算公式如下：
AA2，计算Cxni小波系数中每一小等份系数绝对值的和占Cxni小波系数总体绝对值和的比例，即概率。其中Cxni小波系数第j小等份的概率为Pi，x，n(j)，计算公式如下：
其中，j＝1，2，…，M。
AA3，计算Cxni小波系数中每一小等份占整个N等份的Cxni小波系数中的熵值，其中第j小等份系数占整个N等份的Cxni小波系数中的熵值，计算公式如下：
Hi，x，n(j)＝-Pi，x，n(j)*log2(Pi，x，n(j))；
其中，j＝1，2，…，M。
根据计算出的每一区间小波系数的熵值H，若该小区间标记为低小波系数区间，否，则标记为高小波系数区间；
A2，对Cxni小波系数区间中所有高小波系数区间中的点按照以下公式进行降噪处理：
其中，Hxni(k，w)为第k个高小波系数区间中的第w个小波系数值，q是预设小波熵高调整阈值，H1为Cxni中所有高小波系数区间的小波系数绝对值和的平均数；
A3，对Cxni小波系数区间中所有低小波系数区间中的点按照以下公式进行降噪处理：
其中，Lxni(k，w)为Cxni小波系数区间中第k个低小波系数区间中的第w个小波系数值，p是预设小波熵低调整阈值，L1为Cxni中所有低小波系数区间的小波系数绝对值和的平均数；
A5，将处理过的Cxni小波系数代替原来的Cxni小波系数。
本发明实施例提供的以上膈肌肌电信号降噪方法，首先通过独立成分分析(ICA)方法一次降噪分离出大部分通道的心电噪声，在此基础上，对于残留心电噪声比较强或者心电噪声减弱小的通道采用小波熵的方法将每一尺度的小波系数划分为高小波系数熵与低小波系数熵的区间，对这两个区间采用不同的小波阈值来进行二次降噪，经过二重降噪处理的膈肌肌电信号能够达到良好的降噪效果，有效提升膈肌肌电信号采集和测试的准确性、有效性。
需要说明的是，本发明实施例膈肌肌电信号的降噪方法，是根据膈肌肌电信号及其噪声的特点设计的，并不适用于所有人体测量信号，比如脑电信号的降噪处理。因为体表膈肌肌电属于肌电信号，幅度较大，量级是V，主要的频率范围20～250Hz，其噪声来源主要是心电信号；而脑电信号很微弱，幅度的量级是uV，频率范围是0～60Hz，具有非平稳、非线性、随机性强的特点，容易受到包括心电、眼电、肌电和其他生理性伪迹在内的多种噪声干扰。信号形态的明显差异使得两种信号降噪关注的重点和采用的具体方法都有差异，两者降噪方法不能通用。
实施例二
本实施例提供了一种具体利用本发明膈肌肌电信号降噪方法进行膈肌信号采集和降噪的应用实例。
本实施例中，采集装置采用4组体表电极采集膈肌肌电信号，1组体表电极采集参考心电噪声信号，其中采集膈肌肌电信号的电极在人体体表左右两侧靠近膈肌的位置对称放置，采集参考心电噪声信号的电极放置在人体体表左侧靠近心脏的位置。
各通道数据的采样频率为2.5KHz，将电极采集的数据分别经过截止频率为15Hz的高通滤波，截止频率500Hz的低通滤波，以及相应的频谱插值，得到的初始采集数据，如图2所示，一共包含5道数据，前四道数据从上到下依次对应电极采集的膈肌肌电信号data11、data21、data31、data41，各道数据的横轴为采样的时间点数，纵轴为测得信号的电压幅度，单位为V。其中data21和data31是在人体体表左侧(即靠近心脏的位置)电极采集的信号，心电噪声较强。而data11与data41是在人体体表右侧采集到的信号，离心脏较远，心电噪声相对较弱。第5道数据是心电参考噪声。
将这5道数据作为本发明降噪方法的输入，对图2中的数据分别进行小波变换得到各道数据的小波变换系数，如图3所示。
在图3中可以看到data21与data31的心电噪声成分是比较强的，有很高的尖线，data11与data41的心电噪声相对较弱。
将图3的数据输入到ICA算法中得到各个源成分，如图4所示。
将图4中的各个源成分信号与参考噪声的小波系数序列做相关运算，找到与心电噪声关联较大的源成分并去除。将去除掉噪声源后的源成分数据用混合矩阵进行混合，得到降噪后的各通道小波系数，如图5。
为了验证ICA算法的效果，将ICA处理过的小波系数经过小波反变换观察一下降噪后的时域信号，参考图6可以看到，data21与data31中的心电噪声的幅度大幅度降低，信号的幅度也发生较大的改变，幅度由原来的心电噪声占据主导变为现在膈肌肌电占主导，信噪比有了较大的改善。data11与data41的幅度相对变化较小。
为了评估ICA算法的降噪效果可以采用信号的中心频率和均方值作为评价指标，将ICA降噪前的多通道时域数据与ICA降噪后的数据分别计算信号的中频和均方值计算，得到表1与表2中的数据。
表1降噪前各通道数据的性能指标
| data11 | data21 | data31 | data41 | |
| fm(中频Hz) | 132.05 | 100.10 | 98.59 | 136.81 |
| RMS(均方值) | 0.0974 | 0.1954 | 0.3123 | 0.1138 |
表2降噪后各通道数据的性能指标
| data11 | data21 | data31 | data41 | |
| fm(中频Hz) | 132.99 | 132.89 | 127.80 | 137.98 |
| RMS(均方值) | 0.0949 | 0.0746 | 0.0562 | 0.1115 |
表1是ICA降噪前各通道数据计算得到的均方值和中频值，表2是ICA降噪后各通道数据的均方值与中频值。由于采用了ICA降噪的方法，多通道数据信号整体的均方值在降噪后会降低。其中data31均方值下降幅度大，降噪后的均方值约为原来的18％。data21的均方值下降幅度次之，降噪后的均方值变为原来的38.2％。而data11与data41的均方值下降程度相对较低，均方值下降幅度仅为原来的2％左右。这说明ICA降噪的方法对data21与data31的数据效果明显，对data11与data41的数据效果不是很明显，但也会去除小部分的心电噪声。
由于心电信号的频谱集中在20～50Hz，当采用ICA降噪算法后，原始信号的频谱在该频率范围内大幅下降，会导致其中频值上升。从表1、2的数据来看，仍然是data21与data31的效果显著，中频值变化幅度大。data21降噪后中频值上升为原来的32.76％，data31降噪后中频值上升为原来的29.62％。data11的中频值上升原来的0.7％，data41的中频值上升为原来的0.8％。从这些数据可以看到data21、data31心电去噪效果较好，data21、data31的均方值幅度下降大，中频值上升幅度高，去燥效果佳。而data11与data41的均方值下降幅度小、中频值增幅小，去噪效果弱，但也去除了一部分的心电噪声。之所以data11与data41去燥效果弱一些，是因为data21与data31是电极放置在人体左侧采集的信号，心电噪声幅度大、相关性强，采用ICA方法容易分离出源信号。而data11与data41是位于人体右侧的电极采集的数据，心电噪声幅度小，相比于人体左侧采集的数据而言，心电噪声在幅度和波形上与左侧的信号略有不同，导致ICA对这两道数据中心电噪声分离效果不佳，去除少。
针对这种情况，采用小波熵阈值处理的方式对data11与data41两道数据经过ICA处理后的各尺度系数进行分类，将各尺度系数分为高小波系数区间和低小波系数区间，对不同幅度的小波系数采用不同的阈值进行阈值降噪，得到处理后的结果。
如图7所示，该图是data11采用小波熵阈值去噪的方法得到的去噪结果前后对比图，其中该图下方的波形是data11经过ICA方法处理后的数据，即采用小波熵阈值降噪前的数据，该图上方的波形是data11经过ICA处理后再经过小波熵阈值降噪后得到的结果。可以看到经过小波熵去噪后data11中的心电噪声有了明显的消除，整体效果较好。
表3小波熵阈值去噪后信号的性能指标
| data11 | data41 | |
| fm(中频Hz) | 154.67 | 158.86 |
| RMS(均方值) | 0.0699 | 0.0809 |
表3是经过小波熵阈值处理前后的data11与data41的性能指标的值。ICA方法处理后的data11与data41再经过小波熵阈值降噪后，数据的中频得到明显的提升，data11中频提升为原来的16.3％(与表2中的数据相比)，data41的中频提升为原来的15.13％。且data11与data41的均方值也有明显的下降，data11均方值下降为原来的26.3％，data41均方值下降为原来的27.4％。这说明采用小波熵阈值去噪可以对心电噪声进行较好的去噪处理。
四通道的膈肌肌电信号经过上述二重降噪过程后的数据如图8所示，从图8中可以看到四道膈肌肌电数据经过ICA和小波熵处理后起到了很好的降噪效果。
实施例三
本实施例提供一种计算机装置，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上述实施例一中任一项所述膈肌肌电信号降噪方法的步骤。
实施例四
本实施例提供一种膈肌肌电信号采集装置，包括采集器、处理器和输出单元，所述采集器用于获取所述膈肌肌电信号和所述参考心电噪声信号，所述处理器用于执行上述任一项所述膈肌肌电信号降噪方法；所述输出单元用于输出经过降噪处理的无心电噪声膈肌肌电信号；所述采集器包括4组体表电极，对称放置于人体体表左右两侧靠近膈肌的位置，用于膈肌肌电信号采集，以及1组体表电极，放置于人体体表左侧靠近心脏的位置，用于参考心电噪声信号采集。
由于本发明上述实施例所描述的系统/装置，为实施本发明上述实施例的方法所采用的系统/装置，故而基于本发明上述实施例所描述的方法，本领域所属技术人员能够了解该系统/装置的具体结构及变形，因而在此不再赘述。凡是本发明上述实施例的方法所采用的系统/装置都属于本发明所欲保护的范围。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。
应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。本发明可以借助于包括有若干不同部件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的权利要求中，这些装置中的若干个可以是通过同一个硬件来具体体现。词语第一、第二、第三等的使用，仅是为了表述方便，而不表示任何顺序。可将这些词语理解为部件名称的一部分。
此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。
Claims (10)
1.一种膈肌肌电信号降噪方法，其特征在于，包括：
S1，借助于采集装置获取膈肌肌电信号和参考心电噪声信号；
S2，对所述膈肌肌电信号和所述参考心电噪声信号进行小波变换，得到第一小波系数序列；
S3，对所述第一小波变换系数序列利用独立成分分析法进行第一降噪处理，得到第二小波系数序列；
S4，对所述第二小波系数序列利用小波熵法进行第二降噪处理，得到第三小波系数序列；
S5，对所述第三小波系数序列进行小波逆变换，得到降噪的膈肌肌电信号。
2.根据权利要求1所述的一种膈肌肌电信号降噪方法，其特征在于，所述S2包括：
对所述膈肌肌电信号和所述参考心电噪声信号分别进行5层Harr小波变换，将变换后的各信号的小波系数合并为所述第一小波系数序列。
3.根据权利要求1所述的一种膈肌肌电信号降噪方法，其特征在于，所述S3包括：
S31，将所述第一小波系数序列Wx＝[Wx1；Wx2；…；Wxn]作为独立成分分析算法的输入，经过独立成分分析得到源信号矩阵Y＝[y1；y2；…；yn]，以及解混矩阵W和混合矩阵A，
其中，Wx和Y为n×N的矩阵，W与A都是n×n的矩阵，N为小波系数向量长度，n为所有信号数量，包括若干个通道的膈肌肌电信号和参考心电噪声信号；
S32，将得到的各个源信号序列yi与参考心电噪声对应的小波变换系数序列Wx噪声做相关运算得到相关系数的绝对值θi，计算公式如下：
S33，根据预设相关度判定阈值T1，按照以下公式处理各个源信号序列yi：
S34，将处理过的yi放入对应Y中的第i行，得到处理后的源数据矩阵Y1；
S35，将Y1用混合矩阵A进行独立成分分析反变换，得到第一降噪后的各通道信号的第二小波变换系数序列Wx1，计算公式如下：
Wx1＝A*Y1。
4.根据权利要求1所述的一种膈肌肌电信号降噪方法，其特征在于：
对所述第一小波变换系数序列利用FASTICA算法执行独立成分分析和第一降噪处理。
5.根据权利要求1所述的一种膈肌肌电信号降噪方法，其特征在于，所述S4包括：
S41，对第二小波变换系数序列根据预设噪声指标，选取超过预设噪声指标的通道组成第二降噪处理序列Wx2，其他通道组成Wx11；
S42，将第二降噪处理序列Wx2利用小波熵法进行第二降噪，得到Wx22；
S43，将第二降噪后的Wx22与Wx11合并为第三小波变换系数。
6.根据权利要求5所述的一种膈肌肌电信号降噪方法，其特征在于，所述S42包括：
对所述第二降噪处理序列Wx2每一通道每一尺度小波变换的小波系数序列Cxni进行如下处理：
A1，根据其序列长度L将其平均分为M个区间，每一区间序列的长度为β＝floor(L/M)，计算每一区间小波系数的熵值H，
若该小区间标记为低小波系数区间，否，则标记为高小波系数区间；
A2，对Cxni小波系数区间中所有高小波系数区间中的点按照以下公式进行降噪处理：
其中，Hxni(k，w)为第k个高小波系数区间中的第w个小波系数值，q是预设小波熵高调整阈值，H1为Cxni中所有高小波系数区间的小波系数绝对值和的平均数；
A3，对Cxni小波系数区间中所有低小波系数区间中的点按照以下公式进行降噪处理：
其中，Lxni(k，w)为Cxni小波系数区间中第k个低小波系数区间中的第w个小波系数值，p是预设小波熵低调整阈值，L1为Cxni中所有低小波系数区间的小波系数绝对值和的平均数；
A5，将处理过的Cxni小波系数代替原来的Cxni小波系数。
7.一种计算机设备，其特征在于，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如权利要求1至6任一项所述膈肌肌电信号降噪方法的步骤。
8.一种膈肌肌电信号采集装置，其特征在于，包括采集器、处理器和输出单元：
所述采集器用于获取所述膈肌肌电信号和所述参考心电噪声信号；
所述处理器用于执行权利要求1至6任一项所述膈肌肌电信号降噪方法；
所述输出单元用于输出经过降噪处理的无心电噪声膈肌肌电信号。
9.根据权利要求8所述的一种膈肌肌电信号采集装置，其特征在于，所述采集器采用体表电极获取膈肌肌电信号和参考心电噪声信号。
10.根据权利要求9所述的一种膈肌肌电信号采集装置，其特征在于，所述采集器包括：
4组体表电极，对称放置于人体体表左右两侧靠近膈肌的位置，用于膈肌肌电信号采集；
1组体表电极，放置于人体体表左侧靠近心脏的位置，用于参考心电噪声信号采集。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310353452.7A CN116584960A (zh) | 2023-04-04 | 2023-04-04 | 一种膈肌肌电信号降噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310353452.7A CN116584960A (zh) | 2023-04-04 | 2023-04-04 | 一种膈肌肌电信号降噪方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116584960A true CN116584960A (zh) | 2023-08-15 |
Family
ID=87605166
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310353452.7A Pending CN116584960A (zh) | 2023-04-04 | 2023-04-04 | 一种膈肌肌电信号降噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116584960A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117017323A (zh) * | 2023-09-14 | 2023-11-10 | 中国科学技术大学 | 基于盲源分离的高密度表面膈肌肌电采集与预处理方法 |
-
2023
- 2023-04-04 CN CN202310353452.7A patent/CN116584960A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117017323A (zh) * | 2023-09-14 | 2023-11-10 | 中国科学技术大学 | 基于盲源分离的高密度表面膈肌肌电采集与预处理方法 |
| CN117017323B (zh) * | 2023-09-14 | 2024-03-29 | 中国科学技术大学 | 基于盲源分离的高密度表面膈肌肌电采集与预处理方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108714026B (zh) | 基于深度卷积神经网络和在线决策融合的细粒度心电信号分类方法 | |
| CN106709469B (zh) | 基于脑电和肌电多特征的自动睡眠分期方法 | |
| US7809433B2 (en) | Method and system for limiting interference in electroencephalographic signals | |
| CN109907752B (zh) | 一种去除运动伪影干扰与心电特征检测的心电诊断与监护系统 | |
| Kania et al. | Wavelet denoising for multi-lead high resolution ECG signals | |
| CN110680308B (zh) | 基于改进emd与阈值法融合的心电信号去噪方法 | |
| CN107260166A (zh) | 一种实用化在线脑电伪迹剔除方法 | |
| CN103961092B (zh) | 基于自适应阈值处理的脑电信号去噪方法 | |
| CN102697493A (zh) | 一种快速的脑电信号中眼电伪迹自动识别和去除的方法 | |
| CN102835955A (zh) | 一种无需设定阈值的脑电信号中眼电伪迹自动去除的方法 | |
| CN109589114A (zh) | 基于ceemd和区间阈值的肌电消噪方法 | |
| CN106236080B (zh) | 基于多通道的脑电信号中肌电噪声的消除方法 | |
| CN109948396B (zh) | 一种心拍分类方法、心拍分类装置及电子设备 | |
| Wu et al. | EMGdi signal enhancement based on ICA decomposition and wavelet transform | |
| Abbaspour et al. | Evaluation of wavelet based methods in removing motion artifact from ECG signal | |
| CN106419898A (zh) | 一种去除心电信号基线漂移的方法 | |
| CN109480832A (zh) | 一种单通道的脑电信号中肌电伪迹的消除方法 | |
| CN116584960A (zh) | 一种膈肌肌电信号降噪方法 | |
| Tang et al. | ECG de-noising based on empirical mode decomposition | |
| CN111631710A (zh) | 一种状态相关的动态脑电信号中肌电伪迹的消除方法 | |
| Joy et al. | Wavelet based EMG artifact removal from ECG signal | |
| Haibing et al. | Discrete wavelet soft threshold denoise processing for ECG signal | |
| Ozkaya et al. | Frequency analysis of EMG signals with Matlab sptool | |
| Nougarou et al. | Efficient procedure to remove ECG from sEMG with limited deteriorations: Extraction, quasi-periodic detection and cancellation | |
| Liu et al. | Remove motion artifacts from scalp single channel EEG based on noise assisted least square multivariate empirical mode decomposition |
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 |