Matlab HRV特征提取工具箱：从RR间期到非线性指标的全流程解析

简介：Matlab环境下的HRV（心率变异性）特征提取与非线性计算工具包，面向生物医学工程、运动生理学及心理学领域的研究人员和学生，用于从心电信号中提取RR间期并计算多种HRV指标。资源核心围绕非线性动力学分析展开，覆盖样本熵（SampEn）、近似熵、模糊熵以及分形维数（DFA）等算法，同时包含时间域HRV、频率域HRV、Poincare散点图、滑动窗口处理等模块，并提供RRI提取、异常值替换、小波滤波等预处理脚本，可支撑从原始数据到特征输出的完整流程。压缩包共46个文件，以42个M脚本为主，辅以4个XLS示例数据，整体大小仅299KB，结构清晰、模块化设计便于二次开发。目前已有678人学习下载。代码中配备main.m入口和多个可独立调用的函数文件，并附带心率、HRV-EEG等示例表格，适合具备一定Matlab基础、希望快速搭建或扩展HRV分析流程的用户直接上手。

1. 从RR间期到非线性指标：这套Matlab HRV工具箱为什么值得拆开看

心率变异性分析在外行看来就是算个标准差，但真正跑过数据的人都知道，线性指标在疲劳、压力、睡眠分期这些场景里的区分度往往不够，反而是样本熵、模糊熵、Poincare几何这类非线性特征能拉开差距。这套Matlab源码把HRV特征提取做成了完整闭环：原始ECG或IBI序列进来，经过小波滤波、异常值替换，一路算到时域、频域、时频域和非线性特征，最后统一导出Excel。它按函数拆分，main.m里能看到全部调用次序，适合两类人：一是做心率变异性研究的生医工程方向学生，需要一套能复现、能改参数的特征工程基线；二是要在Matlab环境里快速验证HRV特征有效性的算法工程师，把其中某个函数抽到自己管线里也不费劲。

2. 预处理链路：R波定位、IBI序列与异常值替换的工程细节

2.1 从ECG到RR间期：Extract_RRI与preProcessIBI的分工

HRV分析的质量上限在预处理阶段就决定了，后面算再多的熵也补救不回来。这套代码里Extract_RRI负责从原始ECG波形中定位R波峰值，输出相邻R峰之间的时间间隔序列；preProcessIBI则接受已经存在的RR间期或IBI序列，做重采样、缺失值插值和格式规范化。分工的含义是输入既可以是原始心电信号，也可以是可穿戴设备导出的间期数据——后一种情况在实验数据里非常常见，很多手环手表导出的心率值本身就是厂商滤波后的结果，直接喂给特征函数反而会引入未知延迟。

R波定位的常见做法是基于自适应阈值的Pan-Tompkins变体：先对ECG做带通滤波，再通过移动窗口积分突出QRS能量，最后用回看窗口确认R峰，避免把T波误检成R波。实际使用时要注意采样率，低于500Hz的数据定位误差会直接反映在RR间期抖动上，在计算RMSSD和样本熵时被显著放大。preProcessIBI里还需要处理两类典型伪差：漏检导致的间期翻倍，以及异位搏动造成的短间期。前者表现为RR间期突然变为两倍左右，后者表现为局部间期骤降，这两类都必须标记出来后进入异常值替换环节。

2.2 replaceOutliers：用MAD滑动窗口替换而不是直接删除

RR间期序列是等间隔采样的连续变量，直接删掉异常点会破坏时间连续性，后续做频域分析时插值位置会引入额外的频谱泄漏。replaceOutliers采用的做法是滑动窗口中位数加中位数绝对偏差（MAD）判定，再原位替换，保证序列长度不变。

% 基于滑动中位数与MAD的RR间期异常值替换 function rri_clean = replaceOutliers(rri, threshold) % 输入: rri - 原始RR间期序列, 单位ms % threshold - MAD倍数阈值, 心率变异性研究中常取3~5 win = 21; % 滑动窗口长度, 约20秒心跳数 med_rri = movmedian(rri, win, 'omitnan'); % 窗口内中位数 mad_rri = movmad(rri, win, 'omitnan'); % 窗口内中位数绝对偏差 dev = abs(rri - med_rri); idx_rep = (dev > threshold * mad_rri) | isnan(rri); rri_clean = rri; rri_clean(idx_rep) = med_rri(idx_rep); % 用中位数替换, 不删除 fprintf('替换异常RR间期 %d 个, 占比 %.2f%%\n', ... sum(idx_rep), 100 * sum(idx_rep) / length(rri)); end

这段代码的要点在三个参数上。window取21对应约20秒的滑动窗口，足够覆盖呼吸性窦性心律不齐的周期，又不会把真实的慢变趋势误判成异常；threshold取3时偏激进，适合运动场景下的强伪差数据，取5则更保守，适合静息态研究；用中位数而不是均值做基准，是因为中位数对单点异常天然鲁棒，不会被离群值本身拉偏。替换完成后建议打印占比，如果超过5%，说明原始信号质量或R波检测参数有问题，先回头处理，而不是继续往下算。

替换后的序列还要检查是否有连续多个点被替换，那种情况往往是局部信号丢失，需要回到ECG段重新检测，而不只是替换。

2.3 wavelet_filter与wavelet.m：小波去噪的层数与小波基选择

wavelet_filter负责对ECG或IBI序列做小波去噪，wavelet.m是底层的小波分解与重构实现。ECG去噪的常规选择是db4或sym8小波，分解层数根据采样率决定，目标是只抑制基线漂移和高频肌电干扰，保留QRS能量。

小波去噪最常见的坑是把分解层数设得过大，导致QRS波群被当成细节系数置零，滤波后的ECG里R峰幅值缩水，R波检测漏检率上升。另一个容易忽略的点是软阈值与硬阈值的差异：硬阈值保留峰值但会产生震荡伪迹，软阈值平滑但对R峰幅值有压缩。我这边的做法是在wavelet.m里对细节系数用软阈值、对近似系数不动，只做基线校正，这样后续R波检测的形态信息损失最小。参数调整后，建议用滤波前后RR间期序列的相关性做校验，相关系数低于0.95就该回查阈值设置。

3. 线性HRV特征：时域、频域与联合时频域的计算口径

3.1 timeDomainHRV：SDNN、RMSSD与pNN50的边界条件

时域指标是HRV特征提取的基础盘，timeDomainHRV输出的通常是三组数值：SDNN反映整体变异性，RMSSD反映迷走神经介导的快速变化，pNN50是相邻间期差超过50ms的比例。计算口径上有个容易出错的地方——SDNN既可以是24小时整体标准差，也可以是短时静息段的标准差，两者在文献里都叫SDNN，但参考范围完全不同。这套代码默认按短时分析处理，5分钟静息段的SDNN正常范围在30~60ms，低于20ms需要警惕。

% 时域HRV特征计算核心片段 function feat = timeDomainHRV(rri) diff_rri = diff(rri); % 相邻RR间期差值 feat.SDNN = std(rri, 'omitnan'); % 整体标准差 feat.RMSSD = sqrt(mean(diff_rri.^2, 'omitnan')); % 差值的均方根 feat.pNN50 = 100 * sum(abs(diff_rri) > 50) / numel(diff_rri); feat.meanHR = 60000 / mean(rri); % 平均心率值, 单位bpm end

这里有一个常被忽略的细节：pNN50的分母是差分个数而非间期个数，序列长度为N时差分只有N-1个，样本量小的时候这个偏差会明显影响百分比数值。meanHR用60000除以平均间期得到，注意如果rri单位不是毫秒，这里要同步调整。RMSSD对异常值极度敏感，一个未替换的伪差就能让RMSSD翻倍，所以时域计算必须在replaceOutliers之后执行，顺序颠倒就没有意义了。

3.2 freqDomainHRV：功率谱估计与LF/HF频带划分

频域特征把RR间期序列变换到频率域，用VLF、LF、HF三个频带的功率以及LF/HF比值描述自主神经活动。freqDomainHRV里功率谱估计通常有两种实现：经典周期图法或AR模型法。周期图法直接对去趋势后的间期序列做FFT，窗函数选汉宁窗，零填充到512点以上以获得平滑谱线；AR模型法阶数取16~20，分辨率更高但阶数敏感，不同阶数下LF功率可能差20%以上。

频带划分标准要跟文献对齐：VLF为0.003~0.04Hz，LF为0.04~0.15Hz，HF为0.15~0.4Hz。LF/HF比值在静息态解读为交感与迷走张力的相对平衡，但要注意呼吸频率低于0.15Hz时，呼吸性窦性心律不齐的能量会落入LF频带，此时比值解释需要谨慎。计算HF功率时建议同步输出呼吸率作为参考，我在实际项目中遇到过低频呼吸让LF/HF虚高的情况，不结合呼吸数据根本无法判断。

频域计算前必须对间期序列做去趋势处理，直接用原始RR间期做FFT，低频段会被线性趋势主导，VLF功率失真严重。freqDomainHRV里如果保留过一次差分或多项式去趋势的选项，优先用二阶多项式拟合去除慢漂移。

3.3 timeFreqHRV与slidingWindow：非平稳段的时频联合分析

时长超过几分钟的心率数据往往包含明显的非平稳成分，直接算整段频谱会把瞬时变化平均掉。timeFreqHRV配合slidingWindow做短时傅里叶变换，窗口长度取60~120秒、步长取30秒，逐窗口输出LF和HF功率序列，得到的是频率特征随时间的变化轨迹。

滑动窗口的参数选择直接影响结果形态。窗口太短（低于30秒）时LF频带只有不到2个完整周期，功率估计方差极大；窗口太长则时间分辨率下降，相邻窗口结果几乎重复。实际项目中我会按心率数据的用途来定：运动恢复分析用60秒窗口配30秒步长，睡眠分期用120秒窗口配60秒步长。每个窗口内先做异常值替换再做去趋势，避免单个伪差污染整个窗口的频谱。时频分析输出的特征矩阵可以直接作为后续分类模型的输入，每一行对应一个时间窗的LF、HF、LF/HF和总功率。

4. 非线性特征：样本熵、近似熵、模糊熵与Poincare几何

4.1 SampEn、ApEn与FuzzyEn：三种熵指标的适用边界

非线性特征是这套代码里最有区分度的部分。样本熵（SampEn）衡量时间序列的不规则性，对数据长度不敏感，是当前HRV研究的默认选择；近似熵（ApEn）实现较早但存在对参数m和r的依赖偏置，短序列结果一致性差；模糊熵（FuzzyEn）用指数函数替代阶跃判定，抗噪能力强，对小样本更稳定。三个函数里都有m（嵌入维数）和r（相似容差）两个参数，m通常取2，r取序列标准差的0.15~0.25倍。

% 样本熵核心计算: 统计模板向量匹配对数 function se = Sample_entropy(x, m, r) N = length(x); phi_m = phi(x, m, r); % m维模板匹配对数 phi_m1 = phi(x, m + 1, r); % m+1维模板匹配对数 se = -log(phi_m1 / phi_m); % 比值取负对数 end function p = phi(x, m, r) N = length(x); count = 0; total = 0; for i = 1 : N - m for j = 1 : N - m if i == j, continue; end d = max(abs(x(i:i+m-1) - x(j:j+m-1))); % Chebyshev距离 if d <= r count = count + 1; end total = total + 1; end end p = count / total; end

注意样本熵对数据长度有下限要求，官方指南建议RR间期序列不少于200个点，5分钟静息数据约300~400个点，刚好满足。r取0.15倍标准差时结果偏敏感，健康受试者SampEn通常在1.2~1.8之间；取0.25倍时区分度下降但对噪声更鲁棒。模糊熵与近似熵同样用m=2，模糊熵的梯度参数n取2，隶属度函数选择会使小差异被平滑，弱信号段的表现比样本熵稳定。

4.2 poincareHRV：SD1/SD2与椭圆拟合的几何含义

Poincare散点图把每个RR间期与前一个间期构成二维点，散点云被拟合为椭圆，SD1是垂直于恒等线的离散程度，SD2沿线方向展开，SD1/SD2比值反映短期与长期变异的关系。SD1与RMSSD在数学上高度相关（SD1约等于RMSSD除以根号2），但Poincare图的优势在于可以可视化整体散布形态——病理状态常表现为簇状或离散异常，纯数值无法捕捉。

% 计算Poincare几何指标 function feat = poincareHRV(rri) x = rri(1:end-1); y = rri(2:end); % 相邻间期点对 sd1 = std(x - y, 'omitnan') / sqrt(2); % 短期变异轴 sd2 = std(x + y, 'omitnan') / sqrt(2); % 长期变异轴 feat.SD1 = sd1; feat.SD2 = sd2; feat.SD1SD2 = sd1 / sd2; % 比值, 常作为恢复评估指标 % 椭圆面积, 总面积越大表示整体变异越强 feat.area = pi * sd1 * sd2; end

SD1/SD2比值正常范围在0.25~0.45之间，比值升高常见于短期变异增大，比如呼吸频率加快；比值降低且SD2同步缩小时往往对应长期调控能力下降。椭圆面积是另一个实用指标，运动疲劳状态下面积明显收窄，这个特征在运动员恢复监控里比LF/HF更稳定，因为它不依赖频带划分假设。如果数据里残留未替换的异常点，散点会形成离群尾巴，sd2会被显著拉伸，所以Poincare计算必须放在异常值替换之后。

4.3 非线性特征组合与exportHRV的输出结构

Extract_HRV_Nonlinear_Features把样本熵、近似熵、模糊熵和Poincare指标汇总成特征向量，exportHRV负责把结果写到结构化表格。导出格式上建议按行存样本、按列存特征，第一列是样本ID或时间段，后续每列一个特征，同时写入版本号和时间戳，便于回溯。

组合特征时要注意共线性问题：SD1与RMSSD相关系数经常在0.9以上，样本熵与近似熵之间也有强相关性，全量塞进分类模型会造成冗余。实践当中我先算特征相关矩阵，相关性超过0.85的只保留其中一个，通常保留RMSSD和样本熵，因为两者分别代表线性和非线性的互补信息。导出前最好对每个特征做Z-score归一化，避免样本熵的小数值被SDNN的上百毫秒数值压制。这套代码里PSD.xls和HRV-EEG.xls两个文件是示例输出，新数据跑完后对照它们的列格式检查字段名对齐即可。

5. main.m的调用次序与一个可复现的验收脚本

main.m的调用次序我建议固定为：读取数据、R波检测或IBI导入、异常值替换、线性特征、非线性特征、导出。特征计算之间的依赖关系是单向的，预处理必须先于所有特征，时域和频域可以并行算但结果要合并存储。直接改main.m的代价是每次跑实验都要翻动大量无关代码，更稳妥的方式是写一个独立脚本，只调用这套工具的函数。

% hrv_pipeline_verify.m - 验收管线正确性的最小脚本 clear; clc; fs = 500; % ECG采样率 [ecg, t] = read_ecg('sample_ecg.mat'); % 读取原始ECG rri = Extract_RRI(ecg, fs); % 1. R波定位与间期提取 rri = preProcessIBI(rri, fs); % 2. 重采样与格式统一 rri_clean = replaceOutliers(rri, 3); % 3. MAD异常值替换 td = timeDomainHRV(rri_clean); % 4. 时域特征 fd = freqDomainHRV(rri_clean, fs); % 5. 频域特征 nl = Extract_HRV_Nonlinear_Features(rri_clean); exportHRV(struct('time',td,'freq',fd,'nonlin',nl), 'out_features.xlsx'); fprintf('SDNN=%.2fms RMSSD=%.2fms SampEn=%.3f\n', ... td.SDNN, td.RMSSD, nl.SampEn);

验收时用一段已知质量的数据跑通后检查三点：SDNN与RMSSD数值落在合理区间，poincare的SD1近似等于RMSSD除以根号2，样本熵介于0.5到2.5之间。任何一步出现数量级异常，优先回头检查该函数传入的数据单位——RR间期到底是秒还是毫秒，这个错误在HRV特征提取里反复出现。代码在R2023b及以上版本直接运行即可，movmedian和movmad需要R2016a之后的版本，低版本环境可以手写滑动循环替代。最后把预处理参数记录在输出的Excel附注里，同一个数据集不同实验之间参数不一致，对比结论就没有意义了。

本文还有配套的精品资源，点击获取