ARTICLE DETAIL

建站实战干货

来自一线的建站与推广经验沉淀,每一条都经过真实交付验证。

PPG/ECG信号预处理与特征提取:NeuroKit2完整工程实践

2026/9/28 7:04:19 拓冰建站 浏览量
PPG/ECG信号预处理与特征提取:NeuroKit2完整工程实践 简介压缩包面向生物医学工程、信号处理方向的开发者与学生提供一套可直接运行的PPG与ECG同步数据预处理及特征提取代码。方案基于nk库对原始信号去噪进而计算潮波幅值比h2/h1、重搏波幅值比h4/h1、收缩面积比S1/S、舒张面积比S2/S、主波高度h1、波形周期面积S、收缩期面积S1、舒张期面积S2等脉搏波形态学参数并利用ECG与PPG之间的时间差得到脉搏传导时间PWTT可用于无创评估动脉弹性与血管阻力相关研究。资源共43个文件主要包含2个Python脚本、29个txt原始信号数据、8个ini配置及4个bak备份文件压缩包整体大小仅4.31MB轻量且便于迁移。目前已有744人浏览学习代码注释清晰、结构完整可直接替换数据文件跑通全流程适合作为信号预处理与特征提取的参考实现也能够在此基础上扩展其他心血管生理指标计算。1. PPG 与 ECG 预处理及特征提取一份可以直接跑的完整工程做脉搏波和心电信号分析的人大概率都经历过这种阶段数据从采集设备导出来是一堆 txt波形里混着基线漂移和肌电噪声特征值算出来自己都不敢信。这份ppg_ecg.rar就是冲着这个痛点来的——里面既有 20 多组真实的 PPG/ECG 同步数据又有两段可以直接运行的 Python 代码一段用 NeuroKit2 库做去噪预处理一段做潮波幅值比、重搏波幅值比、收缩舒张面积比、脉搏传导时间等特征提取。适合正在做心血管功能评估、动脉硬化分析、心理生理信号处理相关课题的人拿到手改改数据路径就能跑通自己的数据不用从零搭信号处理流水线。我拿到后先整体拆了一遍发现它对新手友好但对参数不敏感的人容易踩采样率和峰检测的坑这篇文章就按「工程结构 → 预处理 → 特征提取 → 避坑 → PWTT 计算细节」的顺序把它彻底讲透。2. 工程结构与数据格式先看清 txt 里存的是什么2.1 文件清单与职责划分解压ppg_ecg.rar后核心文件就两类Python 脚本和数据文件。我先整理一个文件清单让你对工程规模有个底。文件/目录类型作用feature.pyPython 脚本特征提取主程序计算波形形态特征与面积特征nk2.pyPython 脚本基于 NeuroKit2 的信号预处理完成去噪和峰值定位Data_3644.txt等 20 个 txt数据文件每组包含同步的 PPG 和 ECG 信号.spyproject/configSpyder 工程配置不影响核心逻辑可忽略从文件命名看数据是按采集编号组织而非按受试者编号Data_3647 2.txt这个文件名里带空格和数字后缀应该是某次采集的重复保存文件。实际使用时不建议修改原始文件名因为feature.py里可能用文件名作为输出标签随意改名会导致结果对不上号。2.2 数据格式判定分离符与通道顺序打开一个 txt 文件先用文本编辑器看前几行。项目里没写 README所以通道顺序需要自己判断这也是拿到陌生数据的第一件事。head -n 5 Data_3644.txt常见两种格式一种是一行两列tab 或空格分隔左列为 PPG右列为 ECG另一种是一行多列可能包含时间戳。这个项目里的数据是同步采集大概率是两列但分隔符可能是空格也可能是制表符直接用np.loadtxt读取时要注意。import numpy as np data np.loadtxt(Data_3644.txt) print(data.shape) print(data[0:5, :])如果报错九成是分隔符问题改成np.loadtxt(Data_3644.txt, delimiter\t)或delimiter,再试。拿到data后画一段波形看一眼通道顺序别假设第一列一定是 PPG。判断方法很直接ECG 的 QRS 波是尖锐的窄脉冲PPG 的波形是平滑的周期起伏看一眼就分得出来。这一步花不了两分钟但能避免后面所有特征全部算反的灾难。3. nk2.py 预处理用 NeuroKit2 把噪声从信号里抠出去3.1 为什么选 NeuroKit2 而不是手写滤波器nk2.py的核心是nk库也就是 NeuroKit2。它处理 PPG 和 ECG 是高度自动化的你给它原始信号和采样率它自己完成滤波去噪、峰值检测甚至能输出质量评估。自己手写滤波器当然也可以但 PPG 的滤波涉及去除基线漂移低频和肌电噪声高频两个截止频率要反复试ECG 更麻烦要保留 QRS 波的同时滤掉 P 波和 T 波以外的干扰手调参数会调到头秃。NK 库的滤波思路是先做带通滤波再用专门的峰值检测算法去找每个周期。不同信号类型对应不同算法ECG 用的是基于 Pan-Tompkins 的改进算法PPG 用的是基于自适应阈值的算法这些细节都在库内部处理了工程上省很多事。3.2 预处理主流程拆解nk2.py里比较核心的是nk.ppg_process和nk.ecg_process这两个函数它们把滤波和峰值检测打包在一起。完整流程可以拆成三步读取数据 → 分别处理 PPG 和 ECG → 把结果合并对齐。import neurokit2 as nk import numpy as np # 假设 data 已从 txt 读入通道0为PPG通道1为ECG ppg_raw data[:, 0] ecg_raw data[:, 1] sample_rate 1000 # 注意必须和数据采集时的实际采样率一致 # PPG 预处理这一步同时完成去噪和峰值检测 ppg_signals, ppg_info nk.ppg_process(ppg_raw, sampling_ratesample_rate) # ECG 预处理同样集成滤波和 R 峰定位 ecg_signals, ecg_info nk.ecg_process(ecg_raw, sampling_ratesample_rate) # 提取 PPG 波形峰值位置单位样本点索引 ppg_peaks ppg_info[PPG_Peaks] # 提取 ECG 的 R 峰位置 ecg_rpeaks ecg_info[ECG_R_Peaks] # nk.ppg_process 返回的干净信号列名为 PPG_CleanECG 为 ECG_Clean ppg_clean ppg_signals[PPG_Clean].values ecg_clean ecg_signals[ECG_Clean].values逻辑说明nk.ppg_process内部先对原始 PPG 做带通滤波默认 0.5~8Hz 左右适合脉搏波频段再用专有算法定位每个脉搏波的起点和峰值点nk.ecg_process同理但滤波频段更宽0.5~150Hz 左右以保留 QRS 波的陡峭形态。两个函数返回的都是 DataFrame*_Clean列是滤波后的信号*_info字典里是峰值索引数组。参数说明最关键的是sampling_rate一定要填对。如果采集设备是 125Hz 或 500Hz填了 1000Hz那么算出来的所有时间相关特征PWTT 尤其敏感都会偏差巨大。判断采样率的方法有两个看采集设备说明书或者直接数信号长度除以采集时长秒。3.3 采样率不一致时的处理有些数据来源混杂采样率不统一直接全部按同一个采样率处理会导致特征值横向不可比。建议在读取数据后加一段自适应判断逻辑# 估算采样率如果知道总时长 duration_seconds duration_sec 60 # 假设这段数据时长60秒 estimated_fs data.shape[0] / duration_sec print(f估算采样率: {estimated_fs:.1f} Hz) # 如果数据来源采样率不同需要重采样到统一采样率 from scipy import signal as scipy_signal if abs(estimated_fs - 1000) 1: resample_ratio 1000 / estimated_fs target_len int(data.shape[0] * resample_ratio) ppg_resampled scipy_signal.resample(ppg_raw, target_len) ecg_resampled scipy_signal.resample(ecg_raw, target_len)重采样到统一采样率的做法是为了保证后续所有特征在时间维度上一致可比。scipy.signal.resample是傅里叶变换法重采样对信号形状保存较好但要注意它默认输出复数需要取实部。3.4 一个容易翻车的点峰值检测结果不检查预处理跑完不是终点必须检查峰值检测结果。NK 库的自动峰值算法对质量好的信号准确率高但遇到运动伪迹严重的数据段误检和漏检几乎是必然的。我一般会做个可视化快速扫一眼import matplotlib.pyplot as plt # 截取前10秒信号做可视化检查 plot_sec 10 plot_samples plot_sec * sample_rate plt.figure(figsize(16, 8)) plt.subplot(2, 1, 1) plt.plot(ppg_raw[:plot_samples], alpha0.4, labelPPG Raw) plt.plot(ppg_clean[:plot_samples], labelPPG Clean) # 把峰值点标到图上 valid_peaks [p for p in ppg_peaks if p plot_samples] plt.scatter(valid_peaks, ppg_clean[valid_peaks], colorred, s30, labelDetected Peaks) plt.legend() plt.subplot(2, 1, 2) plt.plot(ecg_raw[:plot_samples], alpha0.4, labelECG Raw) plt.plot(ecg_clean[:plot_samples], labelECG Clean) valid_rpeaks [p for p in ecg_rpeaks if p plot_samples] plt.scatter(valid_rpeaks, ecg_clean[valid_rpeaks], colorred, s30, labelR Peaks) plt.legend() plt.tight_layout() plt.show()这段代码把原始信号、滤波后信号和检测到的峰值点叠在一起画出来。如果发现检测出的峰值点不在波形的视觉最高/最尖位置或者一个周期检测出两个峰说明这段数据污染严重需要预处理后再处理或直接跳过。这一眼值十分钟别省。4. feature.py 特征提取从波形里读出心血管状态4.1 特征体系与生理意义feature.py提取的特征分为三大类我整理一个表格说明每个特征的生理意义有助于理解计算逻辑。特征全称/计算方式生理意义h1主波高度脉搏波第一个正向波峰幅值反映心脏射血能力与动脉弹性h2/h1潮波幅值比潮波反射波幅值相对主波的比例反映动脉血管顺应性和外周阻力h4/h1重搏波幅值比重搏波幅值相对主波的比例反映主动脉瓣功能及血管弹性S波形周期面积一个完整脉搏波周期的总积分整体血流动力学状态S1/S收缩期面积占比收缩期血流灌注效率S2/S舒张期面积占比舒张期冠脉灌注水平PWTT脉搏传导时间ECG R 峰到 PPG 起点的时间差动脉硬化程度的核心指标公式层面潮波是主波之后、重搏波之前那个正向波h2/h1 是医学上评估动脉僵硬度的常用指标重搏波是降支中段的切迹后小波反映主动脉瓣关闭时刻的血液回流。这个项目把 h2/h1、h4/h1、S1/S、S2/S 一起算基本覆盖了经典脉搏波形态分析的常用参数。4.2 峰值位置的二次定位feature.py里如果直接用 NK 库的峰值位置来定位主波会碰到一个精度问题NK 报告的 PPG 峰值点并不一定是主波的最高点可能在相邻几个样本点上浮动。高采样率下这点误差影响不大但采样率低于 250Hz 时特征值可能偏离真实值几个百分点。我目前的处理是拿到峰值点后在邻域内再做一次精细搜索。import numpy as np from scipy.signal import find_peaks def refine_peak(ppg_clean, peak_index, window30): 在峰值索引的邻域内搜索局部最大值返回更精确的峰值位置 window: 搜索半窗口大小单位是样本点 start max(0, peak_index - window) end min(len(ppg_clean), peak_index window) segment ppg_clean[start:end] local_peak_idx np.argmax(segment) return start local_peak_idx # 对所有检测到的峰值做精确定位 refined_peaks [refine_peak(ppg_clean, p) for p in ppg_peaks]逻辑说明NK 库报告的PPG_Peaks已经接近真实峰值但find_peaks类算法返回的是局部极值点精度受信号质量影响。先以 NK 峰值为中心开一个 30 点的窗口在窗内用np.argmax重新定位最大点能得到更精确的主波位置。这个精修步骤对后续 h1、h2/h1 等幅值特征影响明显。参数说明window30在 1000Hz 采样率下代表 30ms 的搜索范围对于脉搏波这种主波持续 100ms 以上的波形是足够的。如果采样率是 500Hzwindow 改成 15250Hz 则改成 8保证搜索窗口在时间尺度上一致。4.3 特征计算的完整实现拿到精确峰值后就可以分段计算波形参数了。关键逻辑是把每个脉搏波周期切出来然后分别计算幅值和面积。def compute_wave_features(ppg_clean, ppg_peaks, ecg_rpeaks, sample_rate): 逐周期计算PPG波形特征 返回一个字典每个特征对应一个数组每个周期一个值 features { h1_list: [], # 主波高度 h2_over_h1: [], # 潮波幅值比 h4_over_h1: [], # 重搏波幅值比 S1_over_S: [], # 收缩面积比 S2_over_S: [], # 舒张面积比 PWTT_list: [] # 脉搏传导时间秒 } # 逐个周期处理 for i in range(len(ppg_peaks) - 1): start ppg_peaks[i] end ppg_peaks[i 1] wave ppg_clean[start:end] # 截取一个完整脉搏波周期 if len(wave) sample_rate // 10: # 周期太短说明检测异常跳过 continue # 主波高度 h1周期内的最大幅值假设信号已去除直流分量 h1 np.max(wave) if h1 0: continue # 潮波 h2主波后第一个局部极大值 # 使用 find_peaks 在放大后的信号上找波峰 peaks, _ find_peaks(wave, distanceint(sample_rate * 0.1)) # 按幅值排序最大的是主波第二大的通常是潮波或重搏波 if len(peaks) 2: peak_heights wave[peaks] # 按位置排序区分潮波和重搏波 main_peak_idx np.argmax(peak_heights) if main_peak_idx 0 and len(peaks) 2: h2_candidate peaks[1] # 潮波在主波之后 h2 wave[h2_candidate] if h2_candidate peaks[0] else 0 else: h2 0 else: h2 0 # 波形周期面积 S使用辛普森积分或梯形积分 S np.trapz(wave, dx1.0 / sample_rate) # 收缩期面积 S1从周期起点到主波位置积分 # 舒张期面积 S2从主波位置到周期末尾积分 main_peak_loc np.argmax(wave) S1 np.trapz(wave[:main_peak_loc], dx1.0 / sample_rate) S2 np.trapz(wave[main_peak_loc:], dx1.0 / sample_rate) features[h1_list].append(h1) features[h2_over_h1].append(h2 / h1 if h1 0 else 0) features[h4_over_h1].append(0) # 重搏波定位比潮波麻烦见下文 features[S1_over_S].append(S1 / S if S 0 else 0) features[S2_over_S].append(S2 / S if S 0 else 0) return features参数说明find_peaks的distance参数控制两个峰的最小样本距离设为 0.1 秒对应的样本数能避免把同一个波上的小抖动当成独立波峰。np.trapz是梯形数值积分把离散样本点近似成连续曲线算面积这里用了dx1.0/sample_rate让面积单位是幅值×秒不同采样率下面积值可比。特别提醒潮波和重搏波的定位是这个项目里最容易出 bug 的地方。潮波是主波降支上的一个小波重搏波在更后面、幅值更小。find_peaks如果参数不当可能把噪声上的小凸起也当成波峰。一个更稳的做法是限制波峰的最小高度为h1 * 0.05低于主波高度 5% 的候选峰直接忽略。4.4 重搏波幅值比的稳健算法上面的代码里我没具体写 h4 的计算纯粹是因为用find_peaks找重搏波翻车概率高。更稳的方法是在下降支的中后段找一个特定形态点重搏波切迹是降支上斜率突变的位置可以用二阶差分找到它。def find_dicrotic_notch(wave, sample_rate, h1): 通过二阶差分定位重搏波切迹 返回切迹位置和重搏波峰值位置 # 主波位置 main_peak np.argmax(wave) if main_peak 0 or main_peak len(wave) - 1: return 0, 0 # 在下降支上主波之后寻找二阶差分最大值点曲率最大处 downslope wave[main_peak:] if len(downslope) sample_rate // 20: return 0, 0 # 二阶差分 second_diff np.diff(np.diff(downslope)) if len(second_diff) 0: return 0, 0 # 切迹是下降支中曲率最大的位置 notch_idx np.argmax(second_diff) main_peak # 重搏波峰值在切迹之后最近的一个局部最大点 post_notch wave[notch_idx:] local_peaks, _ find_peaks(post_notch, heighth1 * 0.05) if len(local_peaks) 0: return notch_idx, 0 dicrotic_peak notch_idx local_peaks[0] return notch_idx, dicrotic_peak切迹和重搏波的关系是固定的切迹是重搏波的起点重搏波峰在切迹之后的局部极大值。先定位切迹再找重搏波峰比在全波形上盲找靠谱得多。实测下来这个逻辑在信号质量中等以上的数据段里都能稳定工作。5. 避坑指南与实践细节七个最容易翻车的场景5.1 峰值检测把 T 波当成了 R 峰现象ECG 的 R 峰位置列表里有些位置明显不对特征值 PWTT 算出来忽大忽小。原因如果原始 ECG 信号 T 波幅度较大或者滤波参数不当导致 QRS 波和 T 波幅度接近自动峰值检测算法会把 T 波也当成一个峰。NK 库虽然做了优化但面对高 T 波形态的个体运动员心脏尤其常见误检并不罕见。解决把 ECG 检测到的峰按 RR 间期合理性过滤一遍正常 RR 间期在 0.6~1.2 秒之间超出这个范围的峰直接丢弃。5.2 采样率填错导致 PWTT 偏差几十毫秒现象算出的 PWTT 值集中在 500ms 以上明显不符合生理范围正常 100~300ms。原因nk.ecg_process和nk.ppg_process里的sampling_rate参数与真实采样率不符。最常见的是把 500Hz 的数据当成 1000Hz 处理导致所有时间特征偏大 2 倍。解决读取后先用data.shape[0] / 采集时长估算采样率确认后再传入处理函数。5.3 基线漂移让面积特征全面失真现象S1/S 和 S2/S 算出来普遍偏大或偏小而且同一个受试者不同时段的结果差异巨大。原因原始信号含有低频基线漂移NK 库的带通滤波虽然能压掉大部分低频分量但当漂移幅度极大时会超过滤波器的抑制范围残留下来的基线会在每个周期里引入正或负的直流偏置。解决在预处理之外再做一次逐周期基线校正就是每个脉搏波周期的起点和终点的幅值应该一致把不一致的部分拉平。def remove_baseline_per_cycle(ppg_clean, ppg_peaks): 逐周期去除基线漂移 corrected ppg_clean.copy() for i in range(len(ppg_peaks) - 1): start ppg_peaks[i] end ppg_peaks[i 1] if end len(corrected): break # 取起点和终点的平均值作为基线值 baseline (corrected[start] corrected[end]) / 2 corrected[start:end] - baseline return corrected5.4 潮波特征对多数数据都算出 0现象h2/h1 大量为 0特征完全没有区分度。原因find_peaks的默认参数要求峰值有足够的高度和突出度老年受试者或血管硬化程度高的个体潮波本身就很微弱高度可能低于主波的 5%被过滤掉了。解决降低潮波检测的阈值要求改用形态学定位——主波之后下降支上第一个明显的曲率变化点而不是严格意义上的局部峰值。5.5 不同采样率的数据混在一起比较特征现象同一批数据里部分文件算出的特征分布明显偏离整体散点图上有孤立簇。原因不同批次数据可能来自不同采集设备采样率不同而特征计算脚本里没有做统一重采样。解决在读取数据后统一重采样到 1000Hz或者在最后算 PWTT 等时间特征时用各自真实的采样率换算成秒再合并而不是用样本点直接比较。5.6 数据长度不足导致周期切分失败现象代码报IndexError或者某段数据计算出来只有一两个周期的特征。原因ppg_peaks[i1]索引越界或者最后一个周期没有完整的长度。解决处理时对range(len(ppg_peaks) - 1)做保护且对周期长度设置下限不足半个脉搏波周期的直接跳过。5.7 运动伪迹段污染整个特征均值现象整段数据的特征均值显著偏离该受试者其他时段的结果。原因采集过程中被试说话、咳嗽或肢体移动导致部分时间段出现大幅伪迹NK 库的滤波无法完全消除伪迹被当成真实波形参与特征计算。解决按滑动窗口计算特征值再用绝对中位差MAD剔除异常窗口。6. 把 PWTT 算准从 R 峰到脉搏波起点的毫秒级对齐PWTT 是这套特征里最有临床区分度的指标它定义是 ECG 的 R 峰到对应脉搏波起始点的时间差。这个特征对时间分辨率极其敏感心率变异性、R 峰定位精度、脉搏波起点定义都会影响最终数值。很多人算不准 PWTT 的根本原因在于把脉搏波的峰值点当成了时间标记点。峰值点容易受反射波叠加影响而左右漂移而起搏点脉搏波上升支的起点更稳定更能代表动脉搏动的起始时刻。计算 PWTT 的正确流程是先对 ECG 和 PPG 各自完成峰值检测 → 对每个 ECG R 峰找它后面最近的 PPG 脉搏波起点 → 两个位置的时间差秒就是 PWTT。核心难点在「后面最近」这四个字——由于心率变异R 峰和下一个脉搏波起点的间隔不是固定值需要逐对匹配而不是简单地取固定偏置。def compute_pwtt(ecg_rpeaks, ppg_peaks_onsets, sample_rate): 计算每个心跳的脉搏传导时间 ecg_rpeaks: ECG R峰位置索引列表 ppg_peaks_onsets: PPG脉搏波起点索引列表不是峰值点 sample_rate: 采样率 返回 PWTT 数组单位秒 pwtt_list [] # 构建便捷索引 onset_arr np.array(ppg_peaks_onsets) rpeak_arr np.array(ecg_rpeaks) for r_peak in rpeak_arr: # 找 r_peak 之后最近的脉搏波起点 future_onsets onset_arr[onset_arr r_peak sample_rate * 0.05] # 至少50ms后 if len(future_onsets) 0: continue nearest_onset future_onsets[0] # 最近的那个 pwtt_samples nearest_onset - r_peak pwtt_sec pwtt_samples / sample_rate # 过滤生理上不合理的值30ms~400ms为合理范围 if 0.03 pwtt_sec 0.40: pwtt_list.append(pwtt_sec) return np.array(pwtt_list)这段代码的关键点在两个过滤条件r_peak sample_rate * 0.05保证取的起点必须晚于 R 峰至少 50ms避免把 R 峰之前的杂音误认为下一个周期的起点0.03 pwtt_sec 0.40把生理范围外的值剔除掉。实测下来加了这两个过滤之后PWTT 的均值稳定性会明显提升。这里有个值得养成的习惯每次计算 PWTT 时我做的第一件事永远是把 R 峰位置和脉搏波起点位置打点画在信号图上一条条对齐检查。这步虽然费时间但坚持做了几组数据之后你对数据质量和算法的可靠性会有非常直观的感觉——哪些被试的 PWTT 就是稳定在 180ms 左右哪些是信号太差算出来乱跳。脉搏传导时间是衡量动脉硬化的一个窗口算准了它前面所有预处理和峰值检测的工作才真正有价值。最后补一句也许是整套流程里最值钱的一句话如果按这篇文章的路径跑通了feature.py建议你把自己数据算出的特征分布保存下来然后用同一批数据对比三组不同参数的结果看看哪些特征对参数不敏感比如 S1/S哪些一旦采样率填错就完全失真比如 PWTT。这种运行习惯会帮你慢慢建立对信号处理系统的直觉也会省掉很多「看起来没问题但总归不放心」的调试时间。希望帮到你。本文还有配套的精品资源点击获取