2. 心率变异性分析

任务目标

从一段PPG（光电容积脉搏波）信号中提取心率变异性（HRV）指标，用于评估自主神经功能。本案例聚焦于时域HRV指标：平均NN间期（mean NN）、SDNN（NN间期标准差）和RMSSD（相邻NN间期差值的均方根），并输出用于后续分析的统计结果。

环境与依赖

需要Python 3.8+，以及numpy、scipy库。若需绘制波形图，可安装matplotlib（本案例不强制要求）。

实现思路

PPG信号中的每个脉搏波峰值对应一次心跳。通过峰值检测获得心跳时间点，计算相邻心跳间隔（NN间期），进而得到HRV时域指标。本案例使用自适应阈值峰值检测，避免因信号幅值漂移导致的漏检或误检。

关键完整代码

import numpy as np
from scipy.signal import find_peaks

def detect_heartbeats(ppg_signal, fs, min_distance=0.3):
    """
    从PPG信号中检测心跳峰值
    ppg_signal: 一维PPG信号数组
    fs: 采样率(Hz)
    min_distance: 两个心跳之间的最小时间间隔(秒)
    返回: 峰值索引数组
    """
    # 自适应阈值：使用滑动窗口的局部均值作为基线
    window = int(fs * 2)  # 2秒窗口
    baseline = np.convolve(ppg_signal, np.ones(window)/window, mode='same')
    corrected = ppg_signal - baseline
    
    # 最小峰值距离（转换为采样点数）
    min_dist_samples = int(min_distance * fs)
    
    # 峰值检测：阈值为修正信号标准差的0.5倍
    threshold = 0.5 * np.std(corrected)
    peaks, _ = find_peaks(corrected, height=threshold, distance=min_dist_samples)
    return peaks

def compute_hrv_metrics(ppg_signal, fs):
    """
    从PPG信号计算HRV时域指标
    返回: (mean_nn_ms, sdnn_ms, rmssd_ms)
    """
    # 1. 检测心跳峰值
    peaks = detect_heartbeats(ppg_signal, fs)
    if len(peaks) < 3:
        raise ValueError("检测到的心跳次数不足，无法计算HRV指标")
    
    # 2. 计算NN间期（单位：毫秒）
    peak_times = peaks / fs  # 转换为秒
    nn_intervals = np.diff(peak_times) * 1000  # 转换为毫秒
    
    # 3. 过滤异常间期（如因运动伪迹导致的异常值）
    # 使用中值偏差(MAD)法剔除超出3倍MAD的间期
    median_nn = np.median(nn_intervals)
    mad = np.median(np.abs(nn_intervals - median_nn))
    if mad > 0:
        lower = median_nn - 3 * 1.4826 * mad
        upper = median_nn + 3 * 1.4826 * mad
        nn_intervals = nn_intervals[(nn_intervals >= lower) & (nn_intervals <= upper)]
    
    if len(nn_intervals) < 2:
        raise ValueError("过滤后NN间期不足，无法计算HRV指标")
    
    # 4. 计算时域HRV指标
    mean_nn = np.mean(nn_intervals)
    sdnn = np.std(nn_intervals, ddof=1)  # 样本标准差
    diff_nn = np.diff(nn_intervals)
    rmssd = np.sqrt(np.mean(diff_nn ** 2))
    
    return mean_nn, sdnn, rmssd

# ===== 最小验证示例 =====
if __name__ == "__main__":
    # 生成模拟PPG信号（60秒，采样率100Hz）
    fs = 100
    t = np.arange(0, 60, 1/fs)
    
    # 模拟心率：基础心率70bpm + 呼吸性窦性心律不齐
    hr_base = 70 / 60  # Hz
    hr_variation = 0.03 * np.sin(2 * np.pi * 0.2 * t)  # 呼吸频率0.2Hz
    instantaneous_hr = hr_base + hr_variation
    # 积分得到心跳相位
    phase = 2 * np.pi * np.cumsum(instantaneous_hr) / fs
    # 构造PPG波形（每个心跳一个高斯脉冲）
    ppg = np.zeros_like(t)
    for i in range(len(t)):
        # 简化：使用正弦波叠加模拟脉搏波
        ppg[i] = 0.5 + 0.4 * np.sin(phase[i]) + 0.1 * np.sin(2 * phase[i])
    
    # 添加少量噪声
    np.random.seed(42)
    ppg += np.random.normal(0, 0.02, len(t))
    
    # 计算HRV指标
    mean_nn, sdnn, rmssd = compute_hrv_metrics(ppg, fs)
    
    print(f"平均NN间期: {mean_nn:.1f} ms")
    print(f"SDNN: {sdnn:.1f} ms")
    print(f"RMSSD: {rmssd:.1f} ms")

运行与验证

将上述代码保存为hrv_analysis.py，在终端运行：

python hrv_analysis.py

预期输出范围（基于模拟信号参数）：平均NN间期约850～870 ms（对应心率约70 bpm），SDNN约15～30 ms，RMSSD约15～35 ms。由于模拟信号包含呼吸性窦性心律不齐，RMSSD应明显大于0。

验证通过条件：程序无报错，输出的三个指标均为正数，且平均NN间期对应的平均心率在68～72 bpm范围内（即NN间期在833～882 ms之间）。若使用真实PPG数据，需确保信号长度至少包含30秒的有效心跳数据，采样率不低于50 Hz。

常见问题

• 峰值检测不准：若信号噪声较大，可先进行带通滤波（0.5～5 Hz）再传入检测函数。
• NN间期异常：运动伪迹会导致间期突变，MAD过滤可剔除大部分异常值；若过滤后数据不足，需延长记录时间。
• 采样率影响：采样率过低（<25 Hz）会降低峰值定位精度，直接影响RMSSD的准确性。