CN113066502A - 基于vmd和多小波的心音分割定位方法 - Google Patents

基于vmd和多小波的心音分割定位方法 Download PDF

Info

Publication number
CN113066502A
CN113066502A CN202110263717.5A CN202110263717A CN113066502A CN 113066502 A CN113066502 A CN 113066502A CN 202110263717 A CN202110263717 A CN 202110263717A CN 113066502 A CN113066502 A CN 113066502A
Authority
CN
China
Prior art keywords
heart sound
signal
decomposition
matrix
interval
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.)
Granted
Application number
CN202110263717.5A
Other languages
English (en)
Other versions
CN113066502B (zh
Inventor
夏侯士戟
梁宇航
马敏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202110263717.5A priority Critical patent/CN113066502B/zh
Publication of CN113066502A publication Critical patent/CN113066502A/zh
Application granted granted Critical
Publication of CN113066502B publication Critical patent/CN113066502B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; SPEECH OR AUDIO CODING OR DECODING
    • G10L19/00—Speech or audio signals analysis-synthesis techniques for redundancy reduction, e.g. in vocoders; Coding or decoding of speech or audio signals, using source filter models or psychoacoustic analysis
    • G10L19/02—Speech or audio signals analysis-synthesis techniques for redundancy reduction, e.g. in vocoders; Coding or decoding of speech or audio signals, using source filter models or psychoacoustic analysis using spectral analysis, e.g. transform vocoders or subband vocoders
    • G10L19/0212—Speech or audio signals analysis-synthesis techniques for redundancy reduction, e.g. in vocoders; Coding or decoding of speech or audio signals, using source filter models or psychoacoustic analysis using spectral analysis, e.g. transform vocoders or subband vocoders using orthogonal transformation
    • G10L19/0216—Speech or audio signals analysis-synthesis techniques for redundancy reduction, e.g. in vocoders; Coding or decoding of speech or audio signals, using source filter models or psychoacoustic analysis using spectral analysis, e.g. transform vocoders or subband vocoders using orthogonal transformation using wavelet decomposition
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B7/00—Instruments for auscultation
    • A61B7/02—Stethoscopes
    • A61B7/04—Electric stethoscopes
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Processing of the speech or voice signal to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Processing of the speech or voice signal to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0272—Voice signal separating
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/03—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/27—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the analysis technique
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/48—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use
    • G10L25/51—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use for comparison or discrimination

Abstract

本发明公开了一种基于VMD和多小波的心音分割定位方法，首先对心音信号样本中心音区间和类型进行标记，采用滑动窗对心音信号样本进行片段划分，使用峭度对VMD分解得到的各个IMF分量进行挑选，对最有效的IMF分量进行GHM多小波包分解，基于分解结果构建心音信号样本的时频域矩阵，根据标记的心音区间从中提取子矩阵并得到时间频率列向量，将其作为训练样本对预设的分类模型进行训练；对待分割定位的心音信号采用相同方法得到时频域矩阵，提取香农能量包络，基于最大类间方差法确定心音区间，提取子矩阵并采用相同方法获取时刻频率列向量，输入训练好的分类模型得到心音信号分割定位结果。采用本发明可以提高高噪声下的心音分割性能。

Description

基于VMD和多小波的心音分割定位方法
技术领域
本发明属于心音处理技术领域，更为具体地讲，涉及一种基于VMD和多小波的心音分割定位方法。
背景技术
目前采用可穿戴心音检测设备进行心音监测已经成为了诊断和控制心脏疾病的重要手段，但是现有的可穿戴心音检测设备所采集的心音信号具有低幅值、低频率、易受干扰等特性，考虑到可穿戴心音检测设备使用的复杂环境，使用过程中可能会引入其他生理音干扰、背景噪声、工频干扰以及运动伪影等，进而会严重影响对所采集到心音信号做进一步的处理。
传统的心音分割定位算法是与心电图相结合，在使用心电图作为参考的情况下，心音分割的效果很好。然而，心电图作为另一种信号源，在医学检查中可能不方便使用。如果嵌入ECG电路会极大地提高电子听诊器的硬件复杂度。因此，使用心音信号作为唯一信号源也是当下的主要研究方向。目前不使用参考信号的主流心音信号处理流程为先对原始心音信号进行去噪之后，对其进行特征提取、分类，从而对心动周期的各组分进行分割定位。
在信号去噪方面，最传统的方法是基于心音信号的频域特征进行滤波去噪，但这种方法仅可抑制心音分量频带外的干扰，心音作为典型的非平稳信号，非常容易受带内噪声或其他生理音干扰，由于时频率(TF)域可以产生更健壮的定位和分类方法，特别是对非平稳信号，因此时频域的分析方法也被应用于了心音处理领域。
目前，研究非平稳信号常用的时频分析工具有：短时傅里叶变换、维格纳-维尔分布、科恩类、小波变换和S变换等。时傅里叶变换的时频分辨率不能自适应，而维格纳-维尔分布和科恩类虽然具有良好的时频特性，但是由于存在交叉干扰项，影响了它们的实际应用范围。作为短时傅里叶变换的继承和发展，S变换采用高斯窗函数且窗宽与频率的倒数成正比，免去了窗函数的选择和改善了窗宽固定的缺陷，并且时频表示中各频率分量的相位谱与原始信号保持直接的联系，使其在信号分析中可以采用更多的特征量。S变换的缺点是在较高频带范围内相对于小波变换来说，对频域的分析上不够精确。
另外，因为可穿戴心音监测平台采集的单通道测量信号是由不同的物理过程产生的，而且它们在统计上是相互独立的，所以可以看作心音信号和背噪两个独立信号的叠加，因此采用基于盲源分离的方法进行处理也是现在的一种新的处理思路。目前使用频率比较高的盲源信号分离的方法主要是经验模态分解(Empirical Mode Decomposition,EMD)、小波分解以及他们的衍生算法如(Ensemble Empirical Mode Decomposition,EEMD)等。然而EMD算法存在模态混叠、模态残差噪声及虚假模态问题，而且EMD算法现实中所产生的本征模函数(Intrinsic Oscillatory Mode,IMF)不会保持完全稳定的频率和振幅。因此在2014年提出了VMD(Variational Mode Decomposition，变分模态分解)的方法，此方法是一种自适应、完全非递归的模态变分和信号处理的方法。该技术具有可以确定模态分解个数的优点它克服了EMD方法存在端点效应和模态分量混叠的问题，并且具有更坚实的数学理论基础。
在心音特征提取方面，传统的方法所提取的特征大多以频域特征和时域特征为主，或者对信号提取包络特征，即利用各种技术构造心音包络来进行心音分割，其中香农包络和希尔伯特包络是最常使用的两种包络。基于包络抽取方法参数的推导直接受到阈值的影响，阈值选择主要是依靠经验值，所以容易因为过学习而导致在不同样本集之间泛化性不强。
总体来说，虽然目前方法可以实现在正常心音信号中识别S1/S2的难度不高，但如果想在可穿戴平台上使用传统的心音处理方法进行处理还是有很多问题：可穿戴平台面对的复杂噪声的使用场景，如杂音、咔嚓声、噼啪声和低信噪比等场景，会极大地提高心音分割定位的难度。
发明内容
本发明的目的在于克服现有技术的不足，提供一种基于VMD和多小波的心音分割定位方法，结合VMD分解、多小波、香农包络和奇异值分解，以提高在高噪声的现实场景中识别和分割心音信号的性能。
为了实现上述发明目的，本发明基于VMD和多小波的心音分割定位方法包括以下步骤：
S1：采集若干心音信号样本Xd(t)，其中d＝1,2,…,D，D表示心音信号样本的数量，t表示时间；对每个心音信号样本Xd(t)进行去噪处理，然后按照采样频率Fs进行重采样得到每个心音信号样本的采样信号Xd(n)，其中n＝1,2,…,Sd，Sd表示采样信号Xd(n)的采样点数；在每个心音信号样本的采样信号Xd(n)中标记第一心音S1和第二心音S2的区间，记标记得到的第k个心音区间为
Figure BDA0002971155300000031
其中
Figure BDA0002971155300000032
分别表示第d个心音信号样本的采样信号Xd(n)中第k个心音区间的起点和终点，k＝1,2,…,Kd，Kd表示采样信号xd(n)中心音区间的数量，并根据心音类型设置其标签flagd,k，flagd,k＝1,2；
S2：采用滑动窗对各个采样信号Xd(n)进行片段划分得到Nd个采样信号片段Xd,i(n′)，i＝1,2,…,Nd，n′＝1,2,…,w，w表示采样信号Xd(n)每个信号片段的长度；
S3：对于每个采样信号Xd(n)，在步骤S2划分得到Nd个信号片段Xd,i(n′)后，对每个信号片段分别进行V层VMD分解，记信号片段Xd,i(n′)经VMD分解得到的第v个分量为
Figure BDA0002971155300000033
S4：对于每个信号片段Xd,i(n′)经过VMD分解得到的V个分量
Figure BDA0002971155300000034
分别计算信号片段Xd,i(n′)和每个分量
Figure BDA0002971155300000035
的峭度值，选择峭度值与信号片段Xd,i(n′)峭度值的L2范数最大的分量作为信号片段Xd,i(n′)对应的IMF分量信号
Figure BDA0002971155300000036
S5：将各个IMF分量信号
Figure BDA0002971155300000037
插值到预设长度W，得到信号xd,i(n″)，n″＝1,2,…,W，W为2λ的整数倍值；然后对每个信号片段xd,i(n″)进行预滤波处理，得到大小为2×W的预滤波信号
Figure BDA0002971155300000038
采用GHM多小波系统对每个预滤波信号
Figure BDA0002971155300000039
进行分解，记分解得到的信号分量数量为K，记第j个信号分量为yd,i,j，j＝1,2,…,K，信号分量yd,i,j的大小为2×W′，W′＝W/2λ，λ表示GHM多小波系统的分解层数；
S6：对于各个信号xd,i(n″)经GHM多小波分解得到K个大小为2×W′信号分量yd,i,j，构建得到大小为2K×W′的信号矩阵Zd,i：
Figure BDA00029711553000000310
对于同一心音信号样本的Nd个信号片段xd,i(n′)，采用其所对应的Nd个信号矩阵Zd,i构建得到大小为2K×Ld的时频域矩阵Md：
Figure BDA0002971155300000041
其中，
Figure BDA0002971155300000042
Figure BDA0002971155300000043
表示时频域矩阵Md中坐标为(α,β)的元素值，其中α＝1,2,…2K，β＝1,2,…,Ld；令
Figure BDA0002971155300000044
元素
Figure BDA0002971155300000045
按照以下公式确定：
Figure BDA0002971155300000046
其中，
Figure BDA0002971155300000047
分别表示向下取整和向上取整；
S7：对于得到的每个时频域矩阵Md对行进行降维，得到降维后大小为P×Ld的时频域矩阵
Figure BDA0002971155300000048
其中P根据实际需要设置；
S8：对步骤S1中每个心音信号样本所标记的心音区间
Figure BDA0002971155300000049
进行展宽得到区间
Figure BDA00029711553000000410
其中γd＝[Ld/Nd]，然后根据区间
Figure BDA00029711553000000411
从时频域矩阵
Figure BDA00029711553000000412
提取对应列构成大小为P×Ud的子矩阵
Figure BDA00029711553000000413
其中
Figure BDA00029711553000000414
对各个子矩阵
Figure BDA00029711553000000415
进行奇异值分解，其分解结果采用如下公式表示：
Figure BDA00029711553000000416
其中，Qd,k表示P×P阶的酉矩阵，Vd,k表示Ud×Ud阶的酉矩阵，∑d,k表示大小为P×Ud的奇异值矩阵；
在奇异值矩阵∑d,k选择最大的前M个奇异值，其中M根据实际需要设置，从矩阵Qd,k中提取出M个坐标为(m,m)的元素值
Figure BDA00029711553000000417
从矩阵Vd,k中提取出M个坐标为(m,m)的元素值
Figure BDA00029711553000000418
其中m＝1,2,…,M，然后拼接得到时间频率列向量Fd,k：
Figure BDA00029711553000000419
S9：将步骤S8得到的各个时间频率列向量Fd,k作为输入，对应的标签flagd,k作为输出，对预设的分类模型进行训练；
S10：采集待分割定位的心音信号
Figure BDA00029711553000000420
然后按照预设采样频率Fs进行采样得到采样信号
Figure BDA00029711553000000421
其中
Figure BDA00029711553000000422
N表示采样信号
Figure BDA00029711553000000423
所包含的采样点数；
S11：采用步骤S2中的相同方法对待分割定位心音信号的采样信号
Figure BDA00029711553000000424
进行片段划分，得到N个采样信号片段
Figure BDA0002971155300000051
采用步骤S3中的相同方法对每个信号片段
Figure BDA0002971155300000052
进行VMD分解，采用步骤S4中的相同方法筛选得到IMF分量信号
Figure BDA0002971155300000053
采用步骤S5中的相同方法对每个IMF分量信号
Figure BDA0002971155300000054
进行插值、预滤波处理和GHM多小波分解，得到K个大小为2×W′的信号分量yi′,j，采用步骤S6中的相同方法构建得到大小为2K×L的时频域矩阵M′，
Figure BDA0002971155300000055
S12：对于待分割定位心音信号的时频域矩阵M′，分别计算每一列的香农能量SSE(b)，从而得到香农能量包络；
S13：对于待分割定位心音信号所对应的香农能量包络，筛选出香农能量最大值SSEmax，然后基于最大类间方差法获取区分噪音和心音的分类阈值TH，然后使用类间方差法确定的分类阈值TH作为较高的门限T1，即T1＝TH，在香农能量包络中提取出连续高于阈值TH的区间，将区间长度小于预设阈值的区间删除，剩余区间即为第一心音S1或第二心音S2所在的粗定心音区间
Figure BDA0002971155300000056
K′表示待分割定位心音信号的心音区间数量；然后再在TH的基础上确定一个较低的阈值T2，令T2＝TH-η，η为预设的大于0的偏差值，在Zk″_temp的基础上向两旁扩展搜索，分别找到香农能量包络与阈值T2相交的两个端点
Figure BDA0002971155300000057
作为心音的正确起始端点，得到心音区间
Figure BDA0002971155300000058
S14：根据步骤S12得到的心音区间Zk′，从时频域矩阵M′提取对应列构成子矩阵Mk″，采用步骤S8中的相同方法，基于奇异值分解得到该心音区间的时间频率列向量Fk″；
S15：将每个心音区间的时间频率列向量Fk″分别输入步骤S9训练好的分类模型，其分类结果即为该心音区间内心音信号的类型；
S16：将待分割定位心音信号的每个心音区间
Figure BDA0002971155300000059
还原得到区间
Figure BDA00029711553000000510
γ＝[L/N]，然后在采样信号x′(n)标记每个区间Zk′，并根据步骤S15中的分类结果确定该区间内心音信号的类型，从而得到心音信号的分割定位结果。
本发明基于VMD和多小波的心音分割定位方法，首先对心音信号样本中第一心音和第二心音的心音区间和类型进行标记，采用滑动窗对心音信号样本进行片段划分，使用峭度对VMD分解得到的各个IMF分量进行挑选，对最有效的IMF分量进行GHM多小波包分解，基于分解结果构建心音信号样本的时频域矩阵，根据标记的心音区间从中提取子矩阵并得到时间频率列向量，将其作为输入，对应类型标签作为输出，对预设的分类模型进行训练；对待分割定位的心音信号采用相同方法得到时频域矩阵，提取香农能量包络，基于最大类间方差法确定心音区间，提取子矩阵并采用相同方法获取时刻频率列向量，输入训练好的分类模型得到分类结果，从而得到心音信号分割定位结果。
本发明具有以下有益效果：
1)本发明可以实现心音信号的自适应分解、自适应降维、全自动提取、自动选择分解层数，无需人工选择，避免人工参数设置干扰；
2)经实验验证，本发明能够在高噪声环境下对心音信号的各组份进行鲁棒的准确定位和分类。
附图说明
图1是本发明基于VMD和多小波的心音分割定位方法的具体实施方式流程图；
图2是心音信号的VMD分解结果示例图；
图3是小波包分解示意图；
图4是本实施例中时频域矩阵的示例图；
图5是本实施例中自行采集心音信号示例图；
图6是本实施例中纯净心音信号叠加白噪声所得到的归一化心音信号波形图；
图7是VMD+EWT算法采用对图6所示心音信号进行处理后的归一化信号波形图；
图8是采用EMD+PCA算法对图5所示心音信号进行处理后的归一化信号波形图；
图9是采用本发明中多小波+PCA算法对图6所示心音信号进行处理后的归一化信号波形图；
图10是本发明和两种对比方法在不同信噪比高斯白噪声下心音信号去噪后的信噪比比较图；
图11是本发明和两种对比方法在不同信噪比高斯白噪声下心音信号去噪后的相关系数比较图；
图12是本发明和两种对比方法在不同信噪比高斯白噪声下心音信号去噪后的均方误差比较图；
图13是本发明和两种对比方法在不同信噪比粉噪声下心音信号去噪后的信噪比比较图；
图14是本发明和两种对比方法在不同信噪粉噪声下心音信号去噪后的相关系数比较图；
图15是本发明和两种对比方法在不同信噪比粉噪声下心音信号去噪后的均方误差比较图；
图16是本发明和两种对比方法在不同信噪比蓝噪声下心音信号去噪后的信噪比比较图；
图17是本发明和两种对比方法在不同信噪比蓝噪声下心音信号去噪后的相关系数比较图；
图18是本发明和两种对比方法在不同信噪比蓝噪声下心音信号去噪后的均方误差比较图；
图19是本发明和两种对比方法在不同信噪比肺音下心音信号去噪后的信噪比比较图；
图20是本发明和两种对比方法在不同信噪比肺音下心音信号去噪后的相关系数比较图；
图21是本发明和两种对比方法在不同信噪比肺音下心音信号去噪后的均方误差比较图；
图22是本发明和两种对比方法在不同信噪比运动伪迹下心音信号去噪后的信噪比比较图；
图23是本发明和两种对比方法在不同信噪比运动伪迹下心音信号去噪后的相关系数比较图；
图24是本发明和两种对比方法在不同信噪比运动伪迹下心音信号去噪后的均方误差比较图；
图25是采用本发明对图6所示心音信号进行处理后提取的香农包络结果图。
图26是采用本发明对图6所示心音信号进行处理后的处理结果图。
具体实施方式
下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。
实施例
图1是本发明基于VMD和多小波的心音分割定位方法的具体实施方式流程图。如图1所示，本发明心音分割定位方法的具体步骤包括：
S101：心音信号样本预处理：
采集若干心音信号样本Xd(t)，其中d＝1,2,…,D，D表示心音信号样本的数量，t表示时间。对每个心音信号样本Xd(t)进行去噪处理，然后按照采样频率Fs进行重采样得到每个心音信号样本的采样信号Xd(n)，其中n＝1,2,…,Sd，Sd表示采样信号Xd(n)的采样点数。在每个心音信号样本的采样信号Xd(n)中标记第一心音S1和第二心音S2的区间，记标记得到的第k个心音区间为
Figure BDA0002971155300000081
其中
Figure BDA0002971155300000082
分别表示第d个心音信号样本的采样信号Xd(n)中第k个心音区间的起点和终点，k＝1,2,…,Kd，Kd表示采样信号xd(n)中心音区间的数量，并根据心音类型设置其标签flagd,k，flagd,k＝1,2。
对可穿戴平台而言，根据心音信号的频谱分布范围主要集中分布在10-200Hz，高频心杂音可达1000Hz的特性，先进行一个10-500Hz的带通预滤波，同时设置一个50Hz的带阻滤波器来滤除工频干扰。本实施例中在进行去噪处理后，根据香农采样定理，将各个样本统一重采样到2000Hz方便后续进行统一处理。
S102：心音信号样本片段划分：
由于本发明中需要对心音信号进行VMD分解，VMD算法最大的局限性是边界效应和突发的信号，直接的解决方案是将信号分解成较短的块，在这些块上信号足够稳定。因此本发明采用滑动窗对各个采样信号Xd(n)进行片段划分得到Nd个采样信号片段Xd,i(n′)，i＝1,2,…,Nd，n′＝1,2,…,w，w表示采样信号Xd(n)每个信号片段的长度，也就是窗宽。
滑动窗的具体类型的参数可以根据需要进行设置，本实施例中采用窗宽为180，步长为90的高斯窗w(n)对心音信号Xd(n)进行了分段操作，其中高斯窗w(n)的表达式为：
Figure BDA0002971155300000091
其中，δ表示窗宽调节参数，本实施例中设置为2.5，T表示信号周期。
S103：对信号片段进行VMD分解：
对于每个采样信号Xd(n)，在步骤S102划分得到Nd个信号片段Xd,i(n′)后，对每个信号片段分别进行V层VMD分解，记信号片段Xd,i(n′)经VMD分解得到的第v个分量为
Figure BDA0002971155300000092
VMD分解是一种非递归分解方法，自适应同步提取复杂信号中的调幅-调频分量。该方法将信号为一组本征模态分量(IMF)，通过这些本征模态分量，可以重构信号，在频域内迭代更新各个本征模态分量，将本征模态分量频谱的重心作为中心频率。假设本征模态分量在频域内紧致聚集在中心频率周围，即具有稀疏性质。每个本征模态分量的稀疏性质通过它的带宽来描述：(1)应用Hilbert变换构造其解析信号，以便获得非负频率的单边频谱；(2)通过和频率调整至各自中心频率处的指数谐波相乘，将频谱平移至基带；(3)通过梯度的平方l2范数估计带宽。为了将约束变分优化问题转化为无约束形式，增加惩罚项和Lagrange乘子，以便快速收敛并增强约束，所以，最小化的目标函数转换为增广Lagrange函数：
Figure BDA0002971155300000093
式中，λ(t)为Lagrange乘子，α为惩罚系数，是数据真实性约束的平衡参数，
Figure BDA0002971155300000094
表示内积，cv(t)为本征模态分量，wv为cv(t)的中心频率。
VMD分解的过程简述如下：
(1)初始化：将
Figure BDA0002971155300000101
和z的初始值均设置为0，预定义收敛阈值，以及设置需要分离的本征模态分量的数量V；
(2)对于v＝1,2,…,V及所有w(w≥0)，更新每个本征模态分量cv(t)和中心频率wv，更新本征模态分量cv(t)的计算公式为：
Figure BDA0002971155300000102
更新中心频率wv的计算公式为：
Figure BDA0002971155300000103
(3)对于所有w≥0，更新Lagrange乘子；
(4)检查收敛条件：
Figure BDA0002971155300000104
其中，ε为预设的阈值。
如果条件满足，令
Figure BDA0002971155300000105
终止分解。否则，令z的值加1，返回步骤(2)。
VMD分解中最主要的参数有两个，一个是分解层数V，另一个是惩罚系数α。
在分解层数V的选择上，如果分解尺度太小，会导致VMD自适应的维纳滤波将重要信息滤除。反之若分解的IMF分量过多，会导致中心频段发生重叠现象。惩罚因子α的选取主要影响各IMF分量的带宽和收敛的速度。本实施例设定惩罚因子α的搜索区间为[20 3000]，搜索步长为20，分解导数V的搜索区间设[2 10]，搜索步长为1，对原始信号实施VMD分解，求出分解后的各IMF分量中最大峭度值并保存，从保存的峭度值中筛选出最大的，这个最大峭度值所采用的参数V值与α值即为最优。根据实验结果可知，对于心音信号中第一心音S1和第二心音S2而言，最优的分解层数V为6，惩罚系数α为2520。
S104：选择IMF分量：
信号的一个重要特性是它的振荡振幅，可以利用每个模态的振幅函数来计算每个IMF分量的概率密度分布函数值(PDF)。基于贝叶斯解释，PDF值表示有关系统或感兴趣信号的知识状态，而不仅仅是频率。因为PDF值包含有关其对应的IMF分量的完整信息，所以可以使用PDF分量的相似性度量来识别、捕捉包含更多有用信息的IMF。对于信号进行VMD分解的结果，可以计算得到各IMF分量之间的PDF值，PDF值包括KL散度、JS散度、互相关系数以及峭度的PDF值，然后计算各IMF分量的PDF值与原信号的PDF值之间的L2范数，L2范数最大的IMF分量即为信号的主模式，比主模式频带要低的IMF内是低频干扰，主要包括运动伪迹中的低频运动信号部分、肺音部分以及有色噪声的红色噪声中的低频部分；比主模式频带要高的IMF内是高频干扰，主要是运动伪迹中听诊头与人体的相对位移所引入高频噪声，以及各色有色噪声。图2是心音信号的VMD分解结果示例图。如图2所示，不难看出，有效心音信号主要集中在第二层内，第一层主要是以运动伪迹为主的低频基线漂移，第二层以下的主要是以高斯白噪声为主的高频噪声干扰。
根据对心音信号的PDF值进行分析，峭度对于心音信号的区分度最高。峭度Kurtosis反映随机变量分布特性的数值统计量，是归一化4阶中心矩。峭度对冲击信号特别敏感，因此，在本发明中，对于每个信号片段Xd,i(n′)经过VMD分解得到的V个分量
Figure BDA0002971155300000111
分别计算信号片段Xd,i(n′)和每个分量
Figure BDA0002971155300000112
的峭度值，选择峭度值与信号片段Xd,i(n′)峭度值的L2范数最大的分量作为信号片段Xd,i(n′)对应的IMF分量信号
Figure BDA0002971155300000113
S105：GHM多小波分解：
在步骤S104得到Nd个信号片段的IMF分量信号
Figure BDA0002971155300000114
后，对每个IMF分量信号
Figure BDA0002971155300000115
分别进行GHM多小波分解。由于多小波分解需要信号的长度为2λ的整数倍，λ表示GHM多小波系统的分解层数。因此需要先将各个IMF分量信号
Figure BDA0002971155300000116
插值到预设长度W，得到信号xd,i(n″)，n″＝1,2,…,W，W为2λ的整数倍值。本实施例中W为2λ和信号片段长度w的最小公倍数。
多小波系统对应的是一个多滤波器，滤波器系数是矩阵，输入具有向量形式，在多小波变换的情况下，需要对每个信号xd,i(n″)进行预滤波处理，得到大小为2×W的预滤波信号
Figure BDA0002971155300000121
本实施例中采用重复行滤波进行预滤波。重复行预滤波是通过对原始信号数据进行复制以得到第二行数据，这样虽然会造成数据的冗余，但是在去噪上能取得不错的效果。采用此方法得到的预滤波信号
Figure BDA0002971155300000122
可以表示为如下公式：
Figure BDA0002971155300000123
然后采用GHM多小波系统对每个预滤波信号
Figure BDA0002971155300000124
进行分解，记分解得到的信号分量数量为K，记第j个信号分量为yd,i,j，j＝1,2,…,K。由于GHM多小波系统中每层分解得到信号分量长度均为被分解信号的一半，那么可知信号分量yd,i,j的大小为2×W′，W′＝W/2λ，λ表示GHM多小波系统的分解层数。
GHM多小波系统对应的是一个多滤波器，其中低通滤波器由四个2×2的矩阵Ha给出，而高通滤波器由四个2×2矩阵Ga给出，a＝0,1,2,3。进行第一层分解后的公式可以表示为：
Figure BDA0002971155300000125
其中：
Figure BDA0002971155300000126
Figure BDA0002971155300000127
本实施例中四个2×2矩阵Ha和四个2×2矩阵Ga的具体取值分别如下：
Figure BDA0002971155300000128
Figure BDA0002971155300000131
第一层分解完成后，
Figure BDA0002971155300000132
表示高频分量，
Figure BDA0002971155300000133
表示低频分量，均为2×(W′/2)的矩阵。
在实际应用中，可以选择对高频分量或低频分量进行连续分解，也可以二者一起进行连续分解，即为小波包分解，本实施例中即采用小波包分解。图3是小波包分解示意图。如图3所示，对每一层分解得到的高频分量和低频分量均进行连续分解，经过λ层分解后，得到2λ个大小为2×W′的信号分量，即得到信号分量数量K＝2λ。
S106：构建时频域矩阵：
对于各个信号xd,i(n″)经GHM多小波分解得到K个大小为2×W′信号分量yd,i,j，构建得到大小为2K×W′的信号矩阵Zd,i：
Figure BDA0002971155300000134
对于同一心音信号样本的Nd个信号片段xd,i(n′)，采用其所对应的Nd个信号矩阵Zd,i构建得到大小为2K×Ld的时频域矩阵Md：
Figure BDA0002971155300000135
其中，
Figure BDA0002971155300000136
Figure BDA0002971155300000137
表示时频域矩阵Md中坐标为(α,β)的元素值，其中α＝1,2,…2K，β＝1,2,…,Ld。令
Figure BDA0002971155300000138
元素
Figure BDA0002971155300000139
按照以下公式确定：
Figure BDA00029711553000001310
其中，
Figure BDA0002971155300000141
分别表示向下取整和向上取整。
根据以上公式可知，在构建时频域矩阵Md时，相邻两个信号矩阵Zd,i各取一半元素进行叠加。
假设采样信号长度Sd＝3960，划分的信号片段数量Nd＝43，长度w＝90，信号片段的IMF分量信号插值后的长度W＝1440，GHM多小波系统的分解层数λ＝4，则W′＝W/2λ＝90。GHM多小波系统采用小波包分解，则时频域矩阵Md的大小为32×1980。图4是本实施例中时频域矩阵的示例图。
S107：时频域矩阵降维：
为了便于后续分析，本发明对于得到的每个时频域矩阵Md对行进行降维，得到降维后大小为P×Ld的时频域矩阵
Figure BDA0002971155300000142
其中P根据实际需要设置。
本实施例中采用主成分分析方法进行降维，主成分分析方法是一种常用的降维方法，其具体过程在此不再赘述。
S108：基于奇异值分解提取心音信号特征：
由于步骤S101中每个心音信号样本所标记的心音区间
Figure BDA0002971155300000143
是基于采样点数为Nd确定的，而时频域矩阵
Figure BDA0002971155300000144
的列数为Ld，因此需要对每个心音信号样本所标记的心音区间
Figure BDA0002971155300000145
进行展宽得到区间
Figure BDA0002971155300000146
其中γd＝[Ld/Nd]，[]表示取整，然后根据区间
Figure BDA0002971155300000147
从时频域矩阵
Figure BDA0002971155300000148
提取对应列构成大小为P×Ud的子矩阵
Figure BDA0002971155300000149
其中
Figure BDA00029711553000001410
对各个子矩阵
Figure BDA00029711553000001411
进行奇异值分解，其分解结果采用如下公式表示：
Figure BDA00029711553000001412
其中，Qd,k表示P×P阶的酉矩阵，Vd,k表示Ud×Ud阶的酉矩阵，∑d,k表示大小为P×Ud的奇异值矩阵。
在奇异值矩阵∑d,k选择最大的前M个奇异值，其中M根据实际需要设置，从矩阵Qd,k中提取出M个坐标为(m,m)的元素值
Figure BDA00029711553000001413
从矩阵Vd,k中提取出M个坐标为(m,m)的元素值
Figure BDA00029711553000001414
其中m＝1,2,…,M，然后拼接得到时间频率列向量Fd,k：
Figure BDA00029711553000001415
S109：训练分类模型：
将步骤S108得到的各个时间频率列向量Fd,k作为输入，对应的标签flagd,k作为输出，对预设的分类模型进行训练。本实施例中分类模型采用SVM(Support VectorMachines，支撑向量机)模型。
S110：获取待分割定位的心音信号：
采集待分割定位的心音信号
Figure BDA0002971155300000151
然后按照预设采样频率Fs进行采样得到采样信号
Figure BDA0002971155300000152
其中
Figure BDA0002971155300000153
N表示采样信号
Figure BDA0002971155300000154
所包含的采样点数。
S111：获取待分割定位心音信号的时频域矩阵：
采用步骤S102中的相同方法对待分割定位心音信号的采样信号
Figure BDA0002971155300000155
进行片段划分，得到N个采样信号片段
Figure BDA0002971155300000156
采用步骤S103中的相同方法对每个信号片段
Figure BDA0002971155300000157
进行VMD分解，采用步骤S104中的相同方法筛选得到IMF分量信号
Figure BDA0002971155300000158
采用步骤S105中的相同方法对每个IMF分量信号
Figure BDA0002971155300000159
进行插值、预滤波处理和GHM多小波分解，得到K个大小为2×W′的信号分量y′i,j，采用步骤S106中的相同方法构建得到大小为2K×L的时频域矩阵M′，
Figure BDA00029711553000001510
S112：提取香农包络：
对于待分割定位心音信号的时频域矩阵M′，分别计算每一列的香农能量SSE(b)，从而得到香农能量包络。香农能量SSE(b)的计算公式如下：
Figure BDA00029711553000001511
其中，M′(p,b)表示时频域矩阵M′中第p行第b列的元素值，b＝1,2,…,L。
接下来还可以对香农能量包络采用均值滤波器进行平滑处理，以实现去噪处理。
S113：心音信号边界估计：
对于待分割定位心音信号所对应的香农能量包络，筛选出香农能量最大值SSEmax，然后基于最大类间方差法获取区分噪音和心音的分类阈值TH。基于最大类间方差法获取区分噪音和心音的分类阈值TH的具体方法为：
1)初始化r＝1，初始化噪音集合
Figure BDA00029711553000001512
心音集合
Figure BDA00029711553000001513
2)对于每一列的香农能量SSE(b)，判断是否
Figure BDA00029711553000001514
R为根据需要设置的常数，本实施例中R＝256，如果是，则将香农能量SSE(b)加入噪音集合φ0，否则将香农能量SSE(b)加入心音集合φ1。
3)根据以下公式计算当前分类的类间方差gr：
gr＝ωr0×ωr1×(μr0-μr1)2
其中，ωr0＝|φ0|/L，ωr1＝|φ1|/L，|φ0|和|φ1|分别表示噪音集合φ0和心音集合φ1中的香农能量数量，μr0和μr1分别表示噪音集合φ0和心音集合φ1中的香农能量平均值。
4)判断是否r＜R，如果是，进入步骤5)，否则进入步骤6)。
5)令r＝r+1，返回步骤2)。
6)筛选出类间方差gr中的最大值gr*，则分类阈值
Figure BDA0002971155300000161
计算得到阈值之后使用单参数的双门限法来定位到心音的边界，具体方法为：使用类间方差法确定的分类阈值TH作为较高的门限T1，即T1＝TH，先进行一次粗判，具体方法为：在香农能量包络中提取出连续高于阈值TH的区间，将区间长度小于预设阈值的区间删除，剩余区间即为第一心音S1或第二心音S2所在的粗定心音区间
Figure BDA0002971155300000162
K′表示待分割定位心音信号的心音区间数量。相当于高于阈值TH的肯定是有效心音，而心音的实际起止点还位于TH和SSE包络交点所对应的时间点之外。然后再在TH的基础上确定一个较低的阈值T2，令T2＝TH-η，η为预设的大于0的偏差值，本实施例中η＝0.15。在Z′k′_temp的基础上向两旁扩展搜索，分别找到香农能量包络与阈值T2相交的两个端点
Figure BDA0002971155300000163
作为心音的正确起始端点，得到心音区间
Figure BDA0002971155300000164
显然，心音区间的起点和终点对应时频域矩阵M′的列序号。
S114：提取待分割定位心音信号的特征：
根据步骤S112得到的心音区间Zk′，从时频域矩阵M′提取对应列构成子矩阵M′k′，采用步骤S108中的相同方法，基于奇异值分解得到该心音区间的时间频率列向量F′k′。
S115：心音信号分类：
将每个心音区间的时间频率列向量F′k′分别输入步骤S109训练好的分类模型，其分类结果即为该心音区间内心音信号的类型。
S116：心音信号分割定位：
由于待分割定位心音信号的时频域矩阵M′的列数是其采样信号x′(n)采样点数的2K倍，为了将心音区间还原到原始信号中，需要先进行心音区间的尺度变换，即将待分割定位心音信号的每个心音区间
Figure BDA0002971155300000171
还原得到区间
Figure BDA0002971155300000172
γ＝[L/N]，然后在采样信号x′(n)标记每个区间Zk′，并根据步骤S115中的分类结果确定该区间内心音信号的类型，从而得到心音信号的分割定位结果。
为了更好地说明本发明的技术效果，采用一个具体实例对本发明进行实验验证。本次实验验证中的心音数据集来自三部分：
数据集1来自于“Classifying Heart Sounds Pascal Challenge Competition”，该数据集包括数据集A和数据集B两部分。其中数据集A为通过iStethoscope Pro iPhone应用程序从公众中提取的数据。数据集B为在巴西累西腓的Real Portugu es医院(RHP)使用数字示波器收集系统采集的孕产妇和胎儿心脏病病房收集的200多个听诊病例组成(Pereira等，2011)，不具备一般性，因此只使用了该数据集得数据集A部分。该数据集得数据采样数据长度1s-30s不等，采样频率为44100Hz，采样位数为16位，本次实验验证中选取其中31例正常心音信号和34例含杂音的心音信号。
数据集2来自于2016年PhysioNet/CinC的心音挑战赛。该数据集的心音数据来自于9个独立的数据库，除去SUFHSDB数据库的数据(来自于胎儿和母亲的心音)，建立了8个独立的心音数据库。其中四个数据库以7比3的比例，被分为训练集与测试集，其余的四个数据库被专门分配为训练集或测试集。该数据集共包括从764例心脏录音/病人身上录取的3153段心音信号，持续时间从5s到120s不等，采样频率为2000Hz。本次实验验证中主要选取了其中超过7S的116个正常心音和290个异常心音进行了处理。
数据集3来自于采用一款可穿戴心音监测平台自行采集的心音信号，图5是本实施例中自行采集心音信号示例图。根据图5可以看到，可穿戴心音监测平台所采集到的信号底噪更大，且容易发生S1和S2的混叠。本数据集共包括20段心音信号，持续时间都被裁剪为4s，采样频率为2000Hz。
本实施例主要从抗噪鲁棒性和抗噪上限两个方面来进行实验验证。
·实验一：抗噪鲁棒性实验
首先对本发明的抗噪鲁棒性进行实验验证。噪声可分为白噪声、粉噪声、蓝噪声、肺音和运动伪迹五种。白噪声在全频段均匀分布；粉噪音为白噪声将高频声音的响度降低，将低频声音的响度提高；蓝色噪声则为低频声音变得更弱，高频声音变得更强；肺音属于人体的生理音，这类噪声的引入是由心音听诊这种采集方式所决定的，现实中听诊获得的心音，或者肺音信号是互为干扰，同时存在的，互相影响，所以听诊得到的心肺音信号在时域上是混叠的，这对听诊的准确性造成很大干扰；运动伪迹是由可穿戴心音监测装置与人体和皮肤之间的相对位移而引起的，噪音信号幅度大，且频谱与心音信号频谱部分重叠，难以滤波消除。为了体现本发明在抗噪鲁棒性方面的优势，采用2种算法作为对比方法，分别为VMD+EWT(经验小波)算法和EMD+PCA(主成分分析)算法，采用信噪比(SNR)、相关系数和均方误差(MSE)作为对于噪声的鲁棒性评估指标。
图6是本实施例中纯净心音信号叠加白噪声所得到的归一化心音信号波形图。图7是VMD+EWT算法采用对图6所示心音信号进行处理后的归一化信号波形图。图8是采用EMD+PCA算法对图5所示心音信号进行处理后的归一化信号波形图。图9是采用本发明中多小波+PCA算法对图6所示心音信号进行处理后的归一化信号波形图。对比图6至图9可知，由于所叠加白噪音强度较大，较为纯净的心音信号已经几乎湮没在了高斯白噪声中，VMD+EWT算法的去噪效果比较差，与原信号相比区别不大，EMD+PCA算法甚至抑制了部分的有效信号并不能肉眼分辨S1和S2的位置，而本发明多小波+PCA算法处理之后，第一心音S1和第二心音S2的位置已经相对明显了。
图10是本发明和两种对比方法在不同信噪比高斯白噪声下心音信号去噪后的信噪比比较图。图11是本发明和两种对比方法在不同信噪比高斯白噪声下心音信号去噪后的相关系数比较图。图12是本发明和两种对比方法在不同信噪比高斯白噪声下心音信号去噪后的均方误差比较图。如图10到图12可知，随着噪声信噪比SNR的增大，三种方法一般符合以下规律：对源心音信号与去噪后心音信号的相关系数和均方误差的提升逐步增大。其中本发明和VMD+EWT算法的效果总体要优于EMD+PCA算法的效果，在噪声较高的时候，本发明和VMD+EWT算法的效果差距不算很大，但是随着噪声的增多，本发明处理的结果就开始展现出他的优势，而且在全过程都能获得比较稳定的增益，说明本发明对高斯白噪声的变化更敏感，能获得更好的去噪效果。
图13是本发明和两种对比方法在不同信噪比粉噪声下心音信号去噪后的信噪比比较图。图14是本发明和两种对比方法在不同信噪粉噪声下心音信号去噪后的相关系数比较图。图15是本发明和两种对比方法在不同信噪比粉噪声下心音信号去噪后的均方误差比较图。如图13到图15可知，EMD+PCA算法对粉噪声并不敏感，处理效果比较差，而本发明方法处理的效果要优于VMD+EWT算法。
图16是本发明和两种对比方法在不同信噪比蓝噪声下心音信号去噪后的信噪比比较图。图17是本发明和两种对比方法在不同信噪比蓝噪声下心音信号去噪后的相关系数比较图。图18是本发明和两种对比方法在不同信噪比蓝噪声下心音信号去噪后的均方误差比较图。如图16到图18可知，添加了较多蓝噪声的情况下，在低噪的情况下，VMD+EWT算法的性能略优一点，但是本发明整体处理效果最好。
图19是本发明和两种对比方法在不同信噪比肺音下心音信号去噪后的信噪比比较图。图20是本发明和两种对比方法在不同信噪比肺音下心音信号去噪后的相关系数比较图。图21是本发明和两种对比方法在不同信噪比肺音下心音信号去噪后的均方误差比较图。如图19到图21可知，添加了肺音的情况下，本发明整体处理效果优于两种对比方法。
图22是本发明和两种对比方法在不同信噪比运动伪迹下心音信号去噪后的信噪比比较图。图23是本发明和两种对比方法在不同信噪比运动伪迹下心音信号去噪后的相关系数比较图。图24是本发明和两种对比方法在不同信噪比运动伪迹下心音信号去噪后的均方误差比较图。如图22到图24可知，添加了运动伪迹的情况下，本发明处理效果与VMD+EWT算法的处理效果相差不大，均优于EMD+PCA算法的处理效果。
综上所述，本发明在不同噪声条件下对心音信号的去噪效果都较好。
·实验二：抗噪上限实验
在抗噪上限实验中，除了VMD+EWT算法和EMD+PCA算法外，还采用了S变换+香农能量包络算法和EWT+香农熵算法作为对比方法。表1是本发明和四种对比方法对于对单种噪声的抗噪上限数据表。
Figure BDA0002971155300000201
表1
从表1可以看出，本发明所提出的方法要明显优于其他几种对比方法，经测试，本发明方法能在各种噪声环境下实现带噪为-8dB的信号中稳定地分割心音信号。
图25是采用本发明对图6所示心音信号进行处理后提取的香农包络结果图。图26是采用本发明对图6所示心音信号进行处理后的处理结果图。可以看到，在添加了较大噪声，信号已经严重失真的情况下本发明依旧能对心音信号进行稳定而准确的分割。

Claims (7)

1.一种心音分割定位方法，其特征在于，包括以下步骤：
S1：采集若干心音信号样本Xd(t)，其中d＝1,2,…,D，D表示心音信号样本的数量，t表示时间；对每个心音信号样本Xd(t)进行去噪处理，然后按照采样频率Fs进行重采样得到每个心音信号样本的采样信号Xd(n)，其中n＝1,2,…,Sd，Sd表示采样信号Xd(n)的采样点数；在每个心音信号样本的采样信号Xd(n)中标记第一心音S1和第二心音S2的区间，记标记得到的第k个心音区间为
Figure FDA0002971155290000011
其中
Figure FDA0002971155290000012
分别表示第d个心音信号样本的采样信号Xd(n)中第k个心音区间的起点和终点，k＝1,2,…,Kd，Kd表示采样信号xd(n)中心音区间的数量，并根据心音类型设置其标签flagd,k，flagd,k＝1,2；
S2：采用滑动窗对各个采样信号Xd(n)进行片段划分得到Nd个采样信号片段Xd,i(n′)，i＝1,2,…,Nd，n′＝1,2,…,w，w表示采样信号Xd(n)每个信号片段的长度；
S3：对于每个采样信号Xd(n)，在步骤S2划分得到Nd个信号片段Xd,i(n′)后，对每个信号片段分别进行V层VMD分解，记信号片段Xd,i(n′)经VMD分解得到的第v个分量为
Figure FDA0002971155290000013
S4：对于每个信号片段Xd,i(n′)经过VMD分解得到的V个分量
Figure FDA0002971155290000014
分别计算信号片段Xd,i(n′)和每个分量
Figure FDA0002971155290000015
的峭度值，选择峭度值与信号片段Xd,i(n′)峭度值的L2范数最大的分量作为信号片段Xd,i(n′)对应的IMF分量信号
Figure FDA0002971155290000016
S5：将各个IMF分量信号
Figure FDA0002971155290000017
插值到预设长度W，得到信号xd,i(n″)，n″＝1,2,…,W，W为2λ的整数倍值；然后对每个信号片段xd,i(n″)进行预滤波处理，得到大小为2×W的预滤波信号
Figure FDA0002971155290000018
采用GHM多小波系统对每个预滤波信号
Figure FDA0002971155290000019
进行分解，记分解得到的信号分量数量为K，记第j个信号分量为yd,i,j，j＝1,2,…,K，信号分量yd,i,j的大小为2×W′，W′＝W/2λ，λ表示GHM多小波系统的分解层数；
S6：对于各个信号xd,i(n″)经GHM多小波分解得到K个大小为2×W′信号分量yd,i,j，构建得到大小为2K×W′的信号矩阵Zd,i：
Figure FDA00029711552900000110
对于同一心音信号样本的Nd个信号片段xd,i(n′)，采用其所对应的Nd个信号矩阵Zd,i构建得到大小为2K×Ld的时频域矩阵Md：
Figure FDA0002971155290000021
其中，
Figure FDA0002971155290000022
Figure FDA0002971155290000023
表示时频域矩阵Md中坐标为(α,β)的元素值，其中α＝1,2,…2K，β＝1,2,…,Ld；令
Figure FDA0002971155290000024
元素
Figure FDA0002971155290000025
按照以下公式确定：
Figure FDA0002971155290000026
其中，
Figure FDA0002971155290000027
分别表示向下取整和向上取整；
S7：对于得到的每个时频域矩阵Md对行进行降维，得到降维后大小为P×Ld的时频域矩阵
Figure FDA0002971155290000028
其中P根据实际需要设置；
S8：对步骤S1中每个心音信号样本所标记的心音区间
Figure FDA0002971155290000029
进行展宽得到区间
Figure FDA00029711552900000210
其中γd＝[Ld/Nd]，然后根据区间
Figure FDA00029711552900000211
从时频域矩阵
Figure FDA00029711552900000212
提取对应列构成大小为P×Ud的子矩阵
Figure FDA00029711552900000213
其中
Figure FDA00029711552900000214
对各个子矩阵
Figure FDA00029711552900000215
进行奇异值分解，其分解结果采用如下公式表示：
Figure FDA00029711552900000216
其中，Qd,k表示P×P阶的酉矩阵，Vd,k表示U×U阶的酉矩阵，∑d,k表示大小为P×Ud的奇异值矩阵；
在奇异值矩阵∑d,k选择最大的前M个奇异值，其中M根据实际需要设置，从矩阵Qd,k中提取出M个坐标为(m,m)的元素值
Figure FDA00029711552900000217
从矩阵Vd,k中提取出M个坐标为(m,m)的元素值
Figure FDA00029711552900000218
其中m＝1,2,…,M，然后拼接得到时间频率列向量Fd,k：
Figure FDA00029711552900000219
S9：将步骤S8得到的各个时间频率列向量Fd,k作为输入，对应的标签flagd,k作为输出，对预设的分类模型进行训练；
S10：采集待分割定位的心音信号
Figure FDA00029711552900000220
然后按照预设采样频率Fs进行采样得到采样信号
Figure FDA0002971155290000031
其中
Figure FDA0002971155290000032
N表示采样信号
Figure FDA0002971155290000033
所包含的采样点数；
S11：采用步骤S2中的相同方法对待分割定位心音信号的采样信号
Figure FDA0002971155290000034
进行片段划分，得到N个采样信号片段
Figure FDA0002971155290000035
采用步骤S3中的相同方法对每个信号片段
Figure FDA0002971155290000036
进行VMD分解，采用步骤S4中的相同方法筛选得到IMF分量信号
Figure FDA0002971155290000037
采用步骤S5中的相同方法对每个IMF分量信号
Figure FDA0002971155290000038
进行插值、预滤波处理和GHM多小波分解，得到K个大小为2×W′的信号分量y′i,j，采用步骤S6中的相同方法构建得到大小为2K×L的时频域矩阵M′，
Figure FDA0002971155290000039
S12：对于待分割定位心音信号的时频域矩阵M′，分别计算每一列的香农能量SSE(b)，从而得到香农能量包络；
S13：对于待分割定位心音信号所对应的香农能量包络，筛选出香农能量最大值SSEmax，然后基于最大类间方差法获取区分噪音和心音的分类阈值TH，然后使用类间方差法确定的分类阈值TH作为较高的门限T1，即T1＝TH，在香农能量包络中提取出连续高于阈值TH的区间，将区间长度小于预设阈值的区间删除，剩余区间即为第一心音S1或第二心音S2所在的粗定心音区间
Figure FDA00029711552900000310
K′表示待分割定位心音信号的心音区间数量；然后再在TH的基础上确定一个较低的阈值T2，令T2＝TH-η，η为预设的大于0的偏差值，在Z′k′_temp的基础上向两旁扩展搜索，分别找到香农能量包络与阈值T2相交的两个端点
Figure FDA00029711552900000311
作为心音的正确起始端点，得到心音区间
Figure FDA00029711552900000312
S14：根据步骤S12得到的心音区间Zk′，从时频域矩阵M′提取对应列构成子矩阵M′k′，采用步骤S8中的相同方法，基于奇异值分解得到该心音区间的时间频率列向量F′k′；
S15：将每个心音区间的时间频率列向量F′k′分别输入步骤S9训练好的分类模型，其分类结果即为该心音区间内心音信号的类型；
S16：将待分割定位心音信号的每个心音区间
Figure FDA00029711552900000313
还原得到区间
Figure FDA00029711552900000314
γ＝[L/N]，然后在采样信号x′(n)标记每个区间Zk′，并根据步骤S15中的分类结果确定该区间内心音信号的类型，从而得到心音信号的分割定位结果。
2.根据权利要求1所述的心音分割定位方法，其特征在于，所述步骤S2中滑动窗采用高斯窗。
3.根据权利要求1所述的心音分割定位方法，其特征在于，所述步骤S3中VMD分解时所采用的分解层数V为6，惩罚系数α为2520。
4.根据权利要求1所述的心音分割定位方法，其特征在于，所述步骤S5中长度W为2λ和信号片段长度w的最小公倍数。
5.根据权利要求1所述的心音分割定位方法，其特征在于，所述步骤S5中预滤波采用重复行滤波，所得到的预滤波信号
Figure FDA0002971155290000041
的表达式为：
Figure FDA0002971155290000042
6.根据权利要求1所述的心音分割定位方法，其特征在于，所述步骤S5中在采用GHM多小波系统对每个预滤波信号
Figure FDA0002971155290000043
进行分解时，采用小波包分解，对每一层分解得到的高频分量和低频分量均进行连续分解，经过λ层分解后，得到信号分量数量K＝2λ。
7.根据权利要求1所述的心音分割定位方法，其特征在于，所述步骤S12中还包括对香农能量包络采用均值滤波器进行平滑处理。
CN202110263717.5A 2021-03-11 2021-03-11 基于vmd和多小波的心音分割定位方法 Active CN113066502B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110263717.5A CN113066502B (zh) 2021-03-11 2021-03-11 基于vmd和多小波的心音分割定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110263717.5A CN113066502B (zh) 2021-03-11 2021-03-11 基于vmd和多小波的心音分割定位方法

Publications (2)

Publication Number Publication Date
CN113066502A true CN113066502A (zh) 2021-07-02
CN113066502B CN113066502B (zh) 2022-07-26

Family

ID=76559968

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110263717.5A Active CN113066502B (zh) 2021-03-11 2021-03-11 基于vmd和多小波的心音分割定位方法

Country Status (1)

Country Link
CN (1) CN113066502B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113705420A (zh) * 2021-08-24 2021-11-26 青岛理工大学 一种风电支撑结构模态响应分量提取方法及系统
CN113749620A (zh) * 2021-09-27 2021-12-07 广州医科大学附属第一医院（广州呼吸中心） 一种睡眠呼吸暂停检测方法、系统、设备及存储介质
CN114869294A (zh) * 2022-05-05 2022-08-09 电子科技大学 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法
CN116631429A (zh) * 2023-07-25 2023-08-22 临沂金诺视讯数码科技有限公司 基于voip呼叫volte通话的语音视频处理方法及系统
CN117351988A (zh) * 2023-12-06 2024-01-05 方图智能(深圳)科技集团股份有限公司 一种基于数据分析的远程音频信息处理方法及系统

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100249629A1 (en) * 2009-03-18 2010-09-30 Coloplast A/S Segmenting a cardiac acoustic signal
CN103948398A (zh) * 2014-04-04 2014-07-30 杭州电子科技大学 适用于Android系统的心音定位分段方法
CN108682433A (zh) * 2018-06-01 2018-10-19 四川长虹电器股份有限公司 基于mfcc的一阶差分系数的心音类型识别方法
CN108814642A (zh) * 2018-05-16 2018-11-16 合肥康聆医疗科技有限公司 一种电子听诊器的心音定位及心率计算方法
CN109646042A (zh) * 2019-01-29 2019-04-19 电子科技大学 一种基于压电传感器的可穿戴心音和肺音监测装置
EP3608918A1 (en) * 2018-08-08 2020-02-12 Tata Consultancy Services Limited Parallel implementation of deep neural networks for classifying heart sound signals
CN110916716A (zh) * 2019-12-30 2020-03-27 龙岩学院 一种可穿戴式心音监测设备
CN111317499A (zh) * 2018-12-17 2020-06-23 天津光电通信技术有限公司 一种基于小波技术的心音信号处理方法
US20200205771A1 (en) * 2015-06-15 2020-07-02 The Research Foundation For The State University Of New York System and method for infrasonic cardiac monitoring
US10729336B1 (en) * 2006-06-30 2020-08-04 Bao Tran Smart watch
CN111528900A (zh) * 2020-05-21 2020-08-14 广东工业大学 基于巴特沃斯滤波器与香农熵法的心音分段方法和装置

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10729336B1 (en) * 2006-06-30 2020-08-04 Bao Tran Smart watch
US20100249629A1 (en) * 2009-03-18 2010-09-30 Coloplast A/S Segmenting a cardiac acoustic signal
CN103948398A (zh) * 2014-04-04 2014-07-30 杭州电子科技大学 适用于Android系统的心音定位分段方法
US20200205771A1 (en) * 2015-06-15 2020-07-02 The Research Foundation For The State University Of New York System and method for infrasonic cardiac monitoring
CN108814642A (zh) * 2018-05-16 2018-11-16 合肥康聆医疗科技有限公司 一种电子听诊器的心音定位及心率计算方法
CN108682433A (zh) * 2018-06-01 2018-10-19 四川长虹电器股份有限公司 基于mfcc的一阶差分系数的心音类型识别方法
EP3608918A1 (en) * 2018-08-08 2020-02-12 Tata Consultancy Services Limited Parallel implementation of deep neural networks for classifying heart sound signals
CN111317499A (zh) * 2018-12-17 2020-06-23 天津光电通信技术有限公司 一种基于小波技术的心音信号处理方法
CN109646042A (zh) * 2019-01-29 2019-04-19 电子科技大学 一种基于压电传感器的可穿戴心音和肺音监测装置
CN110916716A (zh) * 2019-12-30 2020-03-27 龙岩学院 一种可穿戴式心音监测设备
CN111528900A (zh) * 2020-05-21 2020-08-14 广东工业大学 基于巴特沃斯滤波器与香农熵法的心音分段方法和装置

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
AMIR DADASHNIALEHI ET AL: "Deep Learning for Texture Classification Via Multi-wavelet Fusion of Scattering Transforms", 《2017 IEEE INTERNATIONAL CONFERENCE ON MECHATRONICS》 *
MADHUSUDHAN MISHRA ET AL: "Detection of Third Heart Sound Using Variational Mode Decomposition", 《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》 *
SHIJI XIAHOU ET AL: "A strong anti-noise segmentation algorithm based on variational mode decomposition and multi-wavelet for wearable heart sound acquisition system", 《REVIEW OF SCIENTIFIC INSTRUMENTS》 *
周宁: "基于电子听诊器的心音定位及心肺音分离方法研究", 《中国优秀博硕士学位论文全文数据库(硕士)医药卫生科技辑》 *
梁宇航: "面向可穿戴应用的心音信号处理方法研究", 《中国优秀博硕士学位论文全文数据库(硕士)医药卫生科技辑》 *
郭萌: "心音信号自适应阈值小波去噪及先心病诊断方法研究", 《中国优秀博硕士学位论文全文数据库(硕士)医药卫生科技辑》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113705420A (zh) * 2021-08-24 2021-11-26 青岛理工大学 一种风电支撑结构模态响应分量提取方法及系统
CN113749620A (zh) * 2021-09-27 2021-12-07 广州医科大学附属第一医院（广州呼吸中心） 一种睡眠呼吸暂停检测方法、系统、设备及存储介质
CN113749620B (zh) * 2021-09-27 2024-03-12 广州医科大学附属第一医院（广州呼吸中心） 一种睡眠呼吸暂停检测方法、系统、设备及存储介质
CN114869294A (zh) * 2022-05-05 2022-08-09 电子科技大学 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法
CN114869294B (zh) * 2022-05-05 2023-05-30 电子科技大学 基于vmd分解和let模型的粒子滤波运动伪迹抑制方法
CN116631429A (zh) * 2023-07-25 2023-08-22 临沂金诺视讯数码科技有限公司 基于voip呼叫volte通话的语音视频处理方法及系统
CN116631429B (zh) * 2023-07-25 2023-10-10 临沂金诺视讯数码科技有限公司 基于voip呼叫volte通话的语音视频处理方法及系统
CN117351988A (zh) * 2023-12-06 2024-01-05 方图智能(深圳)科技集团股份有限公司 一种基于数据分析的远程音频信息处理方法及系统
CN117351988B (zh) * 2023-12-06 2024-02-13 方图智能(深圳)科技集团股份有限公司 一种基于数据分析的远程音频信息处理方法及系统

Also Published As

Publication number Publication date
CN113066502B (zh) 2022-07-26

Similar Documents

Publication Publication Date Title
CN113066502B (zh) 基于vmd和多小波的心音分割定位方法
CN111107785B (zh) 使用短的单导联ecg记录来检测心房颤动
Lin et al. Heartbeat classification using normalized RR intervals and morphological features
CN108714026B (zh) 基于深度卷积神经网络和在线决策融合的细粒度心电信号分类方法
Safieddine et al. Removal of muscle artifact from EEG data: comparison between stochastic (ICA and CCA) and deterministic (EMD and wavelet-based) approaches
CN111358455B (zh) 一种多数据源的血压预测方法和装置
CN111310570B (zh) 一种基于vmd和wpd的脑电信号情感识别方法及系统
CN108742697B (zh) 心音信号分类方法及终端设备
CN110432895B (zh) 训练数据处理、心电波形检测方法及电子设备
Khan et al. Separating Heart Sound from Lung Sound UsingLabVIEW
CN112971839A (zh) 一种基于前馈卷积神经网络的心音分类方法
CN113723557A (zh) 一种基于多频带时空卷积网络的抑郁脑电分类系统
TWI620546B (zh) 腦波分析方法及其裝置
Thomas et al. Heart sound segmentation using fractal decomposition
CN117357080A (zh) 近红外光谱信号去噪方法及装置、终端设备、存储介质
CN117017297A (zh) 驾驶员疲劳的预测和识别模型建立方法及其应用
CN116211322A (zh) 一种基于机器学习脑电信号的抑郁症识别方法和系统
CN115017960B (zh) 一种基于时空联合mlp网络的脑电信号分类方法及应用
CN113229842B (zh) 一种基于复数深度神经网络的心肺音自动分离方法
Lu et al. Dual temporal convolutional network for single-lead fibrillation waveform extraction
CN112784686A (zh) 一种基于连续多变量变分模态分解的自适应信号分析方法
NSVN et al. Optimal threshold estimation using cultural algorithm for EMD-DWT based ECG denoising
AL-Jobouri et al. From biomedical signal processing techniques to fMRI parcellation
Chieng et al. Qualitative and quantitative performance comparison of ECG noise reduction and signal enhancement method based on various digital filter designs and discrete wavelet transform
Meltzer et al. A clustering approach to construct multi-scale overcomplete dictionaries for ECG modeling

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
GR01 Patent grant
GR01 Patent grant