尧图网页设计 
简介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能量。小波基典型场景分解层数使用要点db4ECG基线漂移去除5~8与QRS形态相似度低不易把R波能量滤掉sym8IBI序列的平滑与重采样6对称性好相位失真小适合保留间期细节coif5平稳段HRV趋势提取4~6对低频段的保真度优于db族小波去噪最常见的坑是把分解层数设得过大导致QRS波群被当成细节系数置零滤波后的ECG里R峰幅值缩水R波检测漏检率上升。另一个容易忽略的点是软阈值与硬阈值的差异硬阈值保留峰值但会产生震荡伪迹软阈值平滑但对R峰幅值有压缩。我这边的做法是在wavelet.m里对细节系数用软阈值、对近似系数不动只做基线校正这样后续R波检测的形态信息损失最小。参数调整后建议用滤波前后RR间期序列的相关性做校验相关系数低于0.95就该回查阈值设置。3. 线性HRV特征时域、频域与联合时频域的计算口径3.1 timeDomainHRVSDNN、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.04HzLF为0.04~0.15HzHF为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通常取2r取序列标准差的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); % m1维模板匹配对数 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:im-1) - x(j:jm-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倍时区分度下降但对噪声更鲁棒。模糊熵与近似熵同样用m2模糊熵的梯度参数n取2隶属度函数选择会使小差异被平滑弱信号段的表现比样本熵稳定。4.2 poincareHRVSD1/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; endSD1/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附注里同一个数据集不同实验之间参数不一致对比结论就没有意义了。本文还有配套的精品资源点击获取 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能量。小波基典型场景分解层数使用要点db4ECG基线漂移去除5~8与QRS形态相似度低不易把R波能量滤掉sym8IBI序列的平滑与重采样6对称性好相位失真小适合保留间期细节coif5平稳段HRV趋势提取4~6对低频段的保真度优于db族小波去噪最常见的坑是把分解层数设得过大导致QRS波群被当成细节系数置零滤波后的ECG里R峰幅值缩水R波检测漏检率上升。另一个容易忽略的点是软阈值与硬阈值的差异硬阈值保留峰值但会产生震荡伪迹软阈值平滑但对R峰幅值有压缩。我这边的做法是在wavelet.m里对细节系数用软阈值、对近似系数不动只做基线校正这样后续R波检测的形态信息损失最小。参数调整后建议用滤波前后RR间期序列的相关性做校验相关系数低于0.95就该回查阈值设置。3. 线性HRV特征时域、频域与联合时频域的计算口径3.1 timeDomainHRVSDNN、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.04HzLF为0.04~0.15HzHF为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通常取2r取序列标准差的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); % m1维模板匹配对数 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:im-1) - x(j:jm-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倍时区分度下降但对噪声更鲁棒。模糊熵与近似熵同样用m2模糊熵的梯度参数n取2隶属度函数选择会使小差异被平滑弱信号段的表现比样本熵稳定。4.2 poincareHRVSD1/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; endSD1/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附注里同一个数据集不同实验之间参数不一致对比结论就没有意义了。本文还有配套的精品资源点击获取