简介面向医疗健康领域研究人员与MATLAB开发者用于对同步采集的PPG、ECG、BP三路生理信号完成滤波、对齐与特征提取可服务于疾病诊断、健康监测等场景。压缩包共30个文件包含28个m脚本、1个mat样例数据与1个README说明脚本覆盖信号预处理、峰值/关键点检测、跨信号对齐、频域特征计算等功能mat文件可供直接运行验证。资源已有1403人学习代码结构清晰支持从原始波形到时域、频域及非线性特征的一站式提取并保留ECG参考对齐、PPG峰值检测等模块便于替换参数或扩展至机器学习输入。整体约33.31MB适合有一定信号处理基础、需要快速搭建生理信号特征提取流程的研究者参考使用。1. bp_features_extraction到底解决什么问题PPG、ECG、BP三路信号一起提特征很多人在搜索引擎里看到 bp_features_extraction 这个 MATLAB 程序名第一反应以为是反向传播Backpropagation相关的工具箱实际上这里的 BP 是 Blood Pressure血压。这个程序的核心任务很直接把同步采集到的 PPG、ECG、BP 三路信号读进来统一做预处理然后分别从每一路里提取生理特征——PPG 提供脉率和灌注信息ECG 提供精确的心搏定位BP 提供逐拍的收缩压、舒张压和趋势。三路信号联合在一起最典型的特征是脉搏传导时间PTT也就是从心脏搏动到外周脉搏波到达之间的时间差它是血压评估、动脉硬化风险筛查里非常常用的指标。这篇文章按“数据读取 → 预处理 → 三路各自提特征 → 联合特征 → 踩坑记录 → 批处理验证”的顺序把整个流程拆到可以直接复现的代码层面。2. 数据读取与预处理三路信号同步性检查与滤波参数选择2.1 readtable 读入三路信号后先做长度与采样率自检拿到一个采集文件之后第一件事不是急着提特征而是回答三个问题采样率是多少、三路信号长度是否一致、时间戳是不是等间隔。采集卡标称的采样率经常和实际有偏差或者某个通道中途丢了几十个点这些都会在后面的 R 波定位、PPG 峰值检测、PTT 计算里变成系统性误差。我一般会让程序先做一次自检把风险在入口处截住。% 读入三路信号时间戳、PPG、ECG、BP 四列 data readtable(record_01.csv); fs 250; % 采集卡标称采样率先按文档填后面会核实 t data{:, 1}; % 第一列是时间戳秒 ppg data{:, 2}; % 第二列是 PPG ecg data{:, 3}; % 第三列是 ECG bp data{:, 4}; % 第四列是 BP % 实际采样率校验时间戳差分取中位数抗单点抖动 dt diff(t); fs_actual 1 / median(dt); fprintf(标称采样率 %.0f Hz实测 %.2f Hz\n, fs, fs_actual);readtable 把 CSV 读成 table用花括号拿到数值矩阵这样比csvread更稳表头带不带都能处理。采样率校验这里的技巧是用median(dt)而不是mean(dt)个别时间戳抖动会让均值偏移但中位数不受少数离群点影响。如果实测采样率和标称值差超过 1%后面所有基于时间的计算都要用fs_actual否则 PTT 会出现几毫秒级别的固定偏移对血压相关分析来说是致命的。长度自检也很关键尤其是多通道采集卡从不同模数转换器出来偶尔会有长度不一致的情况。如果在初始化阶段就统一截断省得后续每个特征提取函数都要处理索引越界。% 长度自检与截断三路必须等长否则逐点运算直接错位 if length(ecg) ~ length(ppg) || length(ecg) ~ length(bp) n min([length(ecg), length(ppg), length(bp)]); ecg ecg(1:n); ppg ppg(1:n); bp bp(1:n); warning(三路长度不一致已统一截断到 %d 点, n); end这里的思路是最小长度对齐取三路最短的公共长度把超出的部分丢掉。注意这种处理只能解决“截断”不能解决“时间偏移”——如果三路信号启动时刻不同单纯截断后 R 波和 PPG 峰在时间轴上仍然是对不齐的这个问题在第 5 章会专门讲。如果只想快速验证流程截断就够了如果要做严肃的 PTT 分析必须把时间偏移估计加进来。2.2 用 filtfilt 做零相位滤波不同信号的频率边界怎么定同步信号预处理里最容易犯的错是用了filter而不是filtfilt。filter是有相位延迟的延迟量随频率变化这会让 R 波位置、PPG 峰值位置发生几个毫秒到几十毫秒的偏移直接毁掉 PTT。MATLAB 里首选filtfilt做零相位滤波虽然计算量稍大但波形特征点位置不会漂移。三路信号频率成分差异很大滤波参数要分开定信号推荐带通范围理由ECG0.5–45 Hz保留 QRS 主能量压低 T 波与肌电干扰PPG0.1–10 Hz保住上升沿与重搏波切迹去掉高频噪声BP0.3–30 Hz保留压力均值与波形形态滤掉机械伪迹ECG 的 0.5 Hz 高通用于压制呼吸引起的基线漂移45 Hz 上限在保留 QRS 高频成分的同时削弱肌电。PPG 的低通取 10 Hz 比较合适再低会把重搏波切迹抹掉。BP 信号比较特殊它本身有真实的直流分量高通下限不能太高0.3 Hz 以上会切掉呼吸性血压波动。% 零相位带通滤波butter 二阶 filtfilt 消除相位偏移 [b1, a1] butter(2, [ 0.5 45] / (fs/2), bandpass); % ECG [b2, a2] butter(2, [ 0.1 10] / (fs/2), bandpass); % PPG [b3, a3] butter(2, [ 0.3 30] / (fs/2), bandpass); % BP ecg_f filtfilt(b1, a1, ecg); ppg_f filtfilt(b2, a2, ppg); bp_f filtfilt(b3, a3, bp);butter的归一化频率是除以奈奎斯特频率fs/2所以采样率变了之后这些参数必须重新计算不能把 250 Hz 采样率下的参数直接套到 500 Hz 数据上。二阶滤波器对于生理信号足够阶数太高会在通带边缘引入振铃反而让 QRS 波群尾部出现假的凹陷。filtfilt是零相位滤波代价是边缘效应信号开头和结尾几十个点会失真尤其是高频成分较少的 PPG。我的习惯是滤波完成之后把每端 0.5 秒的数据标记为“不可用区”在特征统计阶段排除掉这个细节后面还会碰到。3. ECG 特征提取R 波检测与 HRV 时域参数计算3.1 用 findpeaks 做 R 波检测从候选峰值到真实 QRSECG 特征提取的核心是 R 波定位。最直接的思路是对滤波后的信号直接findpeaks但这在真实数据上很容易翻车T 波在某些导联上幅度接近 R 波的 1/3运动伪迹也会制造假峰。常见做法是先对信号做一阶差分平方把 QRS 波群的陡峭能量突出来再用滑动窗积分最后在积分信号上找峰。% R 波检测差分平方 滑动窗积分 阈值找峰 d diff(ecg_f); s d .^ 2; win round(0.12 * fs); % 120ms 滑动窗对应 QRS 典型宽度 m movmean(s, win); thr 0.3 * max(m); % 初始阈值取峰值能量的 30% % 在积分信号上找候选峰250ms 不应期抑制 T 波 [~, locs0] findpeaks(m, ... MinPeakHeight, thr, ... MinPeakDistance, round(0.25 * fs));差分平方会把 QRS 的高斜率特征放大而 P 波、T 波斜率低能量被压缩到一个很小的量级。movmean的作用是把 QRS 能量在 120ms 窗口内累积起来形成一个个突出的峰这样阈值检测更稳定。MinPeakDistance设为 0.25 秒对应最大 240 次/分的心率低于这个间隔出现的候选峰直接丢弃。积分信号上找到的峰位置和真实 R 波位置会有一点点偏移因为滑动窗积分本身是对窗口内能量求和峰值位置偏在窗口末尾一侧。我一般会做一次精修在候选点附近 ±20ms 的原始滤波信号窗口内找最大值把 R 波位置拉回真实波峰。% 在原始滤波信号上精修 R 波位置候选点附近 ±20ms 找最大幅值 fs 250; locs zeros(size(locs0)); for k 1:length(locs0) w max(1, locs0(k)-round(0.02*fs)) : min(length(ecg_f), locs0(k)round(0.02*fs)); [~, j] max(ecg_f(w)); locs(k) w(1) j - 1; end locs locs / fs * 1000; % 转为毫秒时间戳精修这一步看着不起眼但对后续 PTT 计算的影响很大。20ms 的定位误差在 100ms 量级的 PTT 里占比接近 20%不可忽略。locs转成毫秒是因为后面 HRV 和 PTT 统一用毫秒做单位避免混用。3.2 RR 间期与 HRV 参数SDNN、RMSSD、pNN50 的计算口径R 波位置确定之后RR 间期直接做差分就能拿到。HRV 时域参数有三个是最常用的SDNN 反映整体变异性RMSSD 反映短期迷走神经调节pNN50 是相邻 RR 间期差值超过 50ms 的比例。这三个参数的计算口径有讲究尤其是离群值处理如果某个运动伪迹让相邻两个 R 波被误检成一个RR 间期会直接翻倍SDNN 被拉爆。% RR 间期 - HRV 时域参数 rr_ms diff(locs); % 单位毫秒 % 离群值修复中位数滤波替换超过 30% 偏差的间期 med_rr movmedian(rr_ms, 5); bad abs(rr_ms - med_rr) 0.3 * med_rr; rr_ms(bad) med_rr(bad); sdnn std(rr_ms); % 全部 RR 间期标准差 rmssd sqrt(mean(diff(rr_ms).^2)); % 相邻差值均方根 pnn50 sum(abs(diff(rr_ms)) 50) / (length(rr_ms)-1) * 100;diff(rr_ms)会让长度再减一pNN50 的分母用length(rr_ms)-1才准确这个细节经常被忽略。离群值修复用的是 5 点的移动中位数把每个 RR 间期和邻近中位数比较偏差超过 30% 就替换。这个阈值不能太激进否则会把真实的心率突跳也抹掉也不能太宽松否则起不到修复作用。0.3 倍这个数值是经验值我在实际数据上试过对运动伪迹和漏检都有不错的鲁棒性。4. PPG 与 BP 特征提取峰值定位、血压极值与脉搏传导时间4.1 PPG 峰值检测与形态特征MinPeakProminence 比 MinPeakHeight 更可靠PPG 波形比 ECG 平滑得多直接用MinPeakHeight定阈值找峰很容易出问题基线会随呼吸和运动缓慢漂移固定高度阈值在信号偏弱的段落里会把有效峰漏掉。更可靠的做法是用MinPeakProminence突出度衡量的是峰相对两侧谷底的显著程度天然免疫基线漂移。% PPG 峰值检测突出度阈值抗基线漂移 [~, locs_ppg] findpeaks(ppg_f, ... MinPeakProminence, 0.2 * std(ppg_f), ... MinPeakDistance, round(0.3 * fs)); % 至少 0.3 秒一拍 ppg_peak_vals ppg_f(locs_ppg); % 峰幅值 ppg_amp ppg_peak_vals - mean(ppg_f); % 相对基线的幅度 ppi_ms diff(locs_ppg) / fs * 1000; % 脉率间期突出度阈值取0.2 * std(ppg_f)是一个经验起点信号质量好时可以降到 0.15运动伪迹严重时提到 0.3。MinPeakDistance用 0.3 秒对应 200 次/分的极限脉率正常情况下足够。重搏波切迹是 PGG 峰值检测最大的敌人年轻受试者身上切迹明显幅度可达主峰的 40% 左右如果最小间距设得太小切迹会被当成主峰。把MinPeakDistance设在 0.3 秒以上配合突出度阈值能过滤掉绝大多数切迹误检。PPG 的峰-峰间期PPI理论上和 ECG 的 RR 间期一致两者可以互相校验。如果同一段数据里 PPI 的平均值和 RR 间期差超过 5%优先怀疑 PPG 峰值检测有切迹混入或漏检而不是采集问题。这个交叉验证的习惯能省下大量排查时间。4.2 BP 逐拍极值与联合特征从 R 波到 PPG 峰的时间差血压波形的特征提取和 PPG 不一样。收缩压是每拍的峰值舒张压是波谷但问题在于波形形态多变有的段有明显的重搏波有的段波峰平缓。直接对血压信号独立找峰找谷容易出现收缩压和舒张压数量对不上的情况。我一般以收缩压峰为锚点在两个相邻收缩压之间找最小值作为该拍的舒张压这样逐拍配对不会乱。% 血压波形收缩压找峰 [~, idx_sbp] findpeaks(bp_f, ... MinPeakProminence, 0.15 * std(bp_f), ... MinPeakDistance, round(0.3 * fs)); % 以收缩压峰为锚在两峰之间找谷值作为该拍舒张压 dbp_v zeros(length(idx_sbp)-1, 1); sbp_v zeros(length(idx_sbp)-1, 1); for k 1:length(idx_sbp)-1 seg bp_f(idx_sbp(k):idx_sbp(k1)); [~, j] min(seg); dbp_v(k) bp_f(idx_sbp(k) j - 1); sbp_v(k) bp_f(idx_sbp(k)); end map dbp_v (sbp_v - dbp_v) / 3; % 平均动脉压MinPeakProminence设为0.15 * std(bp_f)对收缩压和舒张压的幅度差异比较敏感。这里的idx_sbp(k)到idx_sbp(k1)是一个心拍周期最后一个收缩压峰之后没有完整周期所以循环只到length(idx_sbp)-1。如果用独立找谷的办法收缩压峰和舒张压谷的配对要靠最近邻查找代码更罗嗦还容易错位。三路信号联合的核心计算是脉搏传导时间 PTT定义是同一心拍里 R 波峰值到外周 PPG 峰值的时间差。PTT 和血压呈负相关血管弹性越好、血压越低脉搏波传导越快。计算的难点在于 ECG 的 R 波和 PPG 峰必须一一对应。% 脉搏传导时间每个 R 波之后找最近的 PPG 峰计算时间差 ptt_ms nan(length(locs), 1); for k 1:length(locs) idx find(locs_ppg locs(k), 1, first); if ~isempty(idx) ptt_ms(k) (locs_ppg(idx) - locs(k)); % 已经是毫秒 end end ptt_ms ptt_ms(~isnan(ptt_ms));locs_ppg和locs都已经是毫秒时间戳直接相减。用NaN标记未匹配的拍后续统计时过滤掉而不是用 0 去占位因为 0 会污染均值。匹配不到的常见原因有两个PPG 漏检或者这段信号正好处在滤波边缘区。PTT 的输出不只是一个均值它的标准差和逐拍变化趋势也很有价值——PTT 变异性增大通常意味着呼吸对血压调节的影响变强这在睡眠分期和压力评估里是常用特征。5. 常见问题与排查三个信号特征提取翻车的五个具体坑5.1 时序与同步类坑截断错位导致 PTT 出现负值三路信号长度不一致截断后 R 波和 PPG 峰仍然对不上PTT 匹配率极低。现象是ptt_ms里大量负值或者 NaN画图看 PPG 峰明显出现在 R 波之前。原因是采集卡的 PPG 通道和 ECG 通道启动时间不同步某一路延迟了几百毫秒单纯截断长度无法修正时间偏移。解决方法是先用互相关估计偏移量再对 PPG 做循环移位% 互相关估计两路信号偏移ECG 为基准检测 PPG 延迟 [c, lags] xcorr(ppg_f, ecg_f, round(1 * fs), coeff); [~, im] max(c); shift lags(im); % 正数表示 PPG 晚于 ECG ppg_f circshift(ppg_f, shift);xcorr的第三个参数限制最大搜索范围是 1 秒防止把相邻心拍的相关性误判成全局最优。shift的单位是采样点正值代表 PPG 信号整体滞后circshift把它向前搬。要注意circshift是循环移位会把末尾的数据搬到开头移位后最好把两端各 1 秒的数据标记为无效。滤波边缘效应让文件开头和结尾的特征异常。现象是前几十个点的 R 波幅度骤降、PTT 明显偏大中间部分正常。原因是filtfilt的边界延拓在信号两端制造了假波形。解决方法是提取特征后统一丢弃每端 0.5 秒数据或者在滤波前把信号前后各延长 1 秒计算完再截掉。% 丢弃边缘效应区域滤波后 0.5 秒以内的特征点全部剔除 edge_samples round(0.5 * fs); valid locs edge_samples locs length(ecg_f) - edge_samples; locs locs(valid);这里的阈值 0.5 秒不是拍脑袋定的二阶 butter 滤波器 filtfilt的边界过渡带大约就是 0.3 到 0.5 秒取 0.5 秒留了余量。5.2 检测与形态类坑T 波误检、切迹误判与血压基线漂移ECG 特征提取里最经典的翻车场景是 T 波被当成 R 波。现象是 RR 间期序列里出现规律性的短-长交替短的那个只有正常间期的一半左右。原因是 T 波在某些受试者身上幅度接近 R 波的 1/3阈值取0.3 * max时把关不严高尖 T 波混进候选峰。解决方法是把 ECG 滤波上限从 45 Hz 降到 35 HzT 波能量被压掉一部分同时把不应期从 250ms 放宽到 300ms在不应期内同时出现的两个候选峰只保留幅度大的那个。% 不应期内保留能量更高的峰抑制 T 波 keep true(size(locs0)); for k 2:length(locs0) if locs0(k) - locs0(k-1) round(0.30 * fs) if m(locs0(k)) m(locs0(k-1)) keep(k-1) false; else keep(k) false; end end end locs0 locs0(keep);这段逻辑的关键是“比较能量而不是比较幅度”积分信号m反映的是 QRS 整体能量比单点幅度更可靠。T 波虽然幅度可能接近 R 波但它的持续时间和斜率都比 QRS 小积分能量只有 QRS 的一半左右。PPG 峰值检测把重搏波切迹当主峰导致 PPI 出现双峰交替。现象是 PPI 序列里一半数值明显偏短画出来每个主峰后面紧跟着一个小峰。原因是切迹在部分受试者身上幅度较大MinPeakProminence阈值不够高。解决方法是把突出度阈值从0.2 * std提高到0.25 * std并把低通上限从 10 Hz 降到 7 Hz 压掉切迹的高频分量。切迹在频域上和主峰没有清晰边界7 Hz 是验证过比较稳妥的折中。血压基线漂移导致舒张压均值偏低。现象是算出来的平均舒张压比监护仪显示值低 5-8 mmHg而且整段血压波形随着呼吸周期缓慢上下摆动。原因是高通滤波下限设到了 1 Hz把呼吸性血压波动和真正的低频成分当成噪声滤掉了但这些低频成分本来就是血压的一部分。解决方法是把 BP 高通下限降回 0.3 Hz提取完逐拍收缩压和舒张压之后再对 SBP/DBP 序列做 3 阶中值滤波去毛刺而不是在波形层面过度滤波。% 对逐拍血压序列做中值滤波保留真实趋势 sbp_v medfilt1(sbp_v, 3); dbp_v medfilt1(dbp_v, 3);注意这里的中值滤波只作用于逐拍序列不作用于原始波形这样才能在去毛刺的同时不破坏血压波形的真实形态。6. 进阶做法批量处理、匹配率验证与结果落盘6.1 批量处理把单文件函数跑遍整个目录单文件流程跑通之后面对几十个record_*.csv手点是不现实的。把特征提取封装成函数用循环批量处理每跑完一个文件就把特征向量追加到结果表格里最后一次性落盘。% 批量处理目录下所有 record_*.csv 都跑一遍特征提取 files dir(record_*.csv); resultTable table(); for i 1:numel(files) out extract_features(files(i).name); % 封装好的单文件入口 resultTable [resultTable; out]; %#okAGROW end writetable(resultTable, all_features.csv);extract_features是前面所有步骤的封装输入文件路径输出一行特征包括心率均值、SDNN、RMSSD、PTT 均值、PTT 标准差、SBP 均值、DBP 均值等。这里用普通for而不是parfor是有原因的一是readtable和writetable在并行池里有文件句柄冲突风险二是单文件处理耗时不长并行加速不明显。如果文件数量上百可以先改成parfor但函数内部绝对不能有打印语句或全局变量。6.2 用匹配率与 SQI 做特征质量验证批量跑完不能直接信结果还要做一次质量验证。最实用的指标是“峰匹配率”PPG 峰落在 R 波之后 0.15 到 0.45 秒窗口内的比例。这个时间窗覆盖了正常 PTT 的生理范围匹配率低于 0.9 的文件要么是 PPG 噪声太大要么是两路信号同步有问题。% 峰匹配率PPG 峰是否落在 R 波后 0.15~0.45s 窗口内 win_lo 0.15 * fs; win_hi 0.45 * fs; matchCount 0; for k 1:length(locs_ppg) d1 locs_ppg(k) - locs; if max(d1 win_lo d1 win_hi) matchCount matchCount 1; end end matchRate matchCount / length(locs_ppg);匹配率低于 0.9就直接把这个文件标记为低质量不参与后续统计而不是硬着头皮用被污染的特征。我现在的习惯是每个批次跑完先出一张三色标注图——R 波用红点、PPG 峰用绿点、收缩压峰用蓝点叠加在波形上肉眼扫一遍。图虽然土但几秒钟就能发现代码参数调出来的系统性偏移。项目做多了之后最大的教训就是特征提取这类流程指标再漂亮都不如一张标注图直观它能把 T 波误检、切迹混入、时间轴错位这类问题一眼照出来。希望这篇笔记能帮你少踩几个我当年踩过的坑让三路信号的特征工程一次跑顺。本文还有配套的精品资源点击获取