ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

心音信号处理与分割:Python实现包络提取到S1/S2识别

心音信号处理与分割:Python实现包络提取到S1/S2识别 简介心音信号分析在心脏病早期诊断中具有重要价值尤其适用于瓣膜功能评估与心律失常识别。面向心音信号处理及生物医学工程研究者这套MATLAB代码包聚焦信号分割与包络提取两项核心任务为解决心音时相定位、包络特征分析等提供可直接运行的轻量级实现。第一、第二心音分别对应收缩期与舒张期准确分割与包络分析是识别杂音、评估瓣膜状态的重要前提。包内共4个文件由两个核心脚本与两个自动保存备份组成分别实现心音包络提取与降噪、信号分割整体体积仅2KB结构清晰适合二次开发。已有181人学习/下载适合信号处理方向的学生在课程设计或入门实验中参考。借助这些脚本使用者可快速了解从原始心音数据到包络曲线与分段结果的处理流程为心率计算、瓣膜功能评估等后续分析奠定基础也可作为教学演示中的实用示例。1. 心音信号文件到手之后先别急着解压跑模型拿到一个名为“新建文件夹.zip”的压缩包里面是几十个心音信号文件文件名是一串乱序时间戳或纯数字这就是心音信号分析任务的常规开局。心音信号PCG记录的是心脏瓣膜与血流振动产生的微弱声学信号核心任务通常是把连续录制切成单心搏再进一步区分 S1、S2 对应的收缩期与舒张期。这个“信号分割 包络提取”的前处理链质量直接决定后续分类模型的准确率也让压缩包里的零散文件变成结构化的训练集。这篇内容用 Python 生态把解压、校验、包络提取到分割落盘完整走一遍参数和坑都标注清楚适合生物信号处理工程师和数据工程师照着改。很多人一上来就用wavfile.read读完直接画包络结果发现分割边界飘得离谱。问题不在算法而在上游压缩包里的文件可能是多通道录音某个通道对应心音另一个通道是对照用的心电ECG采样率可能是 2000 Hz 或 8000 Hz滤波参数不能照搬论文。下面这套流程是我处理多个 PCG 数据集时的基线版本每一步都能独立验证也可以按自己的数据调整。2. 心音信号压缩包的解压与文件校验先把文件系统理干净2.1 用 Python zipfile 解压心音 zip先做 CRC 校验再落盘拿到 zip 后的第一步不是直接提取而是先检查这个包是不是完整、内部有没有恶意路径。心音信号文件通常来自医院或采集设备压缩包可能由其他人用 Windows 的“发送到压缩文件夹”生成里面经常混有__MACOSX目录、隐藏的.DS_Store文件甚至文件路径带盘符。如果直接把member.filename拼到解压目录就可能出现路径穿越把文件写到磁盘任意位置。import zipfile from pathlib import Path def extract_pcg_zip(zip_path: str, out_dir: str) - Path: source Path(zip_path) target Path(out_dir) target.mkdir(parentsTrue, exist_okTrue) with zipfile.ZipFile(source) as zf: # testzip() 会逐条读取并校验 CRC返回第一个损坏文件名 corrupt zf.testzip() if corrupt is not None: raise RuntimeError(fzip 内部文件校验失败: {corrupt}) for member in zf.infolist(): if member.is_dir(): continue # 只取文件名最后一段避免 ../../ 等路径穿越 clean_name Path(member.filename).name if not clean_name: continue dest target / clean_name with zf.open(member, r) as src, open(dest, wb) as dst: dst.write(src.read()) return target这个函数里最关键的是testzip()它会把 zip 内的每个压缩成员完整读出来并和中央目录里记录的 CRC32 对比。如果网络下载时丢字节testzip()会返回坏条目这时候整个解压没有继续的意义提示对方重新压缩或换个渠道获取。之后用infolist()而不是namelist()因为infolist()返回ZipInfo对象包含file_size、compress_size、CRC、is_dir()等元数据便于后续做大小统计。取Path(member.filename).name是一种防御式写法zip 内部文件名可能是../../evil.wav直接拼接会让文件写到解压目录之外这一点在公开数据集里不常见但在内部流传的数据包里偶尔会出现。ZipFile默认不处理加密成员打开遇到加密 zip 需要先zf.setpassword(bpassword)所有成员读取时都会用这个密码。需要注意zipfile模块不支持 AES 加密格式只支持传统的 ZipCrypto。如果打开时抛出RuntimeError: File ... is encrypted说明这个包是 AES 加密的常见做法是换成 7-Zip 的命令行工具7z x -p密码 archive.zip来解压或者直接联系数据提供方要一份无密码副本。网上那些“zip 压缩包密码破解工具”大多针对弱密码的传统 ZipCrypto跑去暴力破解心音数据属于浪费时间不如先把密码问清楚。2.2 解压后的文件校验和批量重命名把心音信号文件整理成可追踪清单解压后你会发现文件名基本不能用比如REC001.WAV、新建文件夹 (2).wav或者干脆是20240613_093512.wav这种时间戳。如果不重命名后面的分割结果无法和原始记录关联。我一般的做法是扫描所有 wav 文件读采样率和时长按“研究编号_设备编号_序号”重命名并生成一个 CSV 清单。import soundfile as sf import pandas as pd from pathlib import Path def inspect_and_rename(src_dir: Path, dst_dir: Path, prefix: str PCG): dst_dir.mkdir(parentsTrue, exist_okTrue) records [] files sorted(src_dir.glob(*.wav)) for idx, f in enumerate(files, start1): info sf.info(str(f)) dest_name f{prefix}_{idx:04d}_{info.samplerate:.0f}Hz.wav dest_path dst_dir / dest_name dest_path.write_bytes(f.read_bytes()) # 简化版不保留元数据 records.append({ src_name: f.name, new_name: dest_name, samplerate: info.samplerate, frames: info.frames, duration_s: round(info.frames / info.samplerate, 3), }) manifest pd.DataFrame(records) # 标记重复的采样率时长组合不直接删除 manifest[dup_flag] manifest.duplicated( subset[samplerate, duration_s], keepFalse) manifest.to_csv(dst_dir / manifest.csv, indexFalse, encodingutf-8-sig) return manifest这里用soundfile读取音频信息比wave模块省事得多而且本身能处理samplerate、frames、channels等字段。采样率不同会直接影响后续滤波参数所以放在文件名里。同时生成manifest.csv记录原始文件名和新文件名。检查列表时有一个很容易忽略的点用dup_flag只是提醒你可能有重复录音真正去重需要比较波形相似度而不是只看时长和采样率。时长相同但内容不同的文件很常见尤其是静音段较多的心音记录。文件格式读取库常见采样率典型来源.wavsoundfile / wave8000 / 16000 / 22050 / 44100电子听诊器或麦克风采集.npynumpy.load无固定中间处理结果保存为数组.matscipy.io.loadmat取决于写入时科研设备导出的结构体数据扫描过程里还能顺手发现坏文件。读取sf.info()时如果报Error opening file通常是 wav 文件头损坏原始记录未正常写入应该从 zip 原始包中找对应字节做 CRC 对比而不是直接丢弃。如果 zip 解压时就报了error read zip archive大概率是压缩包没有下载完整重下源文件通常比修复更省时间。3. 心音信号的包络提取从原始波形到可分割的包络线3.1 为什么心音包络要先滤波再提取心音的频率范围在 20~400 Hz 左右但录音时呼吸音、肌肉音、环境噪声会叠加进来。直接用原始幅度做峰值搜索时一次咳嗽就能产生大量伪峰包络也不收敛。所以包络提取必须先做带通滤波限制频段。滤波器的阶数和截止频率要按采样率换算而不是把论文里的数字硬套。包络的常见定义是解析信号的幅度。对一个实信号做希尔伯特变换得到解析信号 (x(t) jH(x(t)))取模就是瞬时包络。它能反映信号的能量走势比原始波形更平滑尤其适合心音这种由 S1、S2 两个主要冲击组成的声音。另一种常见做法是移动平均绝对值先把信号取绝对值再用一个长度为 10~30 ms 的窗口做卷积算法更简单但时间分辨率比希尔伯特变换要差一点。我在处理采样率不齐的多个心音信号文件时默认用希尔伯特变换原因只是参数少截止频率和窗长定了结果基本稳定。3.2 希尔伯特包络的 Python 实现与参数表下面是提取单个心音信号文件包络的完整函数输入是原始波形x输出是滤波平滑后的包络长度与输入一致。import numpy as np from scipy import signal def extract_pcg_envelope(x: np.ndarray, fs: int, lowcut: float 25, highcut: float 400, order: int 4) - np.ndarray: # 1. 带通滤波滤除 DC 漂移和高频噪声 b, a signal.butter(order, [lowcut / (fs / 2), highcut / (fs / 2)], btypebandpass) y signal.filtfilt(b, a, x) # 2. 希尔伯特变换取包络 analytic signal.hilbert(y) env np.abs(analytic) # 3. 再做一次 20 ms 移动平均平滑消除包络上的毛刺 win_len max(1, int(0.02 * fs)) kernel np.ones(win_len) / win_len env np.convolve(env, kernel, modesame) return env代码里第一个参数是lowcut和highcut分别对应带通滤波的上下截止频率。心音 S1 和 S2 的频谱分布在 20~200 Hz某些高频成分能到 400 Hz所以常规设置是 25 Hz 和 400 Hz。如果录音设备是电子听诊器低频噪声较多可以把lowcut提到 50 Hz如果是传感器紧贴皮肤呼吸音明显可以考虑把highcut降到 300 Hz。但注意不能低于 200 Hz否则 S1 和 S2 的瞬态成分都会被滤掉包络变胖分割边界反而更模糊。order是巴特沃斯滤波器阶数四阶是一个折中。阶数越高衰减越陡但filtfilt的相位畸变越小计算开销也增加。在 8000 Hz 采样率下四阶完全够用。scipy.signal.hilbert默认对最后一个轴做全长度变换返回复解析信号。取绝对值得到包络后包络上还会有高频纹波比如 S1 内部两瓣造成的双峰。20 ms 移动平均窗口能把这些毛刺压平但窗口太大会让包络峰变宽导致后面峰值定位不准。经验值采样率 2000 Hz 对应窗长 40 个点采样率 8000 Hz 对应窗长 160 个点这个参数在大多数 PCG 数据上都能稳定工作。采样率 fs20ms 窗长25~400Hz 的归一化频率调参建议2000 Hz40 点0.025 ~ 0.4低频噪声多时提高 lowcut8000 Hz160 点0.00625 ~ 0.1常见电子听诊器16000 Hz320 点0.0031 ~ 0.05先降采样到 8000 Hz 再处理3.3 包络归一化处理来自不同设备的增益差异从不同型号的电子听诊器拿到的心音信号文件幅度增益可能差十倍。如果不做归一化就到分割阶段峰值阈值就得为每个文件单独调没法批量跑。我一般用峰值归一化把包络的最大值映射到 1.0再乘一个比例系数 0.8 留出冗余env_norm env / np.max(env) * 0.8更稳健的做法是用百分位数而不是最大值来归一化比如用np.percentile(env, 99)作为分母这样能避免个别异常大脉冲把整体压得过低。归一化之后包络范围就是 0~0.8后续固定高度阈值例如height0.3 * np.max(env)对每个文件都适用。需要记住的是归一化只影响包络的幅度分布不影响时间位置所以分割点可以同步映射回原始波形。还要记录包络的采样率。后面做峰值分割时所有距离、宽度参数都需要按采样率换算成样本数。很多新人直接把以秒为单位的数值传给find_peaks的distance得到的结果在采样率不同的文件间完全不可用。我的习惯是保存归一化包络时同时保存fs到一个 JSON 元数据文件这样重跑实验时不用回头看原始 wav 头。4. 心音信号分割基于包络的心脏周期切分与校正4.1 从包络峰值中区分 S1 与 S2先估心动周期分割的目标是找到每个心搏的 S1 和 S2。一个心动周期内S1 之后是收缩期一般 0.2~0.4 秒接着是 S2然后舒张期一般 0.4~0.6 秒再进入下一个周期。也就是说S1 与 S2 的间隔明显小于连续两个 S1 的间隔。因此不能只找包络峰值还要识别峰的属性。最先要估计的是心率也就是平均心动周期。对包络做自相关自相关函数的第一个显著峰值对应基频其倒数是平均周期。代码实现def estimate_heart_cycle(env: np.ndarray, fs: int, min_bpm: float 40, max_bpm: float 200) - float: env env - np.mean(env) # 去直流 corr np.correlate(env, env, modefull) corr corr[len(corr)//2:] # 只取正延迟部分 min_lag int(60 * fs / max_bpm) max_lag int(60 * fs / min_bpm) segment corr[min_lag:max_lag] peak_pos np.argmax(segment) lag min_lag peak_pos print(f估计心动周期: {lag / fs:.2f} s, 心率: {60.0 * fs / lag:.1f} bpm) return float(lag)自相关函数的峰值位置表示信号和自身平移后最接近的距离这个距离就是心动周期。用限制心率范围的方式减小出错概率。计算得到 lag 之后接下来用scipy.signal.find_peaks找包络峰值并确保相邻峰距离不小于 0.4 倍 lag这样能滤掉与心动周期无关的伪峰。4.2 用 scipy.signal.find_peaks 切分心搏height、distance、prominence 的联动下面这段代码是核心分割流程。输入是归一化包络env_norm和采样率fs输出是 S1/S2 候选峰的位置和属性。from scipy.signal import find_peaks def detect_s1_s2_peaks(env_norm: np.ndarray, fs: int, est_lag: float): min_distance int(0.4 * est_lag) height float(np.percentile(env_norm, 80)) peaks, props find_peaks( env_norm, heightheight, distancemin_distance, prominence0.15 * np.max(env_norm), width(1, int(0.2 * fs)), ) # 返回峰位置和属性属性里包含 peaks 处的幅值、半宽等 return peaks, propsheight我设为第 80 百分位这比固定 0.5 更贴近数据分布。distance必须设为至少 0.4 倍心动周期防止把同一心搏的 S1 或 S2 里的双峰错误识别成两个独立峰。prominence是突出度衡量峰相对于两侧局部最小值的突出程度对于幅度很小的 S2 但局部仍突出的情况突出度比高度阈值更敏感。width参数限制峰的宽度在 1~0.2 秒之间避免把宽大的呼吸音包络当成心音峰。参数推荐值作用调参方向height第 80 百分位过滤低幅度噪声峰噪声多时提高到 90distance0.4 * est_lag防止双峰被分成两个峰心率快时减小到 0.3prominence0.15 * max(env)识别突出的小峰S2 幅度低时降低width1 ~ 0.2 秒排除宽大呼吸峰呼吸干扰大时收窄得到峰位置后还需要给每个峰标记是 S1 还是 S2。常见做法是以相邻两个 S1 之间包含恰好一个 S2 为约束。从峰序列中每两个连续峰计算间隔如果间隔接近心动周期的一半左右那么前一个峰可能是 S1后一个可能是 S2如果连续三个峰间隔大致相等说明中间被噪声污染或出现早搏此时采用回溯校正先假设第一个峰是 S1然后每隔一个周期找另一个峰作为下一个 S1在两个 S1 之间取包络最大值作为 S2。这段逻辑可以写成一个小函数def assign_s1_s2(peaks: np.ndarray, env_norm: np.ndarray, est_lag: float, fs: int): s1_list, s2_list [], [] prev_peak None for p in peaks: if prev_peak is None: # 第一个峰暂时标记为 S1 s1_list.append(p) else: # 如果当前峰到前一个峰的距离远大于半周期则视为下一 S1 gap (p - prev_peak) / fs if gap 0.75 * (est_lag / fs): s1_list.append(p) else: s2_list.append(p) prev_peak p return np.array(s1_list), np.array(s2_list)这个启发式方法不完美但对平静呼吸、心律齐的常规心音信号文件已经足够。心衰患者或严重心律失常时心搏间期不等S1/S2 的间隔也会变化直接按固定比例划分会出错这时需要做更精细的聚类或引入隐马尔可夫模型。这里不展开因为大多数科研序列从正常数据起步。4.3 分割失败时的三个常见坑双峰、呼吸滑音和缺失心搏第一个坑S1 在包络上出现双峰。心音信号中 S1 本身由二尖瓣和三尖瓣关闭产生两个成分间隔只有 20~30 ms在平滑不充分的包络上会形成两个小峰。我们的distance参数按 0.4 倍心动周期设置后双峰会被强制合并成一个峰因为第二个峰距离太近被忽略。第二个坑呼吸音导致的低频包络隆起。吸气时包络整体抬高但没有明显峰形height阈值可能把它变成一个宽峰。我们依靠width限制峰宽同时额外检查峰两侧的对称性心音峰通常是陡峭上升再快速下降呼吸峰则坡度平缓。可以在峰值位置前后取几毫秒的斜率做判断。第三个坑早搏或漏搏导致的缺失心搏。连续几个心动周期规律后突然出现一个长间歇包络上没有对应峰值这时按照心跳周期外推一个虚拟分割点但这部分数据最好标记为“不确定”不要强行放进训练集。我的处理是所有分割区间写入表格并在quality列标注ok、uncertain和bad训练时只取ok样本。5. 把分割结果批量持久化从心音信号文件到结构化数据集5.1 将心搏片段保存为单独 wav 并生成标签 CSV分割完成后通常要把每个心搏的 S1、S2 片段切出来保存成独立的 wav 文件方便后续做分类或特征工程。下面是我常用的落盘函数import soundfile as sf from pathlib import Path import pandas as pd def export_segments(wave: np.ndarray, fs: int, s1_list: np.ndarray, s2_list: np.ndarray, out_dir: Path, patient_id: str): out_dir.mkdir(parentsTrue, exist_okTrue) rows [] for idx, (s1, s2) in enumerate(zip(s1_list, s2_list)): # 以 S1 起点为中心向两侧扩展 0.2 s start max(0, int(s1 - 0.2 * fs)) # 取到下一 S2 的起点 end min(len(wave), int(s2 0.2 * fs)) seg wave[start:end] name f{patient_id}_{idx:04d}_S1S2.wav sf.write(out_dir / name, seg, fs) rows.append({file: name, start_s: start / fs, end_s: end / fs, duration_s: (end - start) / fs}) pd.DataFrame(rows).to_csv(out_dir / segments.csv, indexFalse, encodingutf-8-sig)时间窗的选取有没有统一标准S1 起点向前 0.2 秒是为了捕获 S1 的起始瞬态向后到 S2 起点之后 0.2 秒是为了完整包含整个收缩期和舒张期。如果是为了做心音分类实际只需要 S1 和 S2 各自的 0.1 秒片段如果是为了心率变异性分析则需要保留完整周期。文件命名里包含patient_id和索引这样即使片段文件很多也能从文件名直接看出属于哪个病例、在哪个位置。CSV 里记录了从原始波形切取的起止时间方便回溯。这样一批心音信号文件就能变成标准的监督学习数据集。5.2 快速验证分割质量的三个指标落盘后别急着删原始文件先用三个指标做一次 sanity check。第一分割出的心搏数量是否在预期范围心搏数 录音时长(秒) / 平均心动周期(秒)误差超过 15% 说明有漏检或过检。第二S1 与 S2 的时间间隔是否符合生理范围收缩期一般在 0.2~0.4 秒若某个片段的 S1-S2 间隔超过 0.5 秒且反复出现大概率是分割错误。第三相邻 S1 间隔的变异系数CV正常呼吸影响下变异系数在 5%~15% 之间如果接近 30%说明文件可能包含大量运动伪影需要重新滤波。验证时可以直接用librosa.play播放随机抽取的 10 个片段听感比任何指标都直观。把env_norm和峰位置画在一张图上缩放显示前 5 秒一旦发现峰位置落在包络谷底检查find_peaks的distance是否与心率估计值用的单位一致。最后一个我常用的技巧把包络向下平移 0.1 后和原始波形叠加再画峰位置这样能同时看到分割点是否切在原始波形的冲击起始位置而不是包络上的滞后位置。如果发现滞后的样本数大约等于移动平均窗口的一半那就是预期的滑动平均相位延迟把峰位置减去这个偏移就能恢复到原始波形上的准确时刻。本文还有配套的精品资源点击获取
返回列表