企业官网建设流程全解析

床垫这种产品，做硬件的人一开始都觉得没什么门槛：一块布、几片传感器、一个采集盒，能出个"在床/离床"就交差了。真做到要输出心率和呼吸率这一步，坑才开始一个接一个冒出来。我手头这个项目做了大概八个月，从最早的PVDF压电薄膜方案，到后来兼容压阻织物传感器，中间换过三版算法。最后稳定下来的核心，是用EMD算法把一路混在一起的体动信号拆成呼吸和心跳两条独立曲线，再分别做频谱估计。这套心跳呼吸分离的流程我完整跑通了，源代码也整理出来了，下面连原理带代码一起讲清楚，做PVDF传感器或者其他压阻人体体征传感器床垫的同行可以直接拿去改。

先说清楚这套东西的适用边界。输入是单通道的体动信号，采样率200Hz左右，16位ADC；输出是逐秒刷新的呼吸率（次/分）和心率（bpm）。传感器可以是PVDF压电薄膜，也可以是压阻式导电织物、压阻橡胶、柔性应变片，只要它能把胸腔的机械形变转成电信号就行。不挑传感器型号，挑的是信号质量和采样链路的设计。适合做智能床垫、智能坐垫、婴儿监护垫的工程师看，也适合做生理信号处理的学生参考，代码全部是Python，依赖只有numpy和scipy。

1. 床垫体征监测的整体设计思路

做这个项目之前我调研过好几条技术路线，光电、雷达、压电、压阻都试过。雷达方案贵、功耗高，光电方案在床垫上根本没法贴合，最后落到压电和压阻这两条路上。它们的共同点是：不主动发射任何东西，纯被动感知机械形变，结构简单，成本可控，塞进床垫里不影响睡感。麻烦的地方在于，传感器出来的是一路信号，呼吸和心跳全叠在里面，怎么把这俩拆开就是整个项目的技术核心。

1.1 PVDF和压阻式传感器在床垫场景下的优劣对比

PVDF压电薄膜的工作原理是：薄膜受力产生形变时，内部偶极子取向变化，两端出现电荷积累。它是典型的动态传感器，输出的是电荷量，形变越快输出越大，静止不动输出就是零。这个特性用在床垫上有两个好处，一是灵敏度高，胸腔那点微米级的起伏它都能感应到；二是天然隔直，缓慢的体温漂移、床垫受压后的静态形变这类干扰不会跑进信号里。

坏处也来自于同一个特性。它对低频响应差，而呼吸正好是低频（0.1到0.5Hz）。如果电荷放大器的反馈电阻不够大，呼吸信号会被削得很厉害，后面EMD再厉害也救不回来。这一点后面会详细算。

压阻式传感器的逻辑完全相反。它的电阻随压力变化，需要外部激励（恒压或恒流），输出是电压变化，能做到直流响应，测静态压力没问题。缺点是灵敏度相对低，而且本身有热漂移，长时间工作基线会跑。好在你测的是呼吸和心跳这种交变成分，把基线去掉就行。

我把两者的对比整理成一张表，方便选型时对照：

实际项目里我推荐的组合是：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.2Hz（132bpm），按照10倍过采样原则，采样率至少22Hz。呼吸0.5Hz，更没压力。但这里有一个容易被忽略的点：工频干扰。

如果采样率取50Hz，奈奎斯特频率是25Hz。50Hz工频高于奈奎斯特频率，会混叠。混叠到哪里？采样频率是50Hz，50Hz正好是采样率的整数倍，混叠到0Hz，表现为基线漂移。更糟的是，如果采样率是49.9Hz，50Hz会混叠到0.1Hz，正好落在呼吸频带里，这种情况你无论怎么滤波都救不回来，因为它在数字域里和呼吸长得一模一样。

结论：采样率不能取50Hz。我推荐200Hz。这样奈奎斯特频率100Hz，50Hz工频落在中间，用数字陷波器可以精确干掉，100Hz的二次谐波也能处理。

那能不能取250Hz或者更高？可以，但没必要。采样率越高，单窗口的数据点越多，EMD的计算量线性增长。200Hz降采样到50Hz后，30秒窗口是1500点，EMD分解一层大概几毫秒，整帧几十毫秒，实时性完全够。

ADC位数选16位。原因在于动态范围。呼吸信号幅度可能是心跳的20倍以上，如果只有12位（4096级），留给心跳的分辨率就只剩200级，波形会明显量化。16位是65536级，心跳能分到3000级左右，足够做谱分析了。

2.3 预处理的四道工序

拿到200Hz的原始数据后，不能直接丢给EMD，中间要做四件事，顺序不能乱。

第一道是工频陷波。用二阶IIR陷波器，中心频率50Hz，Q值取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），步骤是这样的：

1. 找出信号的所有局部极大值点和局部极小值点
2. 用三次样条分别拟合上包络和下包络
3. 计算上下包络的均值曲线m(t)
4. 用原信号减去均值曲线，得到h(t) = x(t) - m(t)
5. 检查h(t)是否满足IMF的两个条件：极值点数和过零点数相差不超过1；上下包络均值处处为零
6. 如果不满足，把h(t)当作新信号重复1到5
7. 如果满足，把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, f0=50.0, q=30.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_sec=2.0): """滑动中值去基线，专治翻身台阶漂移""" w = int(win_sec * fs) if w % 2 == 0: w += 1 if w < 3: return x base = median_filter(x, size=w, mode='nearest') return x - base def bandpass(x, fs, lo, hi, order=4): """零相位带通""" 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], btype='band') return filtfilt(b, a, x) def preprocess(x, fs_in, fs_out=50.0): """完整预处理链：陷波 -> 去基线 -> 带通 -> 重采样""" x = np.asarray(x, dtype=float) x = x - np.mean(x) # 1) 工频陷波 x = notch_filter(x, fs_in, 50.0, q=30.0) if fs_in > 250: x = notch_filter(x, fs_in, 100.0, q=30.0) # 2) 去基线漂移 x = detrend_median(x, fs_in, win_sec=2.0) # 3) 带通限幅 0.05 ~ 10 Hz x = bandpass(x, fs_in, 0.05, 10.0, order=4) # 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_den=200): """把 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, distance=min_gap) min_idx, _ = find_peaks(-x, distance=min_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_index=True) val_u = val[uniq] if len(pos_u) < 4: return None cs = CubicSpline(pos_u, val_u, bc_type='natural') return cs(np.arange(n)) def _sift(x, sd_thresh=0.2, max_iter=50): """单次筛分，返回一个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_imf=8, sd_thresh=0.2, max_iter=50): """经验模态分解，返回IMF数组和残差""" imfs = [] r = np.asarray(x, dtype=float).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_energy=0.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(key=lambda 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_pad=8): """加汉宁窗 + 补零 + 抛物线插值，返回频带内主峰频率""" 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, distance=max(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_in=200.0): """输入一段原始信号，输出呼吸率、心率及置信信息""" x, fs = preprocess(x_raw, fs_in, fs_out=50.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_imf=8, sd_thresh=0.2, max_iter=50) if imfs.shape[0] < 2: return result # --- 呼吸IMF：0.10 ~ 0.60 Hz resp_imf, resp_i, _ = select_imf(imfs, x, fs, (0.10, 0.60), min_energy=0.01) # --- 心跳IMF：0.80 ~ 2.50 Hz card_imf, card_i, _ = select_imf(imfs, x, fs, (0.80, 2.50), min_energy=0.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, order=4) 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(fs=200.0, dur=60.0, rr=0.25, hr=1.20, seed=7): """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(fs=200.0, dur=60.0, rr=0.25, hr=1.20) out = analyze_frame(sig, fs_in=fs) 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, fmin=0.8, fmax=2.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_pad=8) return f, conf

注意包络法需要先降采样，因为包络的有效带宽只有几Hz，用50Hz采样做FFT是浪费，而且高频噪声会污染包络谱。降到25Hz后，2.5Hz的奈奎斯特余量还有5倍，足够了。

4.3 与参考设备的对比记录

我找了12个人做对比测试，其中8个是我同事，4个是找的外部志愿者。测试方法：受试者仰卧在装了PVDF传感器的样机上，同时佩戴指夹式血氧仪测心率，胸前绑一根应变带测呼吸。每人采集10分钟静息数据，取躺下2分钟之后的8分钟数据做统计。

数据分帧方式：60秒窗口，滑动步长1秒，每帧输出一个估计值。把8分钟里所有good质量标记的帧取中位数作为该受试者的最终结果。

统计下来的结果：

这个精度对于非医疗级的睡眠监测产品是够用的。作为对照，我最早用固定带通滤波那一版，心率的平均绝对误差是5.4 bpm，95分位误差超过11 bpm，完全不能用。

还有几个观察值得记录。第一，体重轻的人（BMI低于19）心率误差明显大，因为胸腔传导到床垫的机械能量小，心跳IMF的能量占比往往低于1.5%，接近我设的置信度门槛。第二，呼吸率在浅睡期误差会略大，因为呼吸变得不规则，波形不是周期性的，谱峰会变宽。第三，翻身动作后大约5到8秒内结果不可信，我在实际产品上用加速度计触发一个"运动屏蔽窗口"，这段时间的输出保持上一帧的值不变。

5. 常见问题排查与调参经验谈

写了这么多原理和代码，最后这部分是我最想分享的。因为上面那些东西，你看书看论文都能找到，但下面这些坑，只有真做过床垫项目的人才知道。

5.1 问题速查表

我把这八个月里遇到过的典型问题整理成表，遇到类似现象可以先在这里对一下。

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再做一次带通，这个是最实用的补丁。

实操心得：如果要进一步提升心跳的分离质量，可以试试带掩膜信号的EMD（Masking EMD），也就是在分解前往信号里叠加一个已知频率的高频正弦，把这个正弦的频率设在心跳频带之外，让心跳IMF的特征更明显，分解完再减掉。这个方法比EEMD轻量很多，我在一块低功耗MCU上跑过，单帧耗时大概80ms。

最后再提一个设备端的优化方向。EMD的主要开销在三次样条插值和多次筛分，这两块都可以做定点化。样条的系数计算用查表加线性插值近似，误差在1%以内，但速度能快3倍。筛分次数限制在8次以内，对结果的影响很小。我把这些优化做完之后，200Hz采样、60秒窗口的整帧处理在Cortex-M4上大概需要200ms左右，如果只要心跳和呼吸率这两个数值，每30秒出力一次，CPU占用不到3%。