简介:针对同步采集的PPG、ECG、BP三种生理信号,提供一套完整的MATLAB特征提取方案,适合生物医学工程、健康监测与临床预研场景。程序覆盖信号预处理、时间对齐、特征计算全流程:PPG端采用平滑滤波与峰值检测,提取脉搏周期和波形形态;ECG端实现基线校正、R波定位与PQRS段分析;BP端进行周期切分并计算收缩压、舒张压关键点;同时包含三路信号按参考ECG对齐、频域熵估计、尺度空间关键点提取等高级模块。压缩包共30个文件,以28个m脚本为主,配套一个mat示例数据与README说明,整体大小33MB,函数命名直观便于定位。已有1403人学习下载,代码模块化程度高,可在真实采样数据上直接运行,提取的特征可用于后续机器学习诊断模型。对从事心血管健康分析、无创血压估计或脉搏波研究的开发者,是一套可直接借鉴的算法库;借助这些特征,可进一步开展心率变异性、脉搏波传导时间、动脉硬化风险评估等研究,也可将关键算法迁移到其他生理信号分析中,兼具学习与实用价值。
1. MATLAB 同步信号特征提取:从 PPG、ECG、BP 里拿到你要的指标
做生理信号分析的同行应该都有同感:设备把 PPG、ECG、BP 三路信号同步采回来只是第一步,真正的活儿全在后面的特征提取上。你盯着波形看半天,不如让 MATLAB 把心率、脉搏传导时间、血压趋势这些关键指标一次性算出来。这份 bp_features_extraction 程序就是干这个事的——它把同步采集的 PPG(光电容积脉搏波)、ECG(心电)、BP(血压)三路信号统一处理,提取出可用于后续分析、建模或实验报告的特征参数。适合正在做心血管功能评估、睡眠分期、运动疲劳监测或者可穿戴设备算法验证的开发者,尤其是手里已经有采集好的数据、但不知道怎么下手做特征工程的人。
程序本身不复杂,但涉及信号预处理、峰值检测、特征点定位和参数计算四个环节,每一环都有不少容易翻车的细节。下面我按照自己拆解的思路,把这套程序的核心逻辑和实操要点完整过一遍。
2. 三路信号的特征提取逻辑:为什么是 PPG、ECG、BP 一起算
2.1 三路信号各自能提供什么特征
先搞清楚程序的输入输出关系。PPG 信号来自光电传感器,反映的是外周血管容积随心脏搏动的变化,主要特征包括脉搏周期、波峰幅值、上升沿时间和反射波位置;ECG 信号提供心电活动的准确时相标记,R 波峰值点是最关键的参考点,可以算出心率变异性(HRV)和 RR 间期;BP 信号则是连续血压波形,特征集中在收缩压、舒张压、平均动脉压和脉压差上。
把三路信号放在一起分析的价值在于时间对齐。PPG 的峰值点和 ECG 的 R 波之间存在一个时间延迟,这个延迟就是脉搏传导时间(PTT),它和血压变化有很强的相关性,是很多无袖带血压估计算法的核心特征。程序必须在同一时间轴上完成三路信号的特征点定位,否则后续计算出来的 PTT、HRV 这些指标全是错的。
2.2 程序整体处理流程
这个 MATLAB 程序的处理流程可以拆成四个阶段:
- 数据读取与时间对齐:把三路信号载入工作区,统一采样率,必要时做插值对齐
- 信号预处理:去基线漂移、滤除工频干扰和高频噪声
- 特征点检测:ECG 找 R 波,PPG 找峰值和起始点,BP 找收缩压和舒张压对应位置
- 特征参数计算:根据特征点位置计算心率、PTT、血压均值和变异系数
预处理环节是决定特征提取质量的关键。常见做法是先用移动平均或者高通滤波器去掉基线漂移,再用带通滤波器保留有效频段——ECG 一般用 0.5-45 Hz,PPG 用 0.5-10 Hz,BP 用 0.5-30 Hz。滤波器阶数和截止频率的选择直接影响后面峰值检测的准确率,这里没有统一答案,需要根据实际采集设备的噪声情况调整。
2.3 时间同步问题的处理
我见过不少人在多路信号分析时踩过同一个坑:三路数据明明是同步采集的,但因为在预处理阶段分别做了滤波和重采样,导致信号长度不一致或者时间偏移,最后算出来的特征参数出现系统性偏差。
程序里处理这个问题的常规做法是:
% 假设 fs 是原始采样率,target_fs 是统一采样率 if fs ~= target_fs % 对三路信号分别重采样到统一采样率 ecg_resampled = resample(ecg_raw, target_fs, fs); ppg_resampled = resample(ppg_raw, target_fs, fs); bp_resampled = resample(bp_raw, target_fs, fs); endresample函数会自动设计抗混叠滤波器,比直接用interp1线性插值更可靠。重采样之后务必检查三路信号的长度是否一致,不一致的话要在后续特征提取前处理掉。
2.4 滑动窗口的设计
程序里做特征提取时通常不是对整段数据一次性计算,而是用滑动窗口按时间段切分。窗口长度一般取 30 秒到 5 分钟不等,具体看你要分析的是短期 HRV 还是长期趋势。窗口重叠率通常在 50% 左右,这样可以减少窗边界处的特征突变。
窗口设计对特征序列的平稳性影响很大,窗口太短会导致心率这类低频特征波动剧烈,窗口太长又会把瞬时变化抹平。我一般会先把数据画出来看几段典型波形,确定合理的窗口长度再批量跑。
3. 核心特征点检测:峰值定位算法与参数调优
3.1 ECG R 波检测的经典方法
ECG 信号里最重要的特征点是 R 波峰值。程序里通常使用的检测方法基于 Pam-Tompkins 算法思路:先对 ECG 做带通滤波,然后求差分、平方、移动窗口积分,最后用自适应阈值找峰值。
% ECG R 波检测核心流程 % 1. 带通滤波(5-15 Hz 突出 QRS 波群能量) [b, a] = butter(2, [5 15] / (fs/2), 'bandpass'); ecg_filtered = filtfilt(b, a, ecg_resampled); % 2. 差分 + 平方 + 滑动积分 diff_ecg = diff(ecg_filtered); squared_ecg = diff_ecg .^ 2; window_len = round(0.12 * fs); % 150ms 积分窗口 integrated_ecg = conv(squared_ecg, ones(1, window_len) / window_len, 'same'); % 3. 自适应阈值检测 R 波峰值位置 threshold = 0.3 * max(integrated_ecg); [~, r_peaks] = findpeaks(integrated_ecg, ... 'MinPeakHeight', threshold, ... 'MinPeakDistance', round(0.25 * fs));注意filtfilt是零相位滤波,不会引入相位偏移,这一点对同时计算 PTT 至关重要。如果用普通的filter函数,滤波造成的群延迟会直接污染 R 波和 PPG 峰值之间的时间差计算。MinPeakDistance参数需要根据实际心率范围设置——正常成年人静息心率 60-100 次/分,对应的 R-R 间期是 0.6-1.0 秒,所以最小峰值间距设 0.25 秒既能防止误检又能保证不遗漏。
3.2 PPG 峰值与起始点检测
PPG 信号的特征提取比 ECG 要麻烦一些,因为它的波形形态受外周血管弹性和血流量影响,同一个人的不同时刻波形差异也很大。程序里通常需要提取两个关键点:主波峰值(收缩期峰值)和波形起始点(舒张期最低点)。
% PPG 峰值检测(在主波频段内) [b_ppg, a_ppg] = butter(2, [0.5 8] / (fs/2), 'bandpass'); ppg_filtered = filtfilt(b_ppg, a_ppg, ppg_resampled); % 用 findpeaks 找主波峰值 [ppg_peaks, ppg_locs] = findpeaks(ppg_filtered, ... 'MinPeakHeight', mean(ppg_filtered) + std(ppg_filtered), ... 'MinPeakDistance', round(0.3 * fs));MinPeakHeight这里用均值加标准差作为阈值,是一个比较通用的自适应策略。如果你的 PPG 信号基线漂移比较厉害,可以先用hampel滤波或者中值滤波做一遍基线校正再检测峰值。PPG 起始点检测没有统一算法,常见做法是取相邻两个主波峰值之间的局部最小值,或者用一阶导数的最小值点来定位。
3.3 BP 收缩压和舒张压的定位
BP 连续波形上的收缩压对应波形峰值,舒张压对应波形谷值。程序里的做法通常和 PPG 峰值检测类似,但需要注意 BP 波形的形态在不同测量场景下差异很大——动脉僵硬的患者收缩期峰值和重搏波可能糊在一起,机械通气患者的胸腔压力变化也会影响波形形态。
% BP 特征点:收缩压(峰值)+ 舒张压(谷值) [bp_peaks, bp_peaks_loc] = findpeaks(bp_resampled, ... 'MinPeakHeight', mean(bp_resampled) + 0.5 * std(bp_resampled), ... 'MinPeakDistance', round(0.3 * fs)); % 谷值检测:对信号取反后找峰值 [bp_troughs, bp_troughs_loc] = findpeaks(-bp_resampled, ... 'MinPeakHeight', abs(mean(bp_resampled)) * 0.5, ... 'MinPeakDistance', round(0.3 * fs));收缩压数值就是bp_peaks对应位置的幅值,舒张压数值是bp_troughs对应位置的幅值,平均动脉压可以用积分均值或者经验公式MAP = (SBP + 2 * DBP) / 3来估算。
3.4 特征点质量校验
特征点检测完不能直接拿去计算,要先做一轮合理性校验。常见做法是检查相邻特征点的时间间隔是否落在生理范围内——比如 RR 间期小于 0.3 秒或大于 2 秒,基本可以断定是误检或漏检。程序里可以用一组简单的阈值规则过滤:
% RR 间期合理性检查 rr_intervals = diff(r_peaks) / fs; valid_rr = rr_intervals > 0.3 & rr_intervals < 2.0; r_peaks_clean = r_peaks([true; valid_rr]); % 保留有效点如果异常点比例超过 5%,说明前面的滤波或阈值参数有问题,需要回头调整,而不是硬着头皮用这些疑似异常的特征点继续算参数。这一步我一般会单独做一个检查模块,把异常点位置标记出来画图确认。
4. 特征参数计算:从特征点到可用指标
4.1 心率与 HRV 指标
R 波位置确定后,心率计算是直接取 RR 间期平均值的倒数换算成每分钟搏动次数:
% 心率计算 rr_intervals = diff(r_peaks) / fs; % 单位:秒 heart_rate = 60 / mean(rr_intervals); % 单位:次/分 % HRV 时域指标 sdnn = std(rr_intervals); % 正常 RR 间期标准差 rmssd = sqrt(mean(diff(rr_intervals).^2)); % 相邻差值均方根SDNN 反映整体心率变异性,RMSSD 主要反映副交感神经活性,这两个指标在疲劳评估和自主神经功能研究里用得比较多。如果程序里有频域 HRV 分析需求,需要把 RR 间期序列用插值重采样成等间隔序列再做 FFT,不能在原始不等间隔数据上直接算功率谱,否则频谱会失真。
4.2 脉搏传导时间 PTT 的计算
PTT 是 ECG R 波峰值到 PPG 主波峰值的时间差,这个参数是无袖带血压估算的核心特征。计算时要注意两点:一是 R 波和 PPG 峰值必须在同一个心跳周期内一一对应,二是滤波产生的相位偏移要最小化——这也是前面强调必须用filtfilt的原因。
% PTT 计算:R 波峰值到 PPG 主波峰值的延迟 % 先保证两个序列长度一致且按时间排序 ptt_values = (ppg_locs_interp - r_peaks_interp) / fs * 1000; % 单位:ms % 去除生理上不可能的 PTT(通常 50-400 ms 之间) valid_ptt = ptt_values > 50 & ptt_values < 400; ptt_mean = mean(ptt_values(valid_ptt));如果程序报出奇异的 PTT 值(比如负值或者几千毫秒),多半是 R 波和 PPG 峰值没有配对正确。常见做法是以 R 波位置为参考,在后面的 100-500 ms 窗口内找对应的 PPG 峰值,而不是全局找最近的峰值点。
4.3 血压相关特征
BP 信号的特征除了收缩压、舒张压和平均动脉压之外,还可以计算脉压差(PP = SBP - DBP)和血压变异系数。连续血压波形中的逐搏变异是评估血压稳定性的重要指标,程序里可以这样实现:
% 逐搏血压参数(假设每个心跳周期对应一个 BP 波峰) sbp_values = bp_peaks; % 收缩压序列 dbp_values = -bp_troughs; % 舒张压序列(取反恢复原值) % 如果 BP 信号单位不是 mmHg,需要先做标定转换 pp_values = sbp_values - dbp_values; % 脉压差 map_values = (sbp_values + 2 * dbp_values) / 3; % 平均动脉压4.4 特征参数汇总与导出
程序把所有提取到的特征汇总成表格输出,通常每一行对应一个时间窗口,每一列是一种特征参数。MATLAB 里的输出方式可以是table数据结构,直接导出 CSV 供后续分析:
% 汇总特征表 feature_table = table(... time_window_start, ... heart_rate, ... sdnn, ... rmssd, ... ptt_mean, ... sbp_mean, ... dbp_mean, ... map_mean, ... 'VariableNames', {'TimeStart', 'HR', 'SDNN', 'RMSSD', 'PTT', 'SBP', 'DBP', 'MAP'}); % 导出 CSV writetable(feature_table, 'features_output.csv');导出之前要检查一下数值的范围和缺失值。特征序列里偶尔会冒出NaN或者Inf,多半是某段信号质量太差导致特征点检测失败,这需要在代码里显式处理,比如用上一窗口的有效值填充或者直接标记为缺失。
5. 避坑与常见问题:信号处理里的几个大坑
5.1 滤波导致特征点偏移
现象:计算的 PTT 值系统性偏大或偏小,与文献值差距明显。
原因:用了filter而不是filtfilt,滤波器带来的群延迟污染了特征点的时间标注。或者是带通滤波器的截止频率设置不当,把特征点的频率成分削弱了。
解决:统一改用filtfilt做零相位滤波。PPG 的滤波通带不要压得太窄,建议保留 0.5-8 Hz,ECG 保留 5-20 Hz 的 QRS 主能量区间。如果还是偏,就把滤波前后的信号叠在一起画图对比峰值位置差异。
5.2 三路信号长度不一致
现象:r_peaks和ppg_locs数量差很多,或者下标对不上导致索引越界。
原因:预处理时分别对三路信号做了滤波和重采样,MATLAB 的filter输出长度与输入相同但resample可能因为插值边界产生长度差异。更常见的是原始数据里本来就存在一个或多个通道首尾有裁剪。
解决:在预处理完统一检查length(ecg) == length(ppg) == length(bp),不一致就用min_length截断统一长度。每次加载新数据都强制跑一次这个检查,不要想当然。
5.3 高心率受试者的峰值漏检
现象:心率超过 120 次/分时,R 波检测数量明显少于实际搏动次数。
原因:MinPeakDistance参数设置过大,比如设成了 0.4 秒,高心率时实际 RR 间期只有 0.4 秒左右,部分峰值被程序当成相邻峰的旁瓣丢弃。
解决:MinPeakDistance设为期待最小 RR 间期的 70%。比如最高预期心率 180 次/分,对应 RR 间期 0.333 秒,MinPeakDistance设为 0.25 秒比较安全。
5.4 BP 信号零漂和饱和
现象:BP 波形整体数值漂移,或者峰值削顶,提取的收缩压明显偏离真实值。
原因:连续血压采集设备长时间记录后会有基线漂移,或者传感器贴附压力变化导致波形饱和削顶。
解决:先做基线校正(中值滤波或高通滤波),削顶段如果超过总时长的 10%,这段数据直接放弃不要用来算血压特征,硬算只会污染整批结果。
5.5 特征点配对错位
现象:PTT 值有一半是负的,一半是正的,看起来没规律。
原因:ECG 的某个 T 波被误检成 R 波,或者 PPG 的反射波被误检为主波峰值,导致配对到相邻周期。
解决:先用模态滤波检查 R 波间期是否落在 0.3-2.0 秒范围内,PPG 主波峰值要在 R 波后 50-400 毫秒窗口内匹配。加一个配对逻辑:每个 R 波向下一个最近的 PPG 峰值匹配,但限定搜索窗口。
6. 进阶用法:自动批处理与特征质量评估
程序跑通单条数据之后,真正的价值在于批量处理。我会把整个特征提取流程封装成一个函数,输入是文件路径和参数结构体,输出是特征表格和每段数据的质量评分。这样几十条实验数据可以一口气跑完,不用一遍遍手工操作。
批处理的核心是数据质量评估。我给每段数据打三个维度的分数:波形质量(信噪比估算)、特征点检出率(有效 R 波数占预期搏动数的比例)、特征合理性(心率、PTT 是否落在生理范围内)。综合分数低于阈值的片段自动打标,不进入后续统计。这个做法帮我在处理可穿戴设备长时程数据时省了大量人工翻阅波形的时间,也让实验数据里混入的坏片段没法浑水摸鱼。
另外还要结合项目场景做特征扩展。如果是做睡眠分期,相邻窗口的 PTT 变化率和 HRV 频域特征更重要;如果是做运动疲劳监测,心率恢复时间和血压反应性更有价值。程序本身只给了基础特征提取能力,但它的模块化结构让你可以很方便地在输出表格后面追加新特征列。
用这套程序跑了三批实验数据之后,我才意识到信号质量检查比特征提取算法本身更影响最终结果。从那以后我每次处理新数据都强制先跑一遍质量评估模块,坏数据一律不放过。波形质量不过关的片段,提取出来的特征值再漂亮也不能用——这是我从翻车经历里得出的教训,希望帮到你。
本文还有配套的精品资源,点击获取