床垫这种产品做硬件的人一开始都觉得没什么门槛一块布、几片传感器、一个采集盒能出个在床/离床就交差了。真做到要输出心率和呼吸率这一步坑才开始一个接一个冒出来。我手头这个项目做了大概八个月从最早的PVDF压电薄膜方案到后来兼容压阻织物传感器中间换过三版算法。最后稳定下来的核心是用EMD算法把一路混在一起的体动信号拆成呼吸和心跳两条独立曲线再分别做频谱估计。这套心跳呼吸分离的流程我完整跑通了源代码也整理出来了下面连原理带代码一起讲清楚做PVDF传感器或者其他压阻人体体征传感器床垫的同行可以直接拿去改。先说清楚这套东西的适用边界。输入是单通道的体动信号采样率200Hz左右16位ADC输出是逐秒刷新的呼吸率次/分和心率bpm。传感器可以是PVDF压电薄膜也可以是压阻式导电织物、压阻橡胶、柔性应变片只要它能把胸腔的机械形变转成电信号就行。不挑传感器型号挑的是信号质量和采样链路的设计。适合做智能床垫、智能坐垫、婴儿监护垫的工程师看也适合做生理信号处理的学生参考代码全部是Python依赖只有numpy和scipy。1. 床垫体征监测的整体设计思路做这个项目之前我调研过好几条技术路线光电、雷达、压电、压阻都试过。雷达方案贵、功耗高光电方案在床垫上根本没法贴合最后落到压电和压阻这两条路上。它们的共同点是不主动发射任何东西纯被动感知机械形变结构简单成本可控塞进床垫里不影响睡感。麻烦的地方在于传感器出来的是一路信号呼吸和心跳全叠在里面怎么把这俩拆开就是整个项目的技术核心。1.1 PVDF和压阻式传感器在床垫场景下的优劣对比PVDF压电薄膜的工作原理是薄膜受力产生形变时内部偶极子取向变化两端出现电荷积累。它是典型的动态传感器输出的是电荷量形变越快输出越大静止不动输出就是零。这个特性用在床垫上有两个好处一是灵敏度高胸腔那点微米级的起伏它都能感应到二是天然隔直缓慢的体温漂移、床垫受压后的静态形变这类干扰不会跑进信号里。坏处也来自于同一个特性。它对低频响应差而呼吸正好是低频0.1到0.5Hz。如果电荷放大器的反馈电阻不够大呼吸信号会被削得很厉害后面EMD再厉害也救不回来。这一点后面会详细算。压阻式传感器的逻辑完全相反。它的电阻随压力变化需要外部激励恒压或恒流输出是电压变化能做到直流响应测静态压力没问题。缺点是灵敏度相对低而且本身有热漂移长时间工作基线会跑。好在你测的是呼吸和心跳这种交变成分把基线去掉就行。我把两者的对比整理成一张表方便选型时对照对比维度PVDF压电薄膜压阻式导电织物/压阻橡胶输出类型电荷需电荷放大器电阻变化需激励源低频响应差取决于Rf·Cf好可到DC灵敏度高中等静态压力无法测量可测量功耗极低无源需持续激励布线难度需屏蔽走线敏感相对不敏感成本中等偏高低适合场景呼吸心跳的精细波形在床/离床呼吸粗心率实际项目里我推荐的组合是PVDF做主传感压阻织物做辅助前者负责波形质量后者负责在床判断和姿态区分。如果预算只够一种做心率就选PVDF只做呼吸和在床检测选压阻就够了。1.2 呼吸和心跳为什么会粘在一起很多人第一次看床垫的原始波形会困惑明明波形挺干净的为什么频谱一出来是一堆乱七八糟的峰。原因在于两个信号本身的性质差异极大。呼吸的频率范围大概在0.1到0.5Hz对应每分钟6到30次。它的幅度大胸廓起伏带来的形变可能是心跳的10到50倍。心跳的频率在0.8到2.2Hz对应每分钟48到132次覆盖了绝大部分成年人的静息心率。幅度小但频率高。问题出在两个地方。一是频谱上虽然不重叠但真实信号的呼吸波形不是纯正弦它有二次、三次谐波0.3Hz的呼吸二次谐波就跑到0.6Hz三次谐波接近0.9Hz直接钻进心跳的频带里。二是呼吸过程中胸腔位置在变传感器和心脏之间的耦合强度也跟着变心跳信号的幅度被呼吸调制了这在时域上看就是心跳的包络跟着呼吸一起起伏。这两条决定了单纯用固定参数的带通滤波器拆不干净。滤波器只能按频率切切不掉谐波也处理不了幅度调制带来的畸变。这就是我最后转向EMD的直接原因。1.3 选EMD而不是固定带通滤波的三个理由EMD全称经验模态分解Empirical Mode Decomposition它的核心思想是任何复杂信号都可以分解成有限个本征模态函数IMF每个IMF是一个窄带分量有自己独立的瞬时频率和瞬时幅度。这个过程完全靠数据自己驱动不需要你预先指定任何基函数。选它有三个理由。第一它是自适应的。每个人的呼吸频率不一样同一个人入睡前后也不一样。固定带通滤波器的通带一旦设死遇到呼吸很慢的人比如8次/分0.13Hz滤波器通带如果从0.15Hz起呼吸主频就被削了。EMD不需要预设它自己会把这个频率上的成分单独拆成一个IMF。第二它能处理非平稳信号。人在睡眠中呼吸和心率的频率是缓慢漂移的傅里叶变换假设信号在窗口内平稳时间窗口一长就失真。EMD是基于局部极值点做的天然适应慢变。第三它对谐波的处理更自然。呼吸的三次谐波如果能量够大会被拆到单独的IMF里而不会污染心跳所在的IMF。当然这块也不是完美的后面会讲模态混叠这个坑。代价是计算量比滤波大以及端点效应、停止准则这些需要调。但在30到60秒的分析窗上跑PC端完全无压力嵌入式端需要做定点化和降采样优化。2. 硬件链路与采样参数怎么定算法再好前端信号烂了也是白搭。这一节把我在这块踩过的坑集中讲一下尤其是电荷放大器的低频截止频率计算我见过太多项目在这里翻车呼吸波形被削成一条直线还以为是算法问题。2.1 电荷放大器与低频截止的频率计算PVDF输出的是电荷不能直接接ADC中间必须有一个电荷放大器。基本结构是一个运放反馈回路里挂一个电容Cf和一个电阻Rf。输出电压为V_out -Q / Cf其中Q是传感器产生的电荷量。Cf决定了增益Cf越小增益越大。比如传感器在呼吸时产生1pC的电荷Cf取10nF输出就是0.1mV。这个量级很弱需要ADC前级再做一级放大或者把Cf降到1nF。真正关键的是Rf。它决定了低频截止频率f_c 1 / (2π · Rf · Cf)这个公式很直白Rf和Cf构成一个高通网络低于f_c的成分会被衰减。呼吸是低频信号所以f_c必须设得远低于呼吸的最低频率。我来算一下。假设Cf固定取10nF为了让0.1Hz的呼吸成分衰减不超过3dB即幅值保持在0.707以上f_c需要满足f_c ≤ 0.1 / sqrt(1/0.707² - 1) ≈ 0.1 Hz也就是说f_c至少要压到0.1Hz以下最好压到0.03Hz以下留足余量。取f_c 0.016Hz则Rf 1 / (2π × 0.016 × 10e-9) ≈ 1e9 Ω 1 GΩ1GΩ的电阻不是标准件通常用多颗高阻电阻串联或者用运放的反馈T型网络等效实现。我实际用的方案是两颗500MΩ的高阻电阻串联并联一个小电容做补偿。如果用100MΩ会怎样f_c 1/(2π × 1e8 × 10e-9) 0.159Hz。在0.1Hz处增益衰减为 0.1/sqrt(0.1² 0.159²) 0.53也就是-5.5dB。呼吸幅度直接掉一半而且不同呼吸频率衰减程度还不一样导致波形严重失真。这个坑我是真踩过当时排查了两天才发现是电阻选小了。注意电荷放大器的输入端是高阻节点PCB上必须做防护环guard ring走线尽量短否则漏电流会直接淹没信号。另外传感器的屏蔽层要接放大器地不能两端都接地形成地环路。2.2 采样率和位数的选择推演采样率的选择遵循奈奎斯特准则但工程上要留足余量。心跳最高频率取2.2Hz132bpm按照10倍过采样原则采样率至少22Hz。呼吸0.5Hz更没压力。但这里有一个容易被忽略的点工频干扰。如果采样率取50Hz奈奎斯特频率是25Hz。50Hz工频高于奈奎斯特频率会混叠。混叠到哪里采样频率是50Hz50Hz正好是采样率的整数倍混叠到0Hz表现为基线漂移。更糟的是如果采样率是49.9Hz50Hz会混叠到0.1Hz正好落在呼吸频带里这种情况你无论怎么滤波都救不回来因为它在数字域里和呼吸长得一模一样。结论采样率不能取50Hz。我推荐200Hz。这样奈奎斯特频率100Hz50Hz工频落在中间用数字陷波器可以精确干掉100Hz的二次谐波也能处理。那能不能取250Hz或者更高可以但没必要。采样率越高单窗口的数据点越多EMD的计算量线性增长。200Hz降采样到50Hz后30秒窗口是1500点EMD分解一层大概几毫秒整帧几十毫秒实时性完全够。ADC位数选16位。原因在于动态范围。呼吸信号幅度可能是心跳的20倍以上如果只有12位4096级留给心跳的分辨率就只剩200级波形会明显量化。16位是65536级心跳能分到3000级左右足够做谱分析了。2.3 预处理的四道工序拿到200Hz的原始数据后不能直接丢给EMD中间要做四件事顺序不能乱。第一道是工频陷波。用二阶IIR陷波器中心频率50HzQ值取30。如果采样率足够高比如400Hz以上再加一个100Hz的陷波。陷波器要用零相位滤波filtfilt否则会引入相位失真导致后面的峰间隔计算出现系统性偏差。第二道是滑动中值去基线。翻身、调整睡姿会让床垫压力分布发生阶跃式变化表现为一个大的台阶或缓慢漂移。用宽度2秒的滑动中值滤波估计基线再从原信号里减掉。这个方法比高通滤波好因为中值滤波对阶跃的响应是跟着走不会产生振铃。第三道是带通限幅。通带设0.05到10Hz就够了下限保呼吸上限保心跳的高次谐波并抑制高频噪声。10Hz以上的成分在这个应用里没有任何价值。第四道是重采样到50Hz。用多相滤波重采样resample_poly不要用简单抽取否则会引入混叠。50Hz下心跳2.2Hz有22倍采样足够。实操心得这四道工序里中值滤波的窗口宽度是唯一需要跟着床垫软硬度调的参数。硬床垫的翻身瞬变更快窗口要窄一些1.5秒软床垫形变释放慢窗口可以放到3秒。3. EMD分解的代码实现与关键细节这一节是全文的核心。我会把EMD的实现逻辑拆开讲包括筛分循环、端点延拓、停止准则然后是IMF的自动挑选。最后给出完整代码。代码我自己跑过很多次不是抄来的每个函数都是按实际需求写的。3.1 筛分循环与停止准则EMD的分解过程叫筛分sifting步骤是这样的找出信号的所有局部极大值点和局部极小值点用三次样条分别拟合上包络和下包络计算上下包络的均值曲线m(t)用原信号减去均值曲线得到h(t) x(t) - m(t)检查h(t)是否满足IMF的两个条件极值点数和过零点数相差不超过1上下包络均值处处为零如果不满足把h(t)当作新信号重复1到5如果满足把h(t)作为第一个IMF输出从原信号里减掉它对残差重复整个过程停止准则这里有个常见的做法分歧。Huang在1998年提出的是标准差准则SD准则即连续两次筛分结果的标准差小于阈值就停SD Σ[(h_{k-1}(t) - h_k(t))² / h_{k-1}²(t)]这个阈值一般取0.2到0.3。太大则IMF不够纯太小则迭代次数暴涨而且可能把信号筛成纯粹的调幅波丢失物理意义。另一个是极值点准则如果连续两次筛分后极值点数量和位置基本不变就停止。这个准则在实时的工程实现里更实用因为计算量小而且能避免过度筛分。我实际用的是混合策略SD阈值0.2作为主准则同时设最大迭代次数50次兜底。为什么要有兜底因为在某些噪声段信号几乎没有极值点筛分会陷入死循环。我见过一次没有兜底的实现在一段静默数据上卡了十几秒。3.2 端点效应、包络插值的坑EMD最臭名昭著的问题就是端点效应。原因很简单三次样条插值需要边界条件而信号的第一个极值点和最后一个极值点外面没有数据了。样条在这两段会剧烈发散导致包络在两端严重失真而且这个失真会随着筛分一层层向内传播甚至污染整个分解结果。解决办法有两类。一类是延拓法在信号两端人为补一些数据点让样条有足够的支撑。镜像延拓最简单把第一个极值点关于起点做镜像放到信号左边把最后一个极值点关于终点镜像放到右边。这个方法实现简单效果能接受我用的就是它。另一类是改进样条比如用B样条或者有理样条替代三次样条。效果更好但实现复杂对实时系统不划算。还有一个坑是极值点数量不足。当信号很短或者很平滑时极大值点可能只有一个甚至没有样条无法拟合。这时候必须直接返回把这个分量作为残差输出而不是硬拟合然后崩掉。代码里的len(max_idx) 2判断就是这个作用。另外三次样条插值对极值点的密集程度很敏感。如果信号里高频噪声多极值点会非常密集样条过拟合包络变成锯齿状分解出来的IMF全是噪声。所以预处理那一步的带通限幅非常必要把10Hz以上的噪声干掉极值点数量就正常了。3.3 用频率和能量准则自动挑选IMF分解出一堆IMF之后哪一个是呼吸哪一个是心跳最直接的办法是算每个IMF的平均频率。用零点穿越法估算统计一段信号里的过零点数量除以2再除以时长就是平均频率。这个方法对窄带信号很准实现也简单f_mean (过零点数 / 2) / 窗口时长呼吸IMF的频率应该在0.1到0.6Hz之间心跳IMF在0.8到2.5Hz之间。注意上限给到2.5Hz而不是2.2Hz是为了留余量防止心率高的时候比如运动后上床被漏掉。但光看频率还不够。有时候会有两个IMF都落在呼吸频段里这时候要看能量占比。能量占比定义为该IMF的能量与原始信号能量之比E_ratio Σ(imf²) / Σ(x²)同一个频段里能量占比最大的那个才是主成分其他的是谐波或噪声。这里有一个重要的经验EMD分解出来的第一个IMF最高频通常是噪声最后一个IMF最低频通常是趋势项中间的才是有效成分。我在挑选时会直接从IMF2开始扫跳过IMF1。选中呼吸IMF和心跳IMF后还有一步后处理对心跳IMF再做一次0.7到3.0Hz的带通。为什么因为EMD的模态混叠问题没有完全解决心跳IMF里常常残留呼吸的谐波成分。这一步带通能把这个残留压下去把心率估计的标准差从5bpm降到2bpm左右。3.4 完整可运行的Python源代码下面是完整代码包含预处理、EMD实现、IMF挑选、呼吸率心率估计的全部流程。直接复制到一个.py文件里就能跑。为了演示方便代码末尾生成了一个合成的床垫信号呼吸0.25Hz 心跳1.2Hz 噪声实际使用时把synthesize_mattress_signal换成你的采集数据即可。# -*- coding: utf-8 -*- 智能床垫体征算法EMD分解 心跳呼吸分离 传感器PVDF压电薄膜 / 压阻式导电织物 / 柔性应变片 输出呼吸率(次/分)、心率(bpm)、置信度标记 依赖numpy, scipy import numpy as np from scipy.signal import (butter, filtfilt, iirnotch, find_peaks, resample_poly, hilbert) from scipy.interpolate import CubicSpline from scipy.ndimage import median_filter # # 一、预处理 # def notch_filter(x, fs, f050.0, q30.0): 工频陷波零相位 if f0 fs / 2.0 * 0.95: return x b, a iirnotch(f0, q, fs) return filtfilt(b, a, x) def detrend_median(x, fs, win_sec2.0): 滑动中值去基线专治翻身台阶漂移 w int(win_sec * fs) if w % 2 0: w 1 if w 3: return x base median_filter(x, sizew, modenearest) return x - base def bandpass(x, fs, lo, hi, order4): 零相位带通 nyq fs / 2.0 lo_n max(lo / nyq, 1e-4) hi_n min(hi / nyq, 0.999) if lo_n hi_n: return x b, a butter(order, [lo_n, hi_n], btypeband) return filtfilt(b, a, x) def preprocess(x, fs_in, fs_out50.0): 完整预处理链陷波 - 去基线 - 带通 - 重采样 x np.asarray(x, dtypefloat) x x - np.mean(x) # 1) 工频陷波 x notch_filter(x, fs_in, 50.0, q30.0) if fs_in 250: x notch_filter(x, fs_in, 100.0, q30.0) # 2) 去基线漂移 x detrend_median(x, fs_in, win_sec2.0) # 3) 带通限幅 0.05 ~ 10 Hz x bandpass(x, fs_in, 0.05, 10.0, order4) # 4) 重采样到 50 Hz if abs(fs_in - fs_out) 1e-6: up, down _ratio(fs_out, fs_in) x resample_poly(x, up, down) return x, fs_out def _ratio(fs_out, fs_in, max_den200): 把 fs_out/fs_in 化简成整数比 from fractions import Fraction fr Fraction(fs_out / fs_in).limit_denominator(max_den) return fr.numerator, fr.denominator # # 二、EMD 核心实现 # def _extrema_idx(x): 定位局部极大值和极小值限制最小间隔避免噪声导致极值点爆炸 min_gap max(1, len(x) // 200) max_idx, _ find_peaks(x, distancemin_gap) min_idx, _ find_peaks(-x, distancemin_gap) return max_idx, min_idx def _spline_envelope(x, idx, n): 镜像延拓 三次样条拟合包络 if len(idx) 2: return None # 左端镜像把第一个极值点关于起点镜像 left_pos -idx[0] left_val x[idx[0]] # 右端镜像把最后一个极值点关于终点镜像 right_pos 2 * (n - 1) - idx[-1] right_val x[idx[-1]] pos np.concatenate(([left_pos], idx, [right_pos])) val np.concatenate(([left_val], x[idx], [right_val])) # 去掉重复位置避免样条报错 pos_u, uniq np.unique(pos, return_indexTrue) val_u val[uniq] if len(pos_u) 4: return None cs CubicSpline(pos_u, val_u, bc_typenatural) return cs(np.arange(n)) def _sift(x, sd_thresh0.2, max_iter50): 单次筛分返回一个IMF h x.astype(float).copy() for _ in range(max_iter): max_idx, min_idx _extrema_idx(h) if len(max_idx) 2 or len(min_idx) 2: return h, True # 极值点不足无法继续 upper _spline_envelope(h, max_idx, len(h)) lower _spline_envelope(h, min_idx, len(h)) if upper is None or lower is None: return h, True m 0.5 * (upper lower) h_new h - m denom np.sum(h ** 2) 1e-12 sd np.sum((h - h_new) ** 2) / denom h h_new if sd sd_thresh: break return h, False def emd(x, max_imf8, sd_thresh0.2, max_iter50): 经验模态分解返回IMF数组和残差 imfs [] r np.asarray(x, dtypefloat).copy() for _ in range(max_imf): max_idx, min_idx _extrema_idx(r) if len(max_idx) 2 or len(min_idx) 2: break imf, _ _sift(r, sd_thresh, max_iter) if np.sum(imf ** 2) 1e-12: break imfs.append(imf) r r - imf if len(imfs) 0: return np.zeros((0, len(x))), r return np.array(imfs), r # # 三、IMF 挑选 # def imf_mean_freq(imf, fs): 零点穿越法估平均频率 sign np.sign(imf) sign[sign 0] 1 zc np.sum(np.diff(sign) ! 0) return zc / 2.0 / (len(imf) / fs) def imf_energy_ratio(imf, x): return np.sum(imf ** 2) / (np.sum(x ** 2) 1e-12) def select_imf(imfs, x, fs, band, min_energy0.005): 在指定频带内挑选能量占比最大的IMF lo, hi band cand [] for i, imf in enumerate(imfs): f imf_mean_freq(imf, fs) e imf_energy_ratio(imf, x) if lo f hi and e min_energy: cand.append((i, f, e)) if not cand: return None, None, [] cand.sort(keylambda t: -t[2]) idx cand[0][0] return imfs[idx], idx, cand # # 四、频率估计 # def _parabolic_refine(freqs, amps, k): 抛物线插值细化谱峰位置 if 0 k len(amps) - 1: a, b, c amps[k - 1], amps[k], amps[k 1] denom a - 2 * b c if abs(denom) 1e-12: delta 0.5 * (a - c) / denom if abs(delta) 1.0: return freqs[k] delta * (freqs[1] - freqs[0]) return freqs[k] def spectrum_peak(x, fs, fmin, fmax, zero_pad8): 加汉宁窗 补零 抛物线插值返回频带内主峰频率 n len(x) if n int(fs / fmin * 2): return np.nan, 0.0 w np.hanning(n) xw (x - np.mean(x)) * w nfft 1 int(np.ceil(np.log2(n * zero_pad))) X np.abs(np.fft.rfft(xw, nfft)) freqs np.fft.rfftfreq(nfft, 1.0 / fs) mask (freqs fmin) (freqs fmax) if not np.any(mask): return np.nan, 0.0 sub_f freqs[mask] sub_a X[mask] k int(np.argmax(sub_a)) peak_f _parabolic_refine(sub_f, sub_a, k) # 谱峰突出度主峰 / 频带中位数 prom sub_a[k] / (np.median(sub_a) 1e-12) return peak_f, prom def peak_interval_rate(x, fs, fmin, fmax): 峰值间隔法返回中位频率(Hz)和间隔数 if fmax fmin: return np.nan, 0 min_dist int(fs / fmax * 0.7) pk, _ find_peaks(x, distancemax(1, min_dist)) if len(pk) 3: return np.nan, 0 ibi np.diff(pk) / fs lo, hi 1.0 / fmax, 1.0 / fmin ibi ibi[(ibi lo) (ibi hi)] if len(ibi) 2: return np.nan, 0 return 1.0 / np.median(ibi), len(ibi) # # 五、单帧分析主流程 # def analyze_frame(x_raw, fs_in200.0): 输入一段原始信号输出呼吸率、心率及置信信息 x, fs preprocess(x_raw, fs_in, fs_out50.0) n len(x) duration n / fs result { rr_bpm: np.nan, # 呼吸率 次/分 hr_bpm: np.nan, # 心率 bpm rr_hz: np.nan, hr_hz: np.nan, rr_spec_conf: 0.0, hr_spec_conf: 0.0, hr_energy: 0.0, rr_imf_idx: -1, hr_imf_idx: -1, quality: invalid } if duration 20: return result imfs, res emd(x, max_imf8, sd_thresh0.2, max_iter50) if imfs.shape[0] 2: return result # --- 呼吸IMF0.10 ~ 0.60 Hz resp_imf, resp_i, _ select_imf(imfs, x, fs, (0.10, 0.60), min_energy0.01) # --- 心跳IMF0.80 ~ 2.50 Hz card_imf, card_i, _ select_imf(imfs, x, fs, (0.80, 2.50), min_energy0.002) # ---------- 呼吸率 ---------- if resp_imf is not None: f1, c1 spectrum_peak(resp_imf, fs, 0.10, 0.60) f2, n2 peak_interval_rate(resp_imf, fs, 0.10, 0.60) if np.isfinite(f1) and np.isfinite(f2): if abs(f1 - f2) * 60 3.0: rr_hz 0.5 * (f1 f2) else: rr_hz f1 # 冲突时以谱峰为准 elif np.isfinite(f1): rr_hz f1 elif np.isfinite(f2): rr_hz f2 else: rr_hz np.nan if np.isfinite(rr_hz): result[rr_hz] rr_hz result[rr_bpm] rr_hz * 60.0 result[rr_spec_conf] c1 result[rr_imf_idx] int(resp_i) # ---------- 心率 ---------- if card_imf is not None: # 关键一步对心跳IMF再做一次带通压掉呼吸谐波残留 card_bp bandpass(card_imf, fs, 0.7, 3.0, order4) f1, c1 spectrum_peak(card_bp, fs, 0.80, 2.50) f2, n2 peak_interval_rate(card_bp, fs, 0.80, 2.50) if np.isfinite(f1) and np.isfinite(f2): if abs(f1 - f2) * 60 6.0: hr_hz 0.5 * (f1 f2) else: hr_hz f1 elif np.isfinite(f1): hr_hz f1 elif np.isfinite(f2): hr_hz f2 else: hr_hz np.nan if np.isfinite(hr_hz): result[hr_hz] hr_hz result[hr_bpm] hr_hz * 60.0 result[hr_spec_conf] c1 result[hr_imf_idx] int(card_i) result[hr_energy] imf_energy_ratio(card_imf, x) # ---------- 质量评估 ---------- ok_rr np.isfinite(result[rr_bpm]) and result[rr_spec_conf] 5 ok_hr (np.isfinite(result[hr_bpm]) and result[hr_spec_conf] 4 and result[hr_energy] 0.01) if ok_rr and ok_hr: result[quality] good elif ok_rr or ok_hr: result[quality] partial return result # # 六、演示合成一段床垫信号 # def synthesize_mattress_signal(fs200.0, dur60.0, rr0.25, hr1.20, seed7): rr: 呼吸频率Hz hr: 心率Hz rng np.random.default_rng(seed) t np.arange(0, dur, 1.0 / fs) # 呼吸含二次、三次谐波模拟真实非正弦呼吸 resp (1.00 * np.sin(2 * np.pi * rr * t) 0.22 * np.sin(2 * np.pi * 2 * rr * t 0.7) 0.08 * np.sin(2 * np.pi * 3 * rr * t 1.9)) # 心跳幅度被呼吸调制胸腔耦合强度随呼吸变化 mod 1.0 0.35 * np.sin(2 * np.pi * rr * t 0.4) card 0.045 * mod * np.sin(2 * np.pi * hr * t) # 工频 白噪声 hum 0.01 * np.sin(2 * np.pi * 50.0 * t) noise rng.normal(0, 0.004, len(t)) return resp card hum noise, fs if __name__ __main__: sig, fs synthesize_mattress_signal(fs200.0, dur60.0, rr0.25, hr1.20) out analyze_frame(sig, fs_infs) print(呼吸率: %.1f 次/分 (真值 15.0) % out[rr_bpm]) print(心率 : %.1f bpm (真值 72.0) % out[hr_bpm]) print(呼吸IMF编号: %d, 心跳IMF编号: %d % (out[rr_imf_idx], out[hr_imf_idx])) print(心跳IMF能量占比: %.4f % out[hr_energy]) print(质量标记: %s % out[quality])这段代码里我最想强调的是analyze_frame里对心跳IMF再带通那一步。很多人做完EMD就直接对IMF做FFT结果心率估计一直跳。原因就是模态混叠呼吸的三次谐波有时候会跑进心跳IMF里谱峰就被带偏了。加一道带通成本极低效果立竿见影。4. 呼吸率与心率的提取及实测验证算法跑通只是第一步能不能输出稳定的数值才是产品化的门槛。这一节讲我怎么做频率估计、怎么交叉校验、以及实测下来误差有多大。4.1 呼吸率的两种算法和交叉校验呼吸率的估计我用了两种方法并行。第一种是谱峰法。对呼吸IMF加汉宁窗补零到8倍长度做FFT在0.1到0.6Hz范围内找主峰然后用抛物线插值细化峰位。补零的作用是提高频率轴的采样密度让抛物线插值更有意义。不加补零的话30秒窗口的频率分辨率只有0.033Hz换算成呼吸率是2次/分误差太大。抛物线插值的效果值得说一下。假设真实呼吸是0.23Hz不加插值时谱峰只可能落在0.233Hz离散步长误差0.003Hz也就是0.18次/分。看起来还行但对于更短的分析窗比如20秒离散步长变成0.05Hz误差就上去了。加了抛物线插值能把峰位细化到离散步长的十分之一左右误差稳定在0.1次/分以内。第二种是峰值间隔法。直接在呼吸IMF上找波峰算相邻峰的间隔取中位数然后倒推频率。这个方法的好处是时间分辨率高能快速跟踪呼吸频率的变化。坏处是对噪声敏感如果IMF上有小的杂峰会被误判成呼吸峰。所以我在find_peaks里加了最小间隔限制最小间隔是最高呼吸频率对应周期的0.7倍。两种方法的结果做一致性检查如果差异小于3次/分取平均如果大于3次/分以谱峰法为准同时把这帧标记为低置信度。实测下来呼吸平稳时两者差异通常在1次/分以内翻身时可能出现较大分歧这时候的帧本来也不该被采信。4.2 心率提取从IMF到包络谱心率的估计比呼吸难原因有三信号弱、容易受运动干扰、频带和高频噪声相邻。主路径是对心跳IMF做带通后的谱峰估计和呼吸率的方法一样。但心率的频带更宽0.8到2.5Hz谐波干扰的概率更高。我在这里加了一个额外的可靠性判断谱峰突出度。定义为频带内主峰幅度除以频带内所有谱线幅度的中位数。一个干净的心跳信号这个比值通常在8以上如果低于4说明频带内能量分散没有明显的主频这帧数据不可信。另外还有一个交叉校验用的方法叫包络法。它的思路是当传感器和心脏之间隔了厚被子或者很软的床垫层时高频的搏动细节被平滑掉了心跳表现为一个窄带载波其幅度被心脏每次搏动调制。这时候对心跳IMF做希尔伯特变换取包络包络的起伏频率就是心率。提示包络法不是万能的。如果传感器直接耦合良好心跳在IMF上就是一个准正弦希尔伯特包络会接近常数包络法直接失效。我的做法是两种方法都算谁的谱峰突出度高就用谁。希尔伯特包络的实现很简单def envelope_rate(imf, fs, fmin0.8, fmax2.5): 希尔伯特包络谱法适用于心跳被平滑成载波的场景 env np.abs(hilbert(imf)) env env - np.mean(env) # 包络本身是低频信号重采样到25Hz足够 env_ds resample_poly(env, 1, 2) if fs 50 else env fs_ds fs / 2 if fs 50 else fs f, conf spectrum_peak(env_ds, fs_ds, fmin, fmax, zero_pad8) return f, conf注意包络法需要先降采样因为包络的有效带宽只有几Hz用50Hz采样做FFT是浪费而且高频噪声会污染包络谱。降到25Hz后2.5Hz的奈奎斯特余量还有5倍足够了。4.3 与参考设备的对比记录我找了12个人做对比测试其中8个是我同事4个是找的外部志愿者。测试方法受试者仰卧在装了PVDF传感器的样机上同时佩戴指夹式血氧仪测心率胸前绑一根应变带测呼吸。每人采集10分钟静息数据取躺下2分钟之后的8分钟数据做统计。数据分帧方式60秒窗口滑动步长1秒每帧输出一个估计值。把8分钟里所有good质量标记的帧取中位数作为该受试者的最终结果。统计下来的结果指标平均绝对误差95分位误差呼吸率1.3 次/分2.8 次/分心率2.1 bpm4.6 bpm这个精度对于非医疗级的睡眠监测产品是够用的。作为对照我最早用固定带通滤波那一版心率的平均绝对误差是5.4 bpm95分位误差超过11 bpm完全不能用。还有几个观察值得记录。第一体重轻的人BMI低于19心率误差明显大因为胸腔传导到床垫的机械能量小心跳IMF的能量占比往往低于1.5%接近我设的置信度门槛。第二呼吸率在浅睡期误差会略大因为呼吸变得不规则波形不是周期性的谱峰会变宽。第三翻身动作后大约5到8秒内结果不可信我在实际产品上用加速度计触发一个运动屏蔽窗口这段时间的输出保持上一帧的值不变。5. 常见问题排查与调参经验谈写了这么多原理和代码最后这部分是我最想分享的。因为上面那些东西你看书看论文都能找到但下面这些坑只有真做过床垫项目的人才知道。5.1 问题速查表我把这八个月里遇到过的典型问题整理成表遇到类似现象可以先在这里对一下。现象可能原因排查手段解决方式呼吸率恒为某个固定值不变化呼吸信号被削平算法在拟合噪声画呼吸IMF时域波形看是否有起伏检查电荷放大器Rf是否足够大f_c是否低于0.03Hz心率忽高忽低跳动超过15bpm心跳IMF能量弱谱峰被噪声主导打印每帧心跳IMF的能量占比提高置信度门槛能量低于2%直接丢弃该帧呼吸率是真实值的2倍或3倍呼吸谐波被选为呼吸IMF打印所有IMF的中心频率和能量在候选IMF里强制选最低频的那个或加谐波约束每帧结果差异大重启后不一致端点效应导致分解结果不确定对比不同起点切窗的分解结果增加镜像延拓帧间重叠50%多帧取中位数夜间某个时段数据全废翻身或离床看原始信号RMS是否有突增加运动屏蔽窗口RMS超阈值时冻结输出CPU占用率高实时性跟不上EMD在高采样率下计算量大统计单帧耗时先降采样到25Hz再做EMD心跳2.5Hz仍有10倍余量心率偏高约等于呼吸率的整数倍模态混叠呼吸谐波混进心跳IMF对比呼吸IMF和心跳IMF的时域相关性心跳IMF后接0.7-3.0Hz带通不同人之间精度差异巨大传感器位置和耦合强度不一致记录每个受试者的心跳IMF能量占比产品化时做一次个体校准记录基准能量水平5.2 几个我自己踩出来的调参心得第一个心得关于EMD的分解层数。max_imf参数我设的是8但实际有效的通常只有4到5层。设太大的问题是后面几层全是趋势项白算。设太小的问题是心跳IMF可能出不来。我的经验是先设12层跑一次把每层的中心频率和能量占比打出来看清楚心跳落在第几层然后再把max_imf收敛到那个层数加2。不同床垫不一样硬床垫心跳IMF往往在第3层软床垫因为高频被吸收会落到第2层。第二个心得关于分析窗口长度。论文里常见的是1秒到5秒的窗口那是给心电信号用的。床垫信号里心跳太弱短窗口的谱分辨率根本不够。我试过5秒、10秒、20秒、30秒、60秒最后定在60秒窗口加1秒步长。60秒的频率分辨率是0.017Hz配合抛物线插值能到0.002Hz对应心率误差0.12bpm完全够。代价是心率变化有大约30秒的滞后因为窗口的一半都在过去对于睡眠监测这种慢变场景可以接受。如果要跟踪运动后的心率恢复得用卡尔曼滤波做平滑和延迟补偿。第三个心得关于置信度的设计。一开始我没有置信度算法什么都输出结果用户看到心率从60突然跳到120又跳回来体验极差。后来加了三个门槛谱峰突出度、IMF能量占比、呼吸和心率的合理性检查心率必须大于呼吸率的2倍且心率在40到200之间。三个都过了才输出否则保持上一帧的值。这个改动让数值的平滑度提升了一个数量级用户侧看起来就是偶尔卡一下而不是乱跳。第四个心得关于模态混叠的处理。这是个老大难问题标准EMD解决不了。我试过EEMD加白噪声的集合平均效果好但计算量翻十倍实时系统扛不住。后来用了一个折中方案在预处理阶段把10Hz以上的噪声干掉让极值点分布更均匀模态混叠的概率就明显降低了。另外就是前面反复强调的对心跳IMF再做一次带通这个是最实用的补丁。实操心得如果要进一步提升心跳的分离质量可以试试带掩膜信号的EMDMasking EMD也就是在分解前往信号里叠加一个已知频率的高频正弦把这个正弦的频率设在心跳频带之外让心跳IMF的特征更明显分解完再减掉。这个方法比EEMD轻量很多我在一块低功耗MCU上跑过单帧耗时大概80ms。最后再提一个设备端的优化方向。EMD的主要开销在三次样条插值和多次筛分这两块都可以做定点化。样条的系数计算用查表加线性插值近似误差在1%以内但速度能快3倍。筛分次数限制在8次以内对结果的影响很小。我把这些优化做完之后200Hz采样、60秒窗口的整帧处理在Cortex-M4上大概需要200ms左右如果只要心跳和呼吸率这两个数值每30秒出力一次CPU占用不到3%。 智能床垫呼吸心率监测：基于EMD的PVDF/压阻信号分离与Python实现
床垫这种产品做硬件的人一开始都觉得没什么门槛一块布、几片传感器、一个采集盒能出个在床/离床就交差了。真做到要输出心率和呼吸率这一步坑才开始一个接一个冒出来。我手头这个项目做了大概八个月从最早的PVDF压电薄膜方案到后来兼容压阻织物传感器中间换过三版算法。最后稳定下来的核心是用EMD算法把一路混在一起的体动信号拆成呼吸和心跳两条独立曲线再分别做频谱估计。这套心跳呼吸分离的流程我完整跑通了源代码也整理出来了下面连原理带代码一起讲清楚做PVDF传感器或者其他压阻人体体征传感器床垫的同行可以直接拿去改。先说清楚这套东西的适用边界。输入是单通道的体动信号采样率200Hz左右16位ADC输出是逐秒刷新的呼吸率次/分和心率bpm。传感器可以是PVDF压电薄膜也可以是压阻式导电织物、压阻橡胶、柔性应变片只要它能把胸腔的机械形变转成电信号就行。不挑传感器型号挑的是信号质量和采样链路的设计。适合做智能床垫、智能坐垫、婴儿监护垫的工程师看也适合做生理信号处理的学生参考代码全部是Python依赖只有numpy和scipy。1. 床垫体征监测的整体设计思路做这个项目之前我调研过好几条技术路线光电、雷达、压电、压阻都试过。雷达方案贵、功耗高光电方案在床垫上根本没法贴合最后落到压电和压阻这两条路上。它们的共同点是不主动发射任何东西纯被动感知机械形变结构简单成本可控塞进床垫里不影响睡感。麻烦的地方在于传感器出来的是一路信号呼吸和心跳全叠在里面怎么把这俩拆开就是整个项目的技术核心。1.1 PVDF和压阻式传感器在床垫场景下的优劣对比PVDF压电薄膜的工作原理是薄膜受力产生形变时内部偶极子取向变化两端出现电荷积累。它是典型的动态传感器输出的是电荷量形变越快输出越大静止不动输出就是零。这个特性用在床垫上有两个好处一是灵敏度高胸腔那点微米级的起伏它都能感应到二是天然隔直缓慢的体温漂移、床垫受压后的静态形变这类干扰不会跑进信号里。坏处也来自于同一个特性。它对低频响应差而呼吸正好是低频0.1到0.5Hz。如果电荷放大器的反馈电阻不够大呼吸信号会被削得很厉害后面EMD再厉害也救不回来。这一点后面会详细算。压阻式传感器的逻辑完全相反。它的电阻随压力变化需要外部激励恒压或恒流输出是电压变化能做到直流响应测静态压力没问题。缺点是灵敏度相对低而且本身有热漂移长时间工作基线会跑。好在你测的是呼吸和心跳这种交变成分把基线去掉就行。我把两者的对比整理成一张表方便选型时对照对比维度PVDF压电薄膜压阻式导电织物/压阻橡胶输出类型电荷需电荷放大器电阻变化需激励源低频响应差取决于Rf·Cf好可到DC灵敏度高中等静态压力无法测量可测量功耗极低无源需持续激励布线难度需屏蔽走线敏感相对不敏感成本中等偏高低适合场景呼吸心跳的精细波形在床/离床呼吸粗心率实际项目里我推荐的组合是PVDF做主传感压阻织物做辅助前者负责波形质量后者负责在床判断和姿态区分。如果预算只够一种做心率就选PVDF只做呼吸和在床检测选压阻就够了。1.2 呼吸和心跳为什么会粘在一起很多人第一次看床垫的原始波形会困惑明明波形挺干净的为什么频谱一出来是一堆乱七八糟的峰。原因在于两个信号本身的性质差异极大。呼吸的频率范围大概在0.1到0.5Hz对应每分钟6到30次。它的幅度大胸廓起伏带来的形变可能是心跳的10到50倍。心跳的频率在0.8到2.2Hz对应每分钟48到132次覆盖了绝大部分成年人的静息心率。幅度小但频率高。问题出在两个地方。一是频谱上虽然不重叠但真实信号的呼吸波形不是纯正弦它有二次、三次谐波0.3Hz的呼吸二次谐波就跑到0.6Hz三次谐波接近0.9Hz直接钻进心跳的频带里。二是呼吸过程中胸腔位置在变传感器和心脏之间的耦合强度也跟着变心跳信号的幅度被呼吸调制了这在时域上看就是心跳的包络跟着呼吸一起起伏。这两条决定了单纯用固定参数的带通滤波器拆不干净。滤波器只能按频率切切不掉谐波也处理不了幅度调制带来的畸变。这就是我最后转向EMD的直接原因。1.3 选EMD而不是固定带通滤波的三个理由EMD全称经验模态分解Empirical Mode Decomposition它的核心思想是任何复杂信号都可以分解成有限个本征模态函数IMF每个IMF是一个窄带分量有自己独立的瞬时频率和瞬时幅度。这个过程完全靠数据自己驱动不需要你预先指定任何基函数。选它有三个理由。第一它是自适应的。每个人的呼吸频率不一样同一个人入睡前后也不一样。固定带通滤波器的通带一旦设死遇到呼吸很慢的人比如8次/分0.13Hz滤波器通带如果从0.15Hz起呼吸主频就被削了。EMD不需要预设它自己会把这个频率上的成分单独拆成一个IMF。第二它能处理非平稳信号。人在睡眠中呼吸和心率的频率是缓慢漂移的傅里叶变换假设信号在窗口内平稳时间窗口一长就失真。EMD是基于局部极值点做的天然适应慢变。第三它对谐波的处理更自然。呼吸的三次谐波如果能量够大会被拆到单独的IMF里而不会污染心跳所在的IMF。当然这块也不是完美的后面会讲模态混叠这个坑。代价是计算量比滤波大以及端点效应、停止准则这些需要调。但在30到60秒的分析窗上跑PC端完全无压力嵌入式端需要做定点化和降采样优化。2. 硬件链路与采样参数怎么定算法再好前端信号烂了也是白搭。这一节把我在这块踩过的坑集中讲一下尤其是电荷放大器的低频截止频率计算我见过太多项目在这里翻车呼吸波形被削成一条直线还以为是算法问题。2.1 电荷放大器与低频截止的频率计算PVDF输出的是电荷不能直接接ADC中间必须有一个电荷放大器。基本结构是一个运放反馈回路里挂一个电容Cf和一个电阻Rf。输出电压为V_out -Q / Cf其中Q是传感器产生的电荷量。Cf决定了增益Cf越小增益越大。比如传感器在呼吸时产生1pC的电荷Cf取10nF输出就是0.1mV。这个量级很弱需要ADC前级再做一级放大或者把Cf降到1nF。真正关键的是Rf。它决定了低频截止频率f_c 1 / (2π · Rf · Cf)这个公式很直白Rf和Cf构成一个高通网络低于f_c的成分会被衰减。呼吸是低频信号所以f_c必须设得远低于呼吸的最低频率。我来算一下。假设Cf固定取10nF为了让0.1Hz的呼吸成分衰减不超过3dB即幅值保持在0.707以上f_c需要满足f_c ≤ 0.1 / sqrt(1/0.707² - 1) ≈ 0.1 Hz也就是说f_c至少要压到0.1Hz以下最好压到0.03Hz以下留足余量。取f_c 0.016Hz则Rf 1 / (2π × 0.016 × 10e-9) ≈ 1e9 Ω 1 GΩ1GΩ的电阻不是标准件通常用多颗高阻电阻串联或者用运放的反馈T型网络等效实现。我实际用的方案是两颗500MΩ的高阻电阻串联并联一个小电容做补偿。如果用100MΩ会怎样f_c 1/(2π × 1e8 × 10e-9) 0.159Hz。在0.1Hz处增益衰减为 0.1/sqrt(0.1² 0.159²) 0.53也就是-5.5dB。呼吸幅度直接掉一半而且不同呼吸频率衰减程度还不一样导致波形严重失真。这个坑我是真踩过当时排查了两天才发现是电阻选小了。注意电荷放大器的输入端是高阻节点PCB上必须做防护环guard ring走线尽量短否则漏电流会直接淹没信号。另外传感器的屏蔽层要接放大器地不能两端都接地形成地环路。2.2 采样率和位数的选择推演采样率的选择遵循奈奎斯特准则但工程上要留足余量。心跳最高频率取2.2Hz132bpm按照10倍过采样原则采样率至少22Hz。呼吸0.5Hz更没压力。但这里有一个容易被忽略的点工频干扰。如果采样率取50Hz奈奎斯特频率是25Hz。50Hz工频高于奈奎斯特频率会混叠。混叠到哪里采样频率是50Hz50Hz正好是采样率的整数倍混叠到0Hz表现为基线漂移。更糟的是如果采样率是49.9Hz50Hz会混叠到0.1Hz正好落在呼吸频带里这种情况你无论怎么滤波都救不回来因为它在数字域里和呼吸长得一模一样。结论采样率不能取50Hz。我推荐200Hz。这样奈奎斯特频率100Hz50Hz工频落在中间用数字陷波器可以精确干掉100Hz的二次谐波也能处理。那能不能取250Hz或者更高可以但没必要。采样率越高单窗口的数据点越多EMD的计算量线性增长。200Hz降采样到50Hz后30秒窗口是1500点EMD分解一层大概几毫秒整帧几十毫秒实时性完全够。ADC位数选16位。原因在于动态范围。呼吸信号幅度可能是心跳的20倍以上如果只有12位4096级留给心跳的分辨率就只剩200级波形会明显量化。16位是65536级心跳能分到3000级左右足够做谱分析了。2.3 预处理的四道工序拿到200Hz的原始数据后不能直接丢给EMD中间要做四件事顺序不能乱。第一道是工频陷波。用二阶IIR陷波器中心频率50HzQ值取30。如果采样率足够高比如400Hz以上再加一个100Hz的陷波。陷波器要用零相位滤波filtfilt否则会引入相位失真导致后面的峰间隔计算出现系统性偏差。第二道是滑动中值去基线。翻身、调整睡姿会让床垫压力分布发生阶跃式变化表现为一个大的台阶或缓慢漂移。用宽度2秒的滑动中值滤波估计基线再从原信号里减掉。这个方法比高通滤波好因为中值滤波对阶跃的响应是跟着走不会产生振铃。第三道是带通限幅。通带设0.05到10Hz就够了下限保呼吸上限保心跳的高次谐波并抑制高频噪声。10Hz以上的成分在这个应用里没有任何价值。第四道是重采样到50Hz。用多相滤波重采样resample_poly不要用简单抽取否则会引入混叠。50Hz下心跳2.2Hz有22倍采样足够。实操心得这四道工序里中值滤波的窗口宽度是唯一需要跟着床垫软硬度调的参数。硬床垫的翻身瞬变更快窗口要窄一些1.5秒软床垫形变释放慢窗口可以放到3秒。3. EMD分解的代码实现与关键细节这一节是全文的核心。我会把EMD的实现逻辑拆开讲包括筛分循环、端点延拓、停止准则然后是IMF的自动挑选。最后给出完整代码。代码我自己跑过很多次不是抄来的每个函数都是按实际需求写的。3.1 筛分循环与停止准则EMD的分解过程叫筛分sifting步骤是这样的找出信号的所有局部极大值点和局部极小值点用三次样条分别拟合上包络和下包络计算上下包络的均值曲线m(t)用原信号减去均值曲线得到h(t) x(t) - m(t)检查h(t)是否满足IMF的两个条件极值点数和过零点数相差不超过1上下包络均值处处为零如果不满足把h(t)当作新信号重复1到5如果满足把h(t)作为第一个IMF输出从原信号里减掉它对残差重复整个过程停止准则这里有个常见的做法分歧。Huang在1998年提出的是标准差准则SD准则即连续两次筛分结果的标准差小于阈值就停SD Σ[(h_{k-1}(t) - h_k(t))² / h_{k-1}²(t)]这个阈值一般取0.2到0.3。太大则IMF不够纯太小则迭代次数暴涨而且可能把信号筛成纯粹的调幅波丢失物理意义。另一个是极值点准则如果连续两次筛分后极值点数量和位置基本不变就停止。这个准则在实时的工程实现里更实用因为计算量小而且能避免过度筛分。我实际用的是混合策略SD阈值0.2作为主准则同时设最大迭代次数50次兜底。为什么要有兜底因为在某些噪声段信号几乎没有极值点筛分会陷入死循环。我见过一次没有兜底的实现在一段静默数据上卡了十几秒。3.2 端点效应、包络插值的坑EMD最臭名昭著的问题就是端点效应。原因很简单三次样条插值需要边界条件而信号的第一个极值点和最后一个极值点外面没有数据了。样条在这两段会剧烈发散导致包络在两端严重失真而且这个失真会随着筛分一层层向内传播甚至污染整个分解结果。解决办法有两类。一类是延拓法在信号两端人为补一些数据点让样条有足够的支撑。镜像延拓最简单把第一个极值点关于起点做镜像放到信号左边把最后一个极值点关于终点镜像放到右边。这个方法实现简单效果能接受我用的就是它。另一类是改进样条比如用B样条或者有理样条替代三次样条。效果更好但实现复杂对实时系统不划算。还有一个坑是极值点数量不足。当信号很短或者很平滑时极大值点可能只有一个甚至没有样条无法拟合。这时候必须直接返回把这个分量作为残差输出而不是硬拟合然后崩掉。代码里的len(max_idx) 2判断就是这个作用。另外三次样条插值对极值点的密集程度很敏感。如果信号里高频噪声多极值点会非常密集样条过拟合包络变成锯齿状分解出来的IMF全是噪声。所以预处理那一步的带通限幅非常必要把10Hz以上的噪声干掉极值点数量就正常了。3.3 用频率和能量准则自动挑选IMF分解出一堆IMF之后哪一个是呼吸哪一个是心跳最直接的办法是算每个IMF的平均频率。用零点穿越法估算统计一段信号里的过零点数量除以2再除以时长就是平均频率。这个方法对窄带信号很准实现也简单f_mean (过零点数 / 2) / 窗口时长呼吸IMF的频率应该在0.1到0.6Hz之间心跳IMF在0.8到2.5Hz之间。注意上限给到2.5Hz而不是2.2Hz是为了留余量防止心率高的时候比如运动后上床被漏掉。但光看频率还不够。有时候会有两个IMF都落在呼吸频段里这时候要看能量占比。能量占比定义为该IMF的能量与原始信号能量之比E_ratio Σ(imf²) / Σ(x²)同一个频段里能量占比最大的那个才是主成分其他的是谐波或噪声。这里有一个重要的经验EMD分解出来的第一个IMF最高频通常是噪声最后一个IMF最低频通常是趋势项中间的才是有效成分。我在挑选时会直接从IMF2开始扫跳过IMF1。选中呼吸IMF和心跳IMF后还有一步后处理对心跳IMF再做一次0.7到3.0Hz的带通。为什么因为EMD的模态混叠问题没有完全解决心跳IMF里常常残留呼吸的谐波成分。这一步带通能把这个残留压下去把心率估计的标准差从5bpm降到2bpm左右。3.4 完整可运行的Python源代码下面是完整代码包含预处理、EMD实现、IMF挑选、呼吸率心率估计的全部流程。直接复制到一个.py文件里就能跑。为了演示方便代码末尾生成了一个合成的床垫信号呼吸0.25Hz 心跳1.2Hz 噪声实际使用时把synthesize_mattress_signal换成你的采集数据即可。# -*- coding: utf-8 -*- 智能床垫体征算法EMD分解 心跳呼吸分离 传感器PVDF压电薄膜 / 压阻式导电织物 / 柔性应变片 输出呼吸率(次/分)、心率(bpm)、置信度标记 依赖numpy, scipy import numpy as np from scipy.signal import (butter, filtfilt, iirnotch, find_peaks, resample_poly, hilbert) from scipy.interpolate import CubicSpline from scipy.ndimage import median_filter # # 一、预处理 # def notch_filter(x, fs, f050.0, q30.0): 工频陷波零相位 if f0 fs / 2.0 * 0.95: return x b, a iirnotch(f0, q, fs) return filtfilt(b, a, x) def detrend_median(x, fs, win_sec2.0): 滑动中值去基线专治翻身台阶漂移 w int(win_sec * fs) if w % 2 0: w 1 if w 3: return x base median_filter(x, sizew, modenearest) return x - base def bandpass(x, fs, lo, hi, order4): 零相位带通 nyq fs / 2.0 lo_n max(lo / nyq, 1e-4) hi_n min(hi / nyq, 0.999) if lo_n hi_n: return x b, a butter(order, [lo_n, hi_n], btypeband) return filtfilt(b, a, x) def preprocess(x, fs_in, fs_out50.0): 完整预处理链陷波 - 去基线 - 带通 - 重采样 x np.asarray(x, dtypefloat) x x - np.mean(x) # 1) 工频陷波 x notch_filter(x, fs_in, 50.0, q30.0) if fs_in 250: x notch_filter(x, fs_in, 100.0, q30.0) # 2) 去基线漂移 x detrend_median(x, fs_in, win_sec2.0) # 3) 带通限幅 0.05 ~ 10 Hz x bandpass(x, fs_in, 0.05, 10.0, order4) # 4) 重采样到 50 Hz if abs(fs_in - fs_out) 1e-6: up, down _ratio(fs_out, fs_in) x resample_poly(x, up, down) return x, fs_out def _ratio(fs_out, fs_in, max_den200): 把 fs_out/fs_in 化简成整数比 from fractions import Fraction fr Fraction(fs_out / fs_in).limit_denominator(max_den) return fr.numerator, fr.denominator # # 二、EMD 核心实现 # def _extrema_idx(x): 定位局部极大值和极小值限制最小间隔避免噪声导致极值点爆炸 min_gap max(1, len(x) // 200) max_idx, _ find_peaks(x, distancemin_gap) min_idx, _ find_peaks(-x, distancemin_gap) return max_idx, min_idx def _spline_envelope(x, idx, n): 镜像延拓 三次样条拟合包络 if len(idx) 2: return None # 左端镜像把第一个极值点关于起点镜像 left_pos -idx[0] left_val x[idx[0]] # 右端镜像把最后一个极值点关于终点镜像 right_pos 2 * (n - 1) - idx[-1] right_val x[idx[-1]] pos np.concatenate(([left_pos], idx, [right_pos])) val np.concatenate(([left_val], x[idx], [right_val])) # 去掉重复位置避免样条报错 pos_u, uniq np.unique(pos, return_indexTrue) val_u val[uniq] if len(pos_u) 4: return None cs CubicSpline(pos_u, val_u, bc_typenatural) return cs(np.arange(n)) def _sift(x, sd_thresh0.2, max_iter50): 单次筛分返回一个IMF h x.astype(float).copy() for _ in range(max_iter): max_idx, min_idx _extrema_idx(h) if len(max_idx) 2 or len(min_idx) 2: return h, True # 极值点不足无法继续 upper _spline_envelope(h, max_idx, len(h)) lower _spline_envelope(h, min_idx, len(h)) if upper is None or lower is None: return h, True m 0.5 * (upper lower) h_new h - m denom np.sum(h ** 2) 1e-12 sd np.sum((h - h_new) ** 2) / denom h h_new if sd sd_thresh: break return h, False def emd(x, max_imf8, sd_thresh0.2, max_iter50): 经验模态分解返回IMF数组和残差 imfs [] r np.asarray(x, dtypefloat).copy() for _ in range(max_imf): max_idx, min_idx _extrema_idx(r) if len(max_idx) 2 or len(min_idx) 2: break imf, _ _sift(r, sd_thresh, max_iter) if np.sum(imf ** 2) 1e-12: break imfs.append(imf) r r - imf if len(imfs) 0: return np.zeros((0, len(x))), r return np.array(imfs), r # # 三、IMF 挑选 # def imf_mean_freq(imf, fs): 零点穿越法估平均频率 sign np.sign(imf) sign[sign 0] 1 zc np.sum(np.diff(sign) ! 0) return zc / 2.0 / (len(imf) / fs) def imf_energy_ratio(imf, x): return np.sum(imf ** 2) / (np.sum(x ** 2) 1e-12) def select_imf(imfs, x, fs, band, min_energy0.005): 在指定频带内挑选能量占比最大的IMF lo, hi band cand [] for i, imf in enumerate(imfs): f imf_mean_freq(imf, fs) e imf_energy_ratio(imf, x) if lo f hi and e min_energy: cand.append((i, f, e)) if not cand: return None, None, [] cand.sort(keylambda t: -t[2]) idx cand[0][0] return imfs[idx], idx, cand # # 四、频率估计 # def _parabolic_refine(freqs, amps, k): 抛物线插值细化谱峰位置 if 0 k len(amps) - 1: a, b, c amps[k - 1], amps[k], amps[k 1] denom a - 2 * b c if abs(denom) 1e-12: delta 0.5 * (a - c) / denom if abs(delta) 1.0: return freqs[k] delta * (freqs[1] - freqs[0]) return freqs[k] def spectrum_peak(x, fs, fmin, fmax, zero_pad8): 加汉宁窗 补零 抛物线插值返回频带内主峰频率 n len(x) if n int(fs / fmin * 2): return np.nan, 0.0 w np.hanning(n) xw (x - np.mean(x)) * w nfft 1 int(np.ceil(np.log2(n * zero_pad))) X np.abs(np.fft.rfft(xw, nfft)) freqs np.fft.rfftfreq(nfft, 1.0 / fs) mask (freqs fmin) (freqs fmax) if not np.any(mask): return np.nan, 0.0 sub_f freqs[mask] sub_a X[mask] k int(np.argmax(sub_a)) peak_f _parabolic_refine(sub_f, sub_a, k) # 谱峰突出度主峰 / 频带中位数 prom sub_a[k] / (np.median(sub_a) 1e-12) return peak_f, prom def peak_interval_rate(x, fs, fmin, fmax): 峰值间隔法返回中位频率(Hz)和间隔数 if fmax fmin: return np.nan, 0 min_dist int(fs / fmax * 0.7) pk, _ find_peaks(x, distancemax(1, min_dist)) if len(pk) 3: return np.nan, 0 ibi np.diff(pk) / fs lo, hi 1.0 / fmax, 1.0 / fmin ibi ibi[(ibi lo) (ibi hi)] if len(ibi) 2: return np.nan, 0 return 1.0 / np.median(ibi), len(ibi) # # 五、单帧分析主流程 # def analyze_frame(x_raw, fs_in200.0): 输入一段原始信号输出呼吸率、心率及置信信息 x, fs preprocess(x_raw, fs_in, fs_out50.0) n len(x) duration n / fs result { rr_bpm: np.nan, # 呼吸率 次/分 hr_bpm: np.nan, # 心率 bpm rr_hz: np.nan, hr_hz: np.nan, rr_spec_conf: 0.0, hr_spec_conf: 0.0, hr_energy: 0.0, rr_imf_idx: -1, hr_imf_idx: -1, quality: invalid } if duration 20: return result imfs, res emd(x, max_imf8, sd_thresh0.2, max_iter50) if imfs.shape[0] 2: return result # --- 呼吸IMF0.10 ~ 0.60 Hz resp_imf, resp_i, _ select_imf(imfs, x, fs, (0.10, 0.60), min_energy0.01) # --- 心跳IMF0.80 ~ 2.50 Hz card_imf, card_i, _ select_imf(imfs, x, fs, (0.80, 2.50), min_energy0.002) # ---------- 呼吸率 ---------- if resp_imf is not None: f1, c1 spectrum_peak(resp_imf, fs, 0.10, 0.60) f2, n2 peak_interval_rate(resp_imf, fs, 0.10, 0.60) if np.isfinite(f1) and np.isfinite(f2): if abs(f1 - f2) * 60 3.0: rr_hz 0.5 * (f1 f2) else: rr_hz f1 # 冲突时以谱峰为准 elif np.isfinite(f1): rr_hz f1 elif np.isfinite(f2): rr_hz f2 else: rr_hz np.nan if np.isfinite(rr_hz): result[rr_hz] rr_hz result[rr_bpm] rr_hz * 60.0 result[rr_spec_conf] c1 result[rr_imf_idx] int(resp_i) # ---------- 心率 ---------- if card_imf is not None: # 关键一步对心跳IMF再做一次带通压掉呼吸谐波残留 card_bp bandpass(card_imf, fs, 0.7, 3.0, order4) f1, c1 spectrum_peak(card_bp, fs, 0.80, 2.50) f2, n2 peak_interval_rate(card_bp, fs, 0.80, 2.50) if np.isfinite(f1) and np.isfinite(f2): if abs(f1 - f2) * 60 6.0: hr_hz 0.5 * (f1 f2) else: hr_hz f1 elif np.isfinite(f1): hr_hz f1 elif np.isfinite(f2): hr_hz f2 else: hr_hz np.nan if np.isfinite(hr_hz): result[hr_hz] hr_hz result[hr_bpm] hr_hz * 60.0 result[hr_spec_conf] c1 result[hr_imf_idx] int(card_i) result[hr_energy] imf_energy_ratio(card_imf, x) # ---------- 质量评估 ---------- ok_rr np.isfinite(result[rr_bpm]) and result[rr_spec_conf] 5 ok_hr (np.isfinite(result[hr_bpm]) and result[hr_spec_conf] 4 and result[hr_energy] 0.01) if ok_rr and ok_hr: result[quality] good elif ok_rr or ok_hr: result[quality] partial return result # # 六、演示合成一段床垫信号 # def synthesize_mattress_signal(fs200.0, dur60.0, rr0.25, hr1.20, seed7): rr: 呼吸频率Hz hr: 心率Hz rng np.random.default_rng(seed) t np.arange(0, dur, 1.0 / fs) # 呼吸含二次、三次谐波模拟真实非正弦呼吸 resp (1.00 * np.sin(2 * np.pi * rr * t) 0.22 * np.sin(2 * np.pi * 2 * rr * t 0.7) 0.08 * np.sin(2 * np.pi * 3 * rr * t 1.9)) # 心跳幅度被呼吸调制胸腔耦合强度随呼吸变化 mod 1.0 0.35 * np.sin(2 * np.pi * rr * t 0.4) card 0.045 * mod * np.sin(2 * np.pi * hr * t) # 工频 白噪声 hum 0.01 * np.sin(2 * np.pi * 50.0 * t) noise rng.normal(0, 0.004, len(t)) return resp card hum noise, fs if __name__ __main__: sig, fs synthesize_mattress_signal(fs200.0, dur60.0, rr0.25, hr1.20) out analyze_frame(sig, fs_infs) print(呼吸率: %.1f 次/分 (真值 15.0) % out[rr_bpm]) print(心率 : %.1f bpm (真值 72.0) % out[hr_bpm]) print(呼吸IMF编号: %d, 心跳IMF编号: %d % (out[rr_imf_idx], out[hr_imf_idx])) print(心跳IMF能量占比: %.4f % out[hr_energy]) print(质量标记: %s % out[quality])这段代码里我最想强调的是analyze_frame里对心跳IMF再带通那一步。很多人做完EMD就直接对IMF做FFT结果心率估计一直跳。原因就是模态混叠呼吸的三次谐波有时候会跑进心跳IMF里谱峰就被带偏了。加一道带通成本极低效果立竿见影。4. 呼吸率与心率的提取及实测验证算法跑通只是第一步能不能输出稳定的数值才是产品化的门槛。这一节讲我怎么做频率估计、怎么交叉校验、以及实测下来误差有多大。4.1 呼吸率的两种算法和交叉校验呼吸率的估计我用了两种方法并行。第一种是谱峰法。对呼吸IMF加汉宁窗补零到8倍长度做FFT在0.1到0.6Hz范围内找主峰然后用抛物线插值细化峰位。补零的作用是提高频率轴的采样密度让抛物线插值更有意义。不加补零的话30秒窗口的频率分辨率只有0.033Hz换算成呼吸率是2次/分误差太大。抛物线插值的效果值得说一下。假设真实呼吸是0.23Hz不加插值时谱峰只可能落在0.233Hz离散步长误差0.003Hz也就是0.18次/分。看起来还行但对于更短的分析窗比如20秒离散步长变成0.05Hz误差就上去了。加了抛物线插值能把峰位细化到离散步长的十分之一左右误差稳定在0.1次/分以内。第二种是峰值间隔法。直接在呼吸IMF上找波峰算相邻峰的间隔取中位数然后倒推频率。这个方法的好处是时间分辨率高能快速跟踪呼吸频率的变化。坏处是对噪声敏感如果IMF上有小的杂峰会被误判成呼吸峰。所以我在find_peaks里加了最小间隔限制最小间隔是最高呼吸频率对应周期的0.7倍。两种方法的结果做一致性检查如果差异小于3次/分取平均如果大于3次/分以谱峰法为准同时把这帧标记为低置信度。实测下来呼吸平稳时两者差异通常在1次/分以内翻身时可能出现较大分歧这时候的帧本来也不该被采信。4.2 心率提取从IMF到包络谱心率的估计比呼吸难原因有三信号弱、容易受运动干扰、频带和高频噪声相邻。主路径是对心跳IMF做带通后的谱峰估计和呼吸率的方法一样。但心率的频带更宽0.8到2.5Hz谐波干扰的概率更高。我在这里加了一个额外的可靠性判断谱峰突出度。定义为频带内主峰幅度除以频带内所有谱线幅度的中位数。一个干净的心跳信号这个比值通常在8以上如果低于4说明频带内能量分散没有明显的主频这帧数据不可信。另外还有一个交叉校验用的方法叫包络法。它的思路是当传感器和心脏之间隔了厚被子或者很软的床垫层时高频的搏动细节被平滑掉了心跳表现为一个窄带载波其幅度被心脏每次搏动调制。这时候对心跳IMF做希尔伯特变换取包络包络的起伏频率就是心率。提示包络法不是万能的。如果传感器直接耦合良好心跳在IMF上就是一个准正弦希尔伯特包络会接近常数包络法直接失效。我的做法是两种方法都算谁的谱峰突出度高就用谁。希尔伯特包络的实现很简单def envelope_rate(imf, fs, fmin0.8, fmax2.5): 希尔伯特包络谱法适用于心跳被平滑成载波的场景 env np.abs(hilbert(imf)) env env - np.mean(env) # 包络本身是低频信号重采样到25Hz足够 env_ds resample_poly(env, 1, 2) if fs 50 else env fs_ds fs / 2 if fs 50 else fs f, conf spectrum_peak(env_ds, fs_ds, fmin, fmax, zero_pad8) return f, conf注意包络法需要先降采样因为包络的有效带宽只有几Hz用50Hz采样做FFT是浪费而且高频噪声会污染包络谱。降到25Hz后2.5Hz的奈奎斯特余量还有5倍足够了。4.3 与参考设备的对比记录我找了12个人做对比测试其中8个是我同事4个是找的外部志愿者。测试方法受试者仰卧在装了PVDF传感器的样机上同时佩戴指夹式血氧仪测心率胸前绑一根应变带测呼吸。每人采集10分钟静息数据取躺下2分钟之后的8分钟数据做统计。数据分帧方式60秒窗口滑动步长1秒每帧输出一个估计值。把8分钟里所有good质量标记的帧取中位数作为该受试者的最终结果。统计下来的结果指标平均绝对误差95分位误差呼吸率1.3 次/分2.8 次/分心率2.1 bpm4.6 bpm这个精度对于非医疗级的睡眠监测产品是够用的。作为对照我最早用固定带通滤波那一版心率的平均绝对误差是5.4 bpm95分位误差超过11 bpm完全不能用。还有几个观察值得记录。第一体重轻的人BMI低于19心率误差明显大因为胸腔传导到床垫的机械能量小心跳IMF的能量占比往往低于1.5%接近我设的置信度门槛。第二呼吸率在浅睡期误差会略大因为呼吸变得不规则波形不是周期性的谱峰会变宽。第三翻身动作后大约5到8秒内结果不可信我在实际产品上用加速度计触发一个运动屏蔽窗口这段时间的输出保持上一帧的值不变。5. 常见问题排查与调参经验谈写了这么多原理和代码最后这部分是我最想分享的。因为上面那些东西你看书看论文都能找到但下面这些坑只有真做过床垫项目的人才知道。5.1 问题速查表我把这八个月里遇到过的典型问题整理成表遇到类似现象可以先在这里对一下。现象可能原因排查手段解决方式呼吸率恒为某个固定值不变化呼吸信号被削平算法在拟合噪声画呼吸IMF时域波形看是否有起伏检查电荷放大器Rf是否足够大f_c是否低于0.03Hz心率忽高忽低跳动超过15bpm心跳IMF能量弱谱峰被噪声主导打印每帧心跳IMF的能量占比提高置信度门槛能量低于2%直接丢弃该帧呼吸率是真实值的2倍或3倍呼吸谐波被选为呼吸IMF打印所有IMF的中心频率和能量在候选IMF里强制选最低频的那个或加谐波约束每帧结果差异大重启后不一致端点效应导致分解结果不确定对比不同起点切窗的分解结果增加镜像延拓帧间重叠50%多帧取中位数夜间某个时段数据全废翻身或离床看原始信号RMS是否有突增加运动屏蔽窗口RMS超阈值时冻结输出CPU占用率高实时性跟不上EMD在高采样率下计算量大统计单帧耗时先降采样到25Hz再做EMD心跳2.5Hz仍有10倍余量心率偏高约等于呼吸率的整数倍模态混叠呼吸谐波混进心跳IMF对比呼吸IMF和心跳IMF的时域相关性心跳IMF后接0.7-3.0Hz带通不同人之间精度差异巨大传感器位置和耦合强度不一致记录每个受试者的心跳IMF能量占比产品化时做一次个体校准记录基准能量水平5.2 几个我自己踩出来的调参心得第一个心得关于EMD的分解层数。max_imf参数我设的是8但实际有效的通常只有4到5层。设太大的问题是后面几层全是趋势项白算。设太小的问题是心跳IMF可能出不来。我的经验是先设12层跑一次把每层的中心频率和能量占比打出来看清楚心跳落在第几层然后再把max_imf收敛到那个层数加2。不同床垫不一样硬床垫心跳IMF往往在第3层软床垫因为高频被吸收会落到第2层。第二个心得关于分析窗口长度。论文里常见的是1秒到5秒的窗口那是给心电信号用的。床垫信号里心跳太弱短窗口的谱分辨率根本不够。我试过5秒、10秒、20秒、30秒、60秒最后定在60秒窗口加1秒步长。60秒的频率分辨率是0.017Hz配合抛物线插值能到0.002Hz对应心率误差0.12bpm完全够。代价是心率变化有大约30秒的滞后因为窗口的一半都在过去对于睡眠监测这种慢变场景可以接受。如果要跟踪运动后的心率恢复得用卡尔曼滤波做平滑和延迟补偿。第三个心得关于置信度的设计。一开始我没有置信度算法什么都输出结果用户看到心率从60突然跳到120又跳回来体验极差。后来加了三个门槛谱峰突出度、IMF能量占比、呼吸和心率的合理性检查心率必须大于呼吸率的2倍且心率在40到200之间。三个都过了才输出否则保持上一帧的值。这个改动让数值的平滑度提升了一个数量级用户侧看起来就是偶尔卡一下而不是乱跳。第四个心得关于模态混叠的处理。这是个老大难问题标准EMD解决不了。我试过EEMD加白噪声的集合平均效果好但计算量翻十倍实时系统扛不住。后来用了一个折中方案在预处理阶段把10Hz以上的噪声干掉让极值点分布更均匀模态混叠的概率就明显降低了。另外就是前面反复强调的对心跳IMF再做一次带通这个是最实用的补丁。实操心得如果要进一步提升心跳的分离质量可以试试带掩膜信号的EMDMasking EMD也就是在分解前往信号里叠加一个已知频率的高频正弦把这个正弦的频率设在心跳频带之外让心跳IMF的特征更明显分解完再减掉。这个方法比EEMD轻量很多我在一块低功耗MCU上跑过单帧耗时大概80ms。最后再提一个设备端的优化方向。EMD的主要开销在三次样条插值和多次筛分这两块都可以做定点化。样条的系数计算用查表加线性插值近似误差在1%以内但速度能快3倍。筛分次数限制在8次以内对结果的影响很小。我把这些优化做完之后200Hz采样、60秒窗口的整帧处理在Cortex-M4上大概需要200ms左右如果只要心跳和呼吸率这两个数值每30秒出力一次CPU占用不到3%。