做可穿戴设备心率算法这么久我最大的感受是真正难的不是公式而是手里没有一套完全可控的数据。真实PPG生物信号里混着运动伪迹、呼吸分量、传感器压力变化带来的基线漂移算法算错了你都分不清是模型问题还是数据本身太脏。所以我每次搭心率提取链路都会先模拟一段PPG信号用短时傅里叶变换STFT做时频分析和滤波先把心率干净稳定地拉出来。这套流程跑通之后再拿真实数据去验证你才有底气判断问题是噪声环境还是算法缺陷。这篇就把信号怎么仿、STFT参数怎么调、掩膜怎么设计、心率怎么算一条线全摊开讲适合刚开始接触PPG心率算法或者正被真实数据折磨的同学参考。1. 为什么先搭一套仿真PPG再谈心率提取1.1 真实PPG数据看着简单用起来全是坑直接上真实PPG数据做算法第一反应是“波形挺有规律的”。但你真正处理起来会发现血压计的袖带晃一下、手表的表带松一格、手腕皮肤出汗波形就能变一个样。光电容积脉搏波传感器对光的吸收极其敏感肤色、毛发、汗液、温度都会叠加进去而这些东西根本不在你的算法模型里。更麻烦的是没有真值。你拿到一段60秒的真实PPG想提心率你知道答案是多少吗除非当时同步戴着一套心电图参考设备不然你算出一个72bpm根本无法确认它是准还是不准。我见过不少同学卡在这里算法在仿真数据上表现很好一换真实数据就乱跳然后开始疯狂调阈值、调滤波器的截止频率最后一晚上过去也没搞清楚到底是数据问题还是代码问题。所以我的习惯是仿真先行。先用已知心率生成PPG信号验证STFT滤波链路每一个环节的输出都符合预期再上真实数据。这样任何一步出了问题你都有一个没有歧义的基准能对照。1.2 PPG信号怎么来为什么能反映心率PPG全称光电容积脉搏波描记原理说穿了并不复杂心脏每次泵血时血管里的血容量会周期性变化而组织对光的吸收能力跟血容量直接相关。传感器一侧发光另一侧接收透射光或反射光光强就会跟着心跳周期波动把这束光强信号转成电信号就是PPG波形。典型PPG波形有一个快速上升的主波对应心脏收缩泵血随后是降中峡和重搏波对应动脉瓣关闭和血液回流。虽然波形形态上有这些细节但从频域看它的本质是一个以心率为基频的周期性信号。正常静息心率60次/分钟基频就是1Hz左右波形里的重搏波、谐波成分大多落在基频的整数倍上。所以心率提取的核心任务说到底是找出这个信号里最强的周期性频率成分然后乘60换算成bpm。这也解释了为什么模拟PPG可以用谐波叠加因为只要你把基频、二次谐波、三次谐波按一定幅度和相位叠加起来频谱结构和真实PPG已经非常接近了完全够用来验证STFT滤波的每一个环节。1.3 三种仿真建模思路选哪种我自己试过三种PPG仿真方式各有各的适用场景。第一种是谐波叠加法用正弦函数在心率基频及其倍频上叠加。优点是频率成分完全可控、代码量极少、适合验证频域算法缺点是时域波形跟真实PPG差距大主波和重搏波的形态不够精细。第二种是高斯脉冲序列法在每个心动周期放一个高斯形状的脉冲脉冲宽度模拟心脏泵血的时相。这种波形看起来比正弦更像PPG而且可以通过调节脉冲宽度来控制谐波的衰减速度适合做时域峰值检测类的验证。第三种是物理光学仿真通过建立皮肤、血管、光的相互作用模型生成信号。精度最高但计算量大、参数多做心率算法开发一般都杀鸡用牛刀。做STFT滤波验证我用的是第一种再额外叠加上呼吸分量、低频基线漂移、高斯白噪声。很多同学仿真时只生成一个干净正弦波滤波结果自然很漂亮但真实场景里呼吸的起伏大约在0.2到0.4Hz基线漂移在0.1Hz以下还有传感器本身的白噪声这些背景成分不加上验证出来的算法到了真实数据上会立刻原形毕露。2. 短时傅里叶变换对PPG滤波的三个关键优势2.1 FFT看全局STFT看局部非平稳信号的心率提取必须时频分析单纯用FFT做PPG频谱分析最大的问题是它只有一个时间段的整体统计。比如采集了60秒信号这60秒里心率从72bpm缓慢降到了65bpmFFT算出来只是一个宽带后的峰值你可能得到67bpm但完全看不出心率在变化。心率本质上是非平稳信号。人在说话、呼吸、姿势调整时心率会连续波动运动时还会在几秒内从80bpm冲上150bpm。FFT那种“全局平均”的思路天生不适合跟踪这种变化。而STFT做的事情很简单把信号切成一段一段对每一小段加窗后做FFT然后把这些短时频谱按时间拼成一个二维矩阵。横轴是时间纵轴是频率颜色深浅是能量强度你会看到一条随时间变化的能量亮带那就是心率基频的轨迹。类比一下FFT相当于你把整场电影的所有帧亮度做一个平均只能知道电影总体是亮还是暗STFT相当于按帧播放你能看到每一帧里哪个人在哪、灯是什么颜色。心率提取要的是后者因为心率本身随时在变。2.2 窗长、窗型、重叠率这三个参数决定滤波成败STFT的参数选择直接决定滤波效果这是我踩了最多坑的地方。先说频率分辨率公式是Δf fs / Nfs是采样率N是每一帧FFT的点数。假设PPG信号采样率100Hz窗口取1024点那么频率分辨率就是100/1024≈0.098Hz。心率是频率乘60所以对应的心率分辨率就是0.098×60≈5.9bpm。也就是说这个配置下你理论上能分辨72bpm和78bpm的心率。想要更高分辨率就把N加大比如2048点分辨率能到0.049Hz对应2.9bpm但窗口时间从10.24秒拉长到20.48秒时间分辨率变差心率快速变化时会看不清。窗型同样关键。矩形窗对信号没有幅度加权但它旁瓣非常高频谱泄漏严重会让一个本来很干净的1.2Hz峰值周围全是杂散能量时频图上看起来就是一条粗粗的糊带。汉宁窗的旁瓣明显更低虽然主瓣稍宽一点但整体频谱干净很多。做PPG滤波我一般默认汉宁窗。重叠率解决的是帧与帧之间的连贯性。相邻两帧如果完全不重叠心率轨迹会出现不自然的阶梯感。75%重叠意味着每帧移动N/4个采样点帧际平滑性会好很多STFT矩阵的时间分辨率也更密。2.3 先STFT滤波再做峰值检测为什么顺序不能反有人问我能不能直接用带通滤波器滤完然后再找峰值当然可以但实际调试中你在时域做滤波时很难直观看到“我滤掉的是什么”。而STFT掩膜滤波相当于构造了一个时变带通滤波器你在时频图上画一个掩膜把心率频带以外的所有时频点置零再做逆变换重建时域信号。这个过程的好处是你每一步都能看到时频图上的效果。真实PPG里二次谐波有时候比基频还亮运动伪迹会在某个时间段突然涌进宽频段。如果你直接在时域做固定带通滤波一个频率带从头到尾生效遇到伪迹时段往往无能为力。而STFT掩膜可以做成时变的比如某几帧伪迹很重可以把掩膜在那些时间点压得更窄。它本质上是一种“看你需要保留什么”的透明滤波比闭着眼睛调IIR滤波器系数直观得多。3. 从仿真信号到心率数值的完整实操流程3.1 搭好仿真信号把噪声加足别太小这是完整的仿真信号生成代码我用的是Python的numpy和scipy。import numpy as np from scipy import signal import matplotlib.pyplot as plt fs 100 # 采样率100Hz duration 60 # 60秒信号 t np.arange(duration * fs) / fs # 心率基频设为1.2Hz也就是72bpm f_hr 1.2 # 谐波叠加主体基波 二次谐波 三次谐波 ppg ( 1.0 * np.sin(2 * np.pi * f_hr * t) 0.4 * np.sin(2 * np.pi * 2 * f_hr * t 0.5) 0.2 * np.sin(2 * np.pi * 3 * f_hr * t 1.2) ) # 呼吸分量通常0.2~0.4Hz f_resp 0.25 resp 0.15 * np.sin(2 * np.pi * f_resp * t) # 基线漂移模拟传感器低频干扰 baseline 0.3 * np.sin(2 * np.pi * 0.05 * t) # 高斯白噪声 rng np.random.default_rng(42) noise 0.1 * rng.standard_normal(t.size) ppg_noisy ppg resp baseline noise为什么噪声幅度选0.1而不是0.01因为真实PPG传感器信号的信噪比并不高。0.1的噪声幅度会让波形看起来已经开始毛糙但还不至于淹没心率基频这是最贴近实际调试场景的强度。基线漂移幅度0.3比心率波形还大低频段在频谱上会非常醒目这样刚好检验STFT滤波能不能把它挡住。3.2 STFT时频分析先看谱再动手对仿真信号做STFT我用scipy封装好的函数。f, t_stft, Zxx signal.stft( ppg_noisy, fsfs, windowhann, nperseg1024, noverlap896 ) power np.abs(Zxx) # 只画0~5Hz呼吸和心率能量集中在这个范围 plt.figure(figsize(12, 4)) plt.pcolormesh(t_stft, f, power, shadinggouraud, cmapjet) plt.ylim(0, 5) plt.colorbar(labelMagnitude) plt.xlabel(Time [s]) plt.ylabel(Frequency [Hz]) plt.title(STFT of Simulated PPG) plt.show()我在这段配置下窗口点数1024重叠896帧移是128个采样点相当于时间分辨率约1.28秒。频率分辨率约0.098Hz既能看到心率基频的细微变化又不会丢太多时间细节。时频图上你应该能清楚看到0附近有一条很亮的低频横带那是基线漂移0.25Hz附近有另一条间或可见的细带那是呼吸分量1.2Hz附近是心率基频亮带而且它在整个60秒里保持稳定2.4Hz和3.6Hz处有两条暗一些的带那是二次和三次谐波。整个过程清晰直观。3.3 掩膜滤波与ISTFT重建怎么只留心率成分现在做核心的STFT域滤波。目标是把0.9到1.9Hz以外的所有时频成分压掉这个频带对应心率54到114bpm。fmin, fmax 0.9, 1.9 half_width 0.1 mask np.ones_like(Zxx, dtypefloat) mask[f fmin - half_width, :] 0 mask[f fmax half_width, :] 0 # 边缘斜坡避免硬截止引起振铃 in_low (f fmin - half_width) (f fmin) in_high (f fmax) (f fmax half_width) mask[in_low, :] (f[in_low, np.newaxis] - (fmin - half_width)) / (2 * half_width) mask[in_high, :] (1 - (f[in_high, np.newaxis] - fmax) / (2 * half_width)) _, ppg_filtered signal.istft(Zxx * mask, fsfs)掩膜边缘为什么要做斜坡而不是直接置零如果你在频域用一个矩形窗硬切相当于给频谱乘了一个矩形函数对应到时域会产生振铃重建出来的信号在波形首尾会出现可观的Gibbs现象。用0.1Hz宽的线性过渡带重建信号会平滑得多。逆变换之后ppg_filtered里的波形会变得非常干净低频漂移没了呼吸分量没了高频噪声也没了剩下的几乎是1.2Hz附近一个标准的正弦信号叠加了少量二次谐波成分。这条信号的包络能直观反映心率搏动节奏。3.4 两种心率计算法频谱峰值法和时域峰间隔法滤波之后有两种方式把心率算出来我都推荐试一下因为它们在真实场景下有不同的用处。方法一在STFT谱上找峰值。对每一帧的时间切片在0.9到1.9Hz频带里找幅度最大的频率点这个频率就是该帧的瞬时心率频率。idx_band (f fmin) (f fmax) peak_bin np.argmax(power[idx_band, :], axis0) hr_track f[idx_band][peak_bin] * 60 # 用中值滤波压掉毛刺 hr_median signal.medfilt(hr_track, kernel_size7) print(fEstimated HR: {np.median(hr_median):.1f} bpm)方法二对ISTFT重建后的时域信号做峰值检测。peaks, _ signal.find_peaks( ppg_filtered, distanceint(fs * 60 / 200), # 200bpm对应最小间隔 height0.3 ) hr_beat 60 / np.diff(t[peaks]) print(fBeat-to-beat HR: {np.median(hr_beat):.1f} bpm)两种方式各有特点。谱峰值法在心率变化连续时很稳适合看心率随时间的宏观趋势。时域峰间隔法更精细能得到每一次心跳的间隔用来算心率变异性HRV最合适。我的习惯是两个都跑如果结果相差超过1bpm就去查滤波链路有没有问题。4. 常见问题与排查技巧实录4.1 时频图上心率带“糊”成一片新上手最容易踩的坑是用矩形窗做STFT。我一开始就是图省事window默认起了个矩形窗结果时频图上1.2Hz附近的主峰两侧全是旁瓣泄漏看起来不是一条细线而是一条宽宽的亮带。你从谱上找峰值每帧找出来的频率都差好几Hz心率跟踪结果一跳一跳的。解决方法是把窗换成汉宁窗绝大多数情况下心率带变细立竿见影。如果换完还是糊检查nperseg是不是太小。比如nperseg256频率分辨率是0.39Hz对应心率精度约23bpm那你也就能区分45bpm和68bpm左右根本没法精确提取。把nperseg拉到1024以上情况立刻好很多。4.2 二次谐波被当心率模拟PPG里二次谐波幅度是基波的一半如果掩膜频带设到2.4Hz及以上2.4Hz这条谱线很容易被选成每帧峰值心率会算成144bpm正好是真实心率的2倍。排查的时候先看时频图如果亮带在1.2Hz和2.4Hz各有一条而你算出的心率稳定在144bpm那基本就是谐波误选。解决办法是收窄掩膜频带让上边截止频率低于二次谐波。比如把fmax设为1.9Hz二次谐波2.4Hz被全部滤掉。这也是为什么我前面频带取0.9到1.9而不是0到5Hz的原因。4.3 低频呼吸分量和基线漂移影响呼吸分量0.25Hz基线漂移0.05Hz这些低频成分幅度不小直接做时域峰值检测时波形上的鼓包和下沉可能让检测器把噪声峰当成心跳峰或漏检掉真实的心跳峰。STFT滤波对解决这个问题非常有效只要把掩膜的fmin抬到0.9Hz以上所有低频成分都会消失干净。但要注意一个边界场景如果测试者的心率极低比如50bpm对应0.83Hz而你把下限设成0.9Hz就会把真实心率也滤掉。这种时候把fmin往下调到0.7Hz不要机械沿用通用参数。4.4 参数速查表下面是我在实际项目中不同场景下常用的STFT配置直接抄作业也问题不大。应用场景采样率npersegnoverlap窗函数频率带离线心率复核100Hz20481792Hann0.7~2.5Hz常规心率提取100Hz1024896Hann0.9~1.9Hz运动伪迹较多100Hz512384Hamming1.0~3.0Hz实时低延迟估算100Hz256192Hann0.8~2.0Hz运动场景把通带上限拉到3.0Hz对应180bpm是因为运动心率会冲到160甚至180带外衰减太早会限制测算范围。实时低延迟场景窗口短分辨率低心率数值可以再配合抛物线插值提高精度。4.5 调试顺序与个人心得最后分享一个我自己调试这套流程的顺序希望能让你少走弯路。第一步永远是把原始信号的时频图画出来肉眼看清楚哪些是心率带、哪些是噪声再谈滤波。第二步先把掩膜频带设宽一点确认逆变换之后波形还在再把频带慢慢收窄去找最优边界。第三步跑心率计算之前先确认时频图亮带和滤波后波形峰值节奏对得上这一步能过滤掉80%的调试烦恼。还有一个细节是ISTFT重建信号的首尾几十个采样点通常不可靠因为加窗重建在边界处数据不足。做峰值检测时我习惯把前0.5秒和后0.5秒的数据直接裁掉避免边缘效应对心率统计造成干扰。还有一点想多说一句仿真信号虽然是正弦叠加出来的和真实PPG形态不完全一致但它是验证STFT链路是否正确的有效基准。只要仿真这一步不动摇换到真实数据上遇到问题时你至少能确定问题出在“传感器数据比仿真信号复杂得多”这一层而不是算法框架本身。