心音信号处理与分割:Python实现包络提取到S1/S2识别
2026/9/14 14:43:38 网站建设 项目流程

简介:心音信号分析在心脏病早期诊断中具有重要价值,尤其适用于瓣膜功能评估与心律失常识别。面向心音信号处理及生物医学工程研究者,这套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(parents=True, exist_ok=True) with zipfile.ZipFile(source) as zf: # testzip() 会逐条读取并校验 CRC,返回第一个损坏文件名 corrupt = zf.testzip() if corrupt is not None: raise RuntimeError(f"zip 内部文件校验失败: {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_sizecompress_sizeCRCis_dir()等元数据,便于后续做大小统计。取Path(member.filename).name是一种防御式写法:zip 内部文件名可能是../../evil.wav,直接拼接会让文件写到解压目录之外,这一点在公开数据集里不常见,但在内部流传的数据包里偶尔会出现。

ZipFile默认不处理加密成员打开,遇到加密 zip 需要先zf.setpassword(b"password"),所有成员读取时都会用这个密码。需要注意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(parents=True, exist_ok=True) records = [] files = sorted(src_dir.glob("*.wav")) for idx, f in enumerate(files, start=1): 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"], keep=False) manifest.to_csv(dst_dir / "manifest.csv", index=False, encoding="utf-8-sig") return manifest

这里用soundfile读取音频信息,比wave模块省事得多,而且本身能处理samplerateframeschannels等字段。采样率不同会直接影响后续滤波参数,所以放在文件名里。同时生成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)], btype="bandpass") 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, mode="same") return env

代码里第一个参数是lowcuthighcut,分别对应带通滤波的上下截止频率。心音 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低频噪声多时提高 lowcut
8000 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,后续固定高度阈值(例如height=0.3 * np.max(env))对每个文件都适用。需要记住的是归一化只影响包络的幅度分布,不影响时间位置,所以分割点可以同步映射回原始波形。

还要记录包络的采样率。后面做峰值分割时,所有距离、宽度参数都需要按采样率换算成样本数。很多新人直接把以秒为单位的数值传给find_peaksdistance,得到的结果在采样率不同的文件间完全不可用。我的习惯是保存归一化包络时同时保存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, mode="full") 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, height=height, distance=min_distance, prominence=0.15 * np.max(env_norm), width=(1, int(0.2 * fs)), ) # 返回峰位置和属性,属性里包含 peaks 处的幅值、半宽等 return peaks, props

height我设为第 80 百分位,这比固定 0.5 更贴近数据分布。distance必须设为至少 0.4 倍心动周期,防止把同一心搏的 S1 或 S2 里的双峰错误识别成两个独立峰。prominence是突出度,衡量峰相对于两侧局部最小值的突出程度,对于幅度很小的 S2 但局部仍突出的情况,突出度比高度阈值更敏感。width参数限制峰的宽度在 1~0.2 秒之间,避免把宽大的呼吸音包络当成心音峰。

参数推荐值作用调参方向
height第 80 百分位过滤低幅度噪声峰噪声多时提高到 90
distance0.4 * est_lag防止双峰被分成两个峰心率快时减小到 0.3
prominence0.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列标注okuncertainbad,训练时只取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(parents=True, exist_ok=True) 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", index=False, encoding="utf-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_peaksdistance是否与心率估计值用的单位一致。最后一个我常用的技巧:把包络向下平移 0.1 后和原始波形叠加,再画峰位置,这样能同时看到分割点是否切在原始波形的冲击起始位置,而不是包络上的滞后位置。如果发现滞后的样本数大约等于移动平均窗口的一半,那就是预期的滑动平均相位延迟,把峰位置减去这个偏移就能恢复到原始波形上的准确时刻。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询