ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

心音信号特征提取分析系统完整实现:从预处理到分类识别

心音信号特征提取分析系统完整实现:从预处理到分类识别 简介心音信号作为生物医学工程中重要的生理参数其分析涉及信号处理、特征工程与机器学习等基础技术。在心音采集过程中信号常混有呼吸噪声与环境干扰预处理环节的质量直接影响后续分析效果因此带通滤波与采样率统一是首要步骤。通过包络检测与峰值定位可提取S1、S2心音成分进而从时域、频域及小波包分解中构造多维特征如MFCC、能量比等。针对高维特征易引发的过拟合问题可结合mRMR筛选与PCA降维优化特征空间再采用SVM或随机森林完成分类识别。这些技术广泛应用于医疗辅助诊断、可穿戴设备及远程健康监测场景。本文系统梳理了心音信号特征提取分析系统的完整技术链路涵盖预处理、分割、特征构建、特征选择、模型验证及常见问题排查为工程实践者提供可直接落地的参考方案。 最近好多做生物医学工程毕设的学弟学妹跑来问我说拿到一个“心音信号特征提取分析系统”的课题不知道从哪里下手。这个题目看着很常见年年都有人做但真正要做出一个能跑通、能交差、能答辩的系统远不是把几个函数拼在一起那么简单。它涉及到完整的信号处理链路、分割算法的选型、特征体系的构建还有后期分类验证这一整套流程。我索性把这个课题的完整实现思路、代码细节、踩坑记录整理成一篇长文送给正在被这个project折磨的朋友们。1. 系统整体设计与思路拆解1.1 先把这个系统的功能边界想清楚很多人拿到题目就急着写代码结果写了三天发现自己连“心音信号该用什么采样率”都没搞清楚。我先帮你把问题拉回到系统层面。这个系统本质上要回答三个问题输入是什么、中间要做什么、输出是什么。输入通常是来自公开数据库或自采设备的心音记录例如PhysioNet的CinC Challenge 2016数据集格式一般是.wav或者.mat。输出一般分两种一是对每条心音记录提取出能反映心脏生理状态的特征向量二是基于这些特征做正常/异常的二分类判断。中间的处理环节就是信号预处理、心音分割定位、特征提取、特征选择和分类评估。把这个功能边界画出来之后你会发现整个系统其实是一个标准的“信号处理机器学习”流水线。每一环节都有独立的研究点但完整做下来又相互影响。一个常见的误区是每步都想去创新结果每一步都很糙。一个合格的毕设系统不要求每个模块都是顶级算法但要求链路完整、逻辑自洽、效果可复现。1.2 为什么我依然推荐用MATLAB做这个课题我知道现在Python也很流行深度学习方向的论文基本都拥抱Python生态。但心音信号特征提取这个课题我依然首推MATLAB尤其对于本科生和没系统学过编程的医工学生。MATLAB在信号预处理上真的是几分钟出图。读音频、重采样、滤波、画频谱图都是一行命令的事情不像Python还要装soundfile、scipy、librosa一堆依赖。尤其是调试阶段你对着滤波前后的时域波形和功率谱图观察比对着终端里的数字直观得多。MATLAB的滤波设计工具箱里有filterDesigner这种图形化工具拉一拉参数就能看幅度响应对确定截止频率这种问题非常友好。另一个现实原因是毕设资源。网上的心音信号处理开源代码MATLAB版本占了相当大的比例。无论是参考还是复现MATLAB都更省力。当然你完全可以用Python重写但我见过太多同学把时间浪费在环境配置上结果核心算法没时间打磨这就本末倒置了。1.3 技术路线的完整流程整个系统的技术路线我个人会分成七个环节来推进数据读取与统一采样率解决不同来源的原始信号采样率不一致问题重采样到统一频率。信号预处理带通滤波去除呼吸音、环境噪声等干扰成分并进行归一化。心音分段定位利用包络检测算法找出每一个心动周期中的S1、S2心音位置。特征提取从时域、频域、时频域多个维度计算特征参数。特征筛选用mRMR等方法剔除冗余特征保留对分类最有用的维度。分类模型训练与验证选择SVM或随机森林等分类器用交叉验证评估性能。结果可视化与评估绘制ROC曲线、混淆矩阵等用于论文或汇报展示。这里每一环后面都有大量细节。下面我从预处理开始逐步把每一步的原理和代码都讲透。2. 心音信号的预处理决定成败的第一道关口2.1 数据读入与采样率统一很多人直接读.wav文件就开始滤波这是第一个坑。不同医院、不同设备采集的心音采样率可能完全不同常见的有2000Hz、4000Hz、8000Hz甚至11025Hz。你如果混在一起做统一的带通滤波截止频率设置就会出问题。比如你想滤掉25Hz以下的干扰如果某条信号采样率是2000HzNyquist频率就是1000Hz你设置的滤波器边界针对的是归一化频率。当不同采样率混在一起同样的归一化截止频率对应的实际频率就变了处理结果自然不可比。所以第一步永远是统一采样率。不用太高心音的主要能量集中在20到200Hz之间就算考虑高频杂音和瓣膜开闭音1000到2000Hz已经足够。我习惯统一重采样到1000Hz既保留有效频带又能减少后续计算量。代码很简单% 假设 rawSig 是原始心音信号fs 是原始采样率 targetFs 1000; sigResampled resample(rawSig, targetFs, fs);resample函数内部会自动做抗混叠滤波所以直接用就行。如果你有多条信号最好把所有信号统一成同一个采样率再存成统一的mat结构体方便后续批量处理。2.2 带通滤波的参数设计与零相位滤波器心音信号中最重要的是S1和S2两者的频率范围基本在20到150Hz之间。但舒张期杂音、收缩期喷射音等病理成分频率能到400Hz甚至更高所以我的经验是带通范围取25到400Hz比较稳妥。低于25Hz的基本是呼吸音、肌肉活动伪迹、基线漂移高于400Hz的多为环境噪声和传感器电子噪声。滤波器阶数选4阶就够了。太高阶不仅引入相位失真还会造成计算量浪费。这里有个关键细节滤波一定要用零相位滤波也就是filtfilt函数千万不要用filter。filter是有因果性的会带来相位延迟导致滤波后的波形在时间轴上发生偏移心音峰值的位置就不准了。心音分割靠的就是峰值位置位置错了整条检测链路就全错了。[b, a] butter(4, [25 400]/(targetFs/2), bandpass); sigFiltered filtfilt(b, a, sigResampled);filtfilt会先把信号正向滤波一遍再反向滤波一遍相位延迟刚好抵消输出信号的峰值位置和原始信号保持一致。代价是边界效应所以滤波前最好先截去开头和结尾各一小段。2.3 归一化与整段信号处理不同设备采集的增益不同信号幅值差异很大。有的人直接拿原始幅值去做包络和阈值检测结果同一个算法在不同记录上要么阈值太高检不出要么阈值太低误检一片。所以滤波之后一定要做归一化。我常用的归一化方式不是简单的除以最大值而是减去均值再除以标准差也就是z-score标准化。这样处理后信号的均值为0方差为1不管原始幅值范围是±0.1还是±10后续处理都在一个可比的尺度上。sigNorm (sigFiltered - mean(sigFiltered)) / std(sigFiltered);预处理到这一步信号基本可以进入分割阶段了。不过在做分割前建议先画一下滤波前后的频谱图对照确认没有把心音的有效频段滤掉。很多时候你以为自己在做降噪实际上滤波器参数设计错了把S1、S2的低频能量消掉大半后面再怎么调算法都救不回来。3. 心音分割与特征提取系统的核心发力点3.1 基于包络法的心音定位心音分割最经典的做法是先提取心音信号的包络然后对包络做峰值检测。常见的包络提取方式有绝对值包络、Hilbert包络、香农能量包络等。我强烈推荐香农能量包络它对低幅值心音比较敏感也能抑制高幅值噪声从而让真正的S1、S2峰值更容易显现出来。香农能量的计算公式是对每个采样点的能量先做归一化然后取对数再取负。简单理解就是把信号能量映射到一个对数量级上让微弱的心音成分不至于被强噪声完全淹没。实现时一般用滑动窗口窗口长度取20ms比较合适。窗口太长会抹平S1和S2之间的间隔窗口太短则包络毛刺太多影响峰值检测。windowLen round(0.02 * targetFs); se zeros(size(sigNorm)); halfWin floor(windowLen / 2); for i 1:length(sigNorm) idxLow max(1, i - halfWin); idxHigh min(length(sigNorm), i halfWin); seg sigNorm(idxLow:idxHigh); energy seg.^2; se(i) -mean(energy .* log(energy eps)); end % 对包络做平滑 seSmooth movmean(se, round(0.01 * targetFs));得到平滑后的包络曲线后峰值检测可以用MATLAB的findpeaks函数。需要调两个参数一是MinPeakDistance这个值根据心率范围估算。正常心率60到100次/分钟时单次心跳周期约0.6到1秒所以峰值最小间隔至少要设为0.3秒否则会把一个心音的双峰误判成两个独立峰值。二是MinPeakHeight这个需要根据包络幅值自适应设置我通常取包络最大值的0.2到0.3倍作为阈值。[pks, locs] findpeaks(seSmooth, MinPeakDistance, round(0.3*targetFs), ... MinPeakHeight, 0.2 * max(seSmooth), MinPeakProminence, 0.05 * max(seSmooth));这里有个实操技巧MinPeakProminence这个参数也很有用它能过滤掉那些虽然高度高但周围也是高值的平缓峰防止把高幅值包络上的抖动误检成峰值。3.2 S1与S2的区分策略找到所有峰值位置后下一步是把这些峰值区分为S1或S2。S1是二尖瓣和三尖瓣关闭产生的心音出现在心室收缩期开始频率相对较低持续时间稍长S2是主动脉瓣和肺动脉瓣关闭产生的心音出现在心室舒张期开始频率相对较高持续时间较短。单靠单个峰值本身的特征很难稳定区分更可靠的是利用心动周期的时间结构。在有同步心电图的情况下R波前的峰是S1T波后的峰是S2。但这个课题通常只给心音信号所以要用纯心音信号进行区分。常见思路是利用收缩期比舒张期短的特点S1到下一个S2的间隔小于S2到下一个S1的间隔。如果你的应用场景允许这已经足够了。实现时可以用一个简单的聚类思路。把所有相邻峰值间隔分成两类短间隔归为“收缩期间隔”长间隔归为“舒张期间隔”。有了这些间隔规律就能逐步给每个峰打上S1或S2的标签。还有一个备选方法S1的能量通常比S2略集中低频分量占比更高可以对比S1和S2的频谱重心做辅助判断。但这个方法在病理心音上会失效只能作为补充。3.3 特征提取的三层体系分割好了真正的主角才上场特征提取。特征至少要覆盖三个层面我用一个表格来总结特征类别特征名称说明时域特征心音峰值幅度、S1-S2间隔、心率、包络统计量反映心音强度和节律特点计算简单频域特征功率谱重心、频带能量比、谱熵、Mel频率倒谱系数MFCC反映心音频谱构成区分病理杂音时频域特征小波包分解各子带能量比、样本熵同时反映时间与频率变化对非平稳信号更有效时域特征很好理解直接对分割后的S1段和S2段做统计计算。比如S1段幅值均值、标准差、上升时间、持续时间、心率变异性等。这些特征虽然简单但在区分某些瓣膜疾病时依然很有价值。频域特征中最实用的是功率谱重心和频带能量比。功率谱重心就是频谱的质心反映信号能量集中在高频还是低频。频带能量比可以把20到400Hz分成几个子带比如20到80Hz、80到150Hz、150到400Hz分别计算各子带能量占全带能量的比例。这个特征对收缩期杂音的检测特别有效因为正常心音的能量主要落在低频区域高频区域的能量占比很低。MFCC特征是从语音识别领域借用过来的计算步骤涉及预加重、分帧、加窗、FFT、Mel滤波器组、DCT最后得到若干维系数。心音的频率范围虽然比语音低但MFCC的倒谱思想依然适用它能提取出信号频谱的包络形状信息。时频域特征我推荐用小波包分解它比短时傅里叶变换能更好地兼顾时间分辨率和频率分辨率。以db3小波做三层小波包分解为例第3层有8个子频带可以计算每个子频带的能量占整个信号能量的比例作为能量分布特征。这个特征对识别连续性杂音和暂态性杂音有很好的区分度。wname db3; level 3; wpTree wpdec(sigNorm, level, wname); % 提取第3层8个子频带的能量比 energyVector []; for i 0:7 nodeCoeffs wpcoef(wpTree, [level i]); energyVector(i1) sum(nodeCoeffs.^2); end energyRatio energyVector / sum(energyVector);这套特征体系下来单条心音记录通常能提取出50到100维特征。特征多了不一定好我后面专门讲特征筛选。4. 特征选择与降维别让冗余特征拖垮分类器4.1 为什么要做特征选择很多初学者有个误区特征提得越多越好分类器就能学得越好。实际上高维小样本数据集上特征越多越容易过拟合。你只有一两百条样本却提取了上百维特征分类器会把训练集里的噪声都记住测试集上一塌糊涂。这就是所谓的“维度灾难”。此外很多特征之间高度相关。比如时域特征里的峰值幅度和包络统计量可能高度相关频域里的各子带能量比加起来恒等于1本身就携带冗余信息。把这些冗余特征直接送进分类器不仅增加计算量还会干扰分类边界。所以特征选择不是可选项而是必选项。我通常的做法分两步先用过滤式方法粗筛再做PCA降维或进一步精筛。4.2 mRMR筛选兼顾相关性和冗余性特征选择里最有名的方法之一是mRMR也就是最大相关最小冗余。它的核心思想是选择的特征要和类别标签高度相关但彼此之间的冗余性尽量小。这个思路比单纯按方差或者卡方排序更靠谱因为它专门考虑了特征间的相互干扰。MATLAB里没有内置的mRMR函数但可以自己写也可以用File Exchange上的mRMRe实现。核心代码其实不长% 假设 featMat 是样本×特征矩阵label 是样本×1标签 % 计算每个特征与标签的相关性用F统计量或互信息 % 迭代选择特征每次选择与标签相关最大、与已选特征冗余最小的特征我没法在文章里贴一整套mRMR代码因为实现版本很多。但你可以从最简单的思路入手第一轮选F统计量最大的特征之后每一轮对每个候选特征算一个打分打分等于它与标签的互信息减去它与已选特征的平均互信息然后选打分最高的加入。这个思路足以支撑论文里的算法描述。实操中我会先用mRMR筛出20个左右的候选特征再做一次卡方检验或方差分析剔除掉那些p值过高、区分度不明显的特征。这样最终保留的特征组合通常能稳定地达到不错的分类效果。4.3 PCA降维的实际操作与累计贡献率筛选完特征后如果维度还在20左右而样本只有几十条我还会再做一次PCA降维。PCA通过线性变换把高维特征映射到彼此正交的少数几个主成分上本质上是把原来的相关特征提取成新的综合特征。% featSelected 是筛选后的特征矩阵 [coeff, score, latent, tsquared, explained] pca(featSelected); % 查看累计贡献率 cumsum(explained) % 选择累计贡献率达到95%的主成分 cumContribution cumsum(explained); nComps find(cumContribution 95, 1); featPCA score(:, 1:nComps);这里有一个重要的提醒PCA会改变特征的可解释性降维后的主成分不再是原始物理特征而是一种数学组合。如果你的论文重点是心音病理机理分析那就应该多依赖mRMR筛选出的物理特征而不是直接用PCA结果。如果重点是构建一个分类系统PCA完全可以用。还记得最重要的一条PCA必须在训练集上做然后把同样的变换矩阵应用在测试集上。千万不要在全部数据上一起做PCA这样会把测试集信息泄露给训练过程得到的评估指标会虚高答辩时评委稍微一追问就露馅了。正确做法是用cvpartition划分训练测试后在训练集上计算PCA系数再通过center和score变换测试集。5. 分类识别与系统评估5.1 分类器选型对比到这一步你已经得到了一个干净的特征矩阵接下来就是正常和异常心音的分类。心音分类常用的分类器有这么几个支持向量机SVM、k最近邻kNN、随机森林、朴素贝叶斯等。我的建议是优先试SVM和随机森林。SVM对中小样本分类效果很好尤其配合RBF核函数能拟合非线性边界。随机森林的好处是它能输出特征重要性你可以借它反向验证前面特征选择的合理性。kNN虽然简单但在特征维数不高、样本量不大时效果也还行缺点是分类边界不够平滑对噪声敏感。MATLAB里用fitcecoc做多分类SVM但心音分类通常是二分类直接用fitcsvm就够了。% featPCA 是降维后的特征label 是二分类标签 SVMModel fitcsvm(featPCA, label, ... KernelFunction, rbf, ... BoxConstraint, 1, ... Standardize, true); % 十折交叉验证评估 cvModel crossval(SVMModel, KFold, 10); predLabel kfoldPredict(cvModel);这里的BoxConstraint是SVM的惩罚参数默认是1。你可以用bayesopt或fitcsvm配合OptimizeHyperparameters做自动调参但要注意样本量太少时超参数寻优很容易过拟合所以不要一味追求交叉验证精度。随机森林在MATLAB里对应的是TreeBagger也可以叫fitcensemble选Bag方法。训练完以后可以直接查看特征重要性rfModel TreeBagger(100, featSelected, label, Method, classification, ... OOBPrediction, on); importance rfModel.OOBPermutedPredictorDeltaError;5.2 评估指标的解读分类模型做好以后不能只看一个准确率。心音分类有个特点正常样本通常多于异常样本如果异常样本只占20%那就算你全部预测成正常准确率也有80%。但这种模型没有任何临床价值。所以必须同时关注敏感度、特异度、阳性预测值和F1分数。我建议在测试阶段至少输出以下指标指标公式意义准确率(TPTN)/(TPTNFPFN)整体判对的比例敏感度TP/(TPFN)异常样本有多少被正确检出特异度TN/(TNFP)正常样本有多少被正确识别F1分数2PrecisionRecall/(PrecisionRecall)精确率和召回率的平衡这些指标可以通过confusionmat先求出混淆矩阵再手动计算也可以用perfcurve画ROC曲线并计算AUC。ROC曲线和AUC几乎是答辩时的必展示图跟频谱图和包络图一样重要。6. 常见问题与排查技巧实录6.1 常见问题速查表我把这几年做这个课题以及帮别人翻车后总结的常见问题整理成一张速查表建议收藏之后慢慢对照。问题表现可能原因解决方案滤波后心音完全没波形带通范围太窄或截止频率设错画频谱图确认心音频段调整带通范围峰值检测到一大堆冗余峰值MinPeakDistance设置太小按心率下限估算最小峰间距至少0.3秒S1和S2总区分不开舒张期和收缩期间隔差异不明显引入频谱特征辅助判断同一段信号每次运行结果不同随机划分或随机种子没有固定用rng固定随机种子设置cvpartition分类准确率很高但曲线图上表现为过拟合数据泄漏例如全数据PCA后再交叉验证确保所有降维、归一化都在训练集上完成特征矩阵里出现NaN或Inf信号截取长度不一致或除零检查每段特征提取的窗口长度统一样本长度6.2 我在实际项目中反复踩过的坑分享两个最有代表性的坑。第一个是包络双峰问题。正常心音的S1和S2本身由多个瓣膜分量叠加在包络上可能出现双峰甚至三峰。如果你用简单的峰值检测同一个心音会被识别成两三个相邻峰值后续间隔特征全乱。我一开始用单纯的高度阈值怎么调都调不干净。后来发现两个参数组合能解决大部分问题MinPeakDistance设成0.3倍采样率加上MinPeakProminence设为最大包络高度的0.05倍。前者保证两个峰之间必须有足够时间间隔后者确保即使峰高于阈值如果它周围都是高包络区域也会被过滤掉。这个细节是跑几十条样本之后才总结出来的直接决定检测结果稳定性。第二个是样本长度不一致的问题。公开数据集里有些心音记录只有10秒有些有30秒乃至1分钟但你提取特征时如果按整段信号做统计那特征根本没有可比性。最开始我就踩了这个坑把长短不一的信号丢进特征提取函数结果特征矩阵里混合了不同数量心跳周期的统计量。后来我改为统一提取每个心动周期内的特征再取平均或者只截取固定长度比如10秒的稳定段做分析。这个方法让特征更规范也大幅度提升了分类指标。6.3 后续可扩展的思路如果你完成基础版本后还有余力可以考虑往两个方向扩展。一个是把SVM换成深度学习模型比如用一维卷积神经网络直接对心音波形做端到端分类不需要手工提取特征。但深度学习需要大量数据公开的心音数据集样本量有限直接用CNN很容易过拟合迁移学习的空间也不大所以我不建议作为毕设首选。另一个方向是做成实时系统用音频采集设备连续采集心音信号在滑动窗口上实时提取特征并判断是否有异常这个方向适合用MATLAB的App Designer做个GUI界面展示实时波形和分类结果答辩效果会很好。最后再给一句个人建议做的心音课题一定要以能完整复现为底线。把分割的结果画出来把特征提取的代码封装成函数把分类器的参数记录好。这样就算答辩时被老师问到了自己不熟的地方也能靠扎实的调试过程和可视化作支撑。毕竟心音信号处理这个方向水很深但骨架就摆在这里把它老老实实做扎实就已经超出大多数人一截了。本文还有配套的精品资源点击获取
返回列表