房颤信号特征提取：从RR间期到样本熵的MATLAB实现全流程

简介：基于MATLAB的房颤信号特征提取项目，面向生物医学工程、信号处理专业的研究者与学生，围绕心电图（ECG）中房颤这一常见心律失常，系统实现从原始信号导入、滤波预处理、R波检测到特征参数提取与分类评估的完整链路。压缩包共含7个文件，其中4个M脚本承担算法实现与流程控制，2个CSV文件提供正常与房颤心电样本数据，1个MAT文件保存处理后的数据矩阵，整体大小约2.03MB，结构清晰便于对照学习。目前已有202人学习下载。项目中预处理环节采用带通滤波抑制肌电干扰与工频噪声，R波检测通过峰值定位划分心动周期，进而计算RR间期和心率变异性指标，提取房颤相关特征；必要时还可接入统计、时频或非线性分析，并辅以分类模型完成模式识别。代码注释详细，逻辑明确，既适合初学者理解房颤信号处理的每个步骤，也可作为工程基础，进一步扩展至其他心律失常分析或更复杂的生物医学信号处理任务。

1. 房颤信号特征提取：先抓住RR间期，再谈P波

房颤信号特征提取这个任务，落到工程上要做的事情很明确：把一段体表心电信号变成一组数字，让这组数字能稳定地区分房颤心律和正常窦性心律。很多人上手第一反应是去检测P波，但P波幅值小、易受基线漂移和肌电干扰，在实际数据里稳定性很差。而QRS波幅值大、规律清晰，由相邻R波位置差分得到的RR间期序列，恰好能呈现房颤最典型的临床特征——绝对不规则。这篇文章给出的路径是数据读入、滤波、R波检测、RR间期序列清理，再到时域统计、频域估计和非线性熵特征计算，每一步都给出可运行的MATLAB代码和参数依据。适合研究心电信号处理的在校学生，也适合在医疗AI团队里想给分类模型补充可解释特征的工程师。

2. 数据读入与预处理：R波位置准确，特征才有意义

2.1 两条特征提取路线的输入差异

房颤信号特征提取在工程上分成两条路线。第一条是HRV统计路线，输入是一维RR间期序列，后续计算SDNN、RMSSD、样本熵等指标，这条路线的核心假设是房颤的心室反应绝对不规则，所有特征都在刻画“不规则”的不同侧面；第二条是心房活动分析路线，输入是消除QRS-T后的心电图残余信号，目标是提取f波主导频率、f波幅度等直接反映心房电活动的特征。两条路线的预处理方式不同，但起点一致，都需要先获得准确的R波时间位置。

拿到数据的第一步是确定格式。公共数据库如PhysioNet的MIT-BIH房颤数据库，可以用WFDB Toolbox的rdsamp直接读取，也可以从网页导出CSV或MAT格式再load进MATLAB。如果数据来自自己采集的设备，需要额外确认采样率、导联位置和A/D转换位深，这三项直接决定后面带通滤波器边界和R波检测阈值该怎么设。

2.2 带通滤波：把0.5-40 Hz之外的东西先切掉

原始ECG不能直接用来检测R波。0.5 Hz以下主要是基线漂移，也就是呼吸和身体移动让整条信号线缓慢上下浮动；40 Hz以上主要是肌电噪声。二者都会干扰QRS波的形态判断。我习惯用designfilt生成一个4阶Butterworth带通滤波器，再用filtfilt做零相位滤波：

fs = 250; % 采样率，单位Hz ecgRaw = ...; % 读入的原始心电向量 bpFilt = designfilt('bandpassiir', ... 'FilterOrder', 4, ... 'HalfPowerFrequency1', 0.5, ... 'HalfPowerFrequency2', 40, ... 'SampleRate', fs, ... 'DesignMethod', 'butter'); ecgFilt = filtfilt(bpFilt, ecgRaw);

用filtfilt而不是filter是有原因的。filtfilt对信号做一次正向滤波后再反向滤波一次，相位响应为零，R波峰值在滤波前后不会发生时间偏移。如果换成filter，滤波器会给每个频率分量引入不同时间延迟，R波坐标系统性偏移几个采样点，后面算RR间期时每一个值都会带上恒定误差。

2.3 R波检测：用自适应阈值替代固定阈值

R波检测最常见的稳定方案是Pan-Tompkins算法，流程是带通滤波、差分、平方、滑动窗口积分、自适应阈值。MATLAB的findpeaks函数可以省去其中一部分手工步骤，但直接写一个固定MinPeakHeight在长时程数据里基本会失败，因为不同片段的信号幅度差异很大。我一般先把信号做局部分段归一化，再用局部极值作为阈值基准：

localMax = movmax(abs(ecgFilt), round(fs * 1.5)); thr = 0.5 * localMax; [pks, locs] = findpeaks(ecgFilt, ... 'MinPeakHeight', mean(ecgFilt) + 0.3 * thr, ... 'MinPeakDistance', round(0.2 * fs));

这段代码的要点在MinPeakDistance。0.2秒对应300次/分钟的心率上限，低于这个间距的峰值直接丢弃，保证T波不会被当成R波。MinPeakHeight不是固定值，而是跟随1.5秒滑动窗口内的局部幅度变化，这样即使前10分钟信号幅值很大、后10分钟明显变小，阈值也会自动调整。R波检测完成后需要目视抽查，把前20秒的信号和检测点画在一起快速扫一眼，能解决大部分漏检和误检问题。

注意：MinPeakDistance的单位跟随采样率，500 Hz数据要改为round(0.2 * 500) = 100个采样点。

2.4 用diff生成RR间期序列并清理异常值

得到R波位置后，相邻两个位置差就是RR间期。换算成毫秒之后，先用生理范围做一次硬过滤，再用中位数滑动窗口处理残留的跳变点：

rrMs = diff(locs) / fs * 1000; validMask = (rrMs > 300) & (rrMs < 2000); rrMs = rrMs(validMask); rrClean = filloutliers(rrMs, 'linear', 'movmedian', 21, ... 'ThresholdFactor', 5);

参数说明：300 ms到2000 ms对应30到200 bpm的心率范围，超出这个范围的间期要么是检测错误，要么是极端的病理性长间期；filloutliers使用的是邻域中位数加5倍中位数绝对偏差的判断规则，超过阈值的点用线性插值替换而不是直接删除，这样序列长度保持不变，后面的频谱分析不会出现时间轴断裂。以下是我常用的RR间期清理参数表：

如果后续只算时域特征，直接删除异常点更省事；但要继续做频谱分析就建议保留长度，避免时间戳出现空洞。

3. 时域与HRV统计特征：房颤最直接的可计算证据

3.1 SDNN、RMSSD、pNN50为什么适合房颤

房颤发生时，心房率可以达到每分钟350到600次，但因为房室结不应期的随机过滤，心室反应完全没有规律。这种绝对不规则在RR间期序列上的表现，第一是整体离散程度增大，第二是相邻间期差异增大，第三是长短间期交替频繁。SDNN捕捉第一个特征，计算全部RR间期的标准差；RMSSD捕捉第二个特征，计算相邻RR间期差值的均方根；pNN50捕捉第三个特征，计算相邻RR间期差值超过50 ms的比例。

这三个指标单独看都存在干扰项。窦性心律合并频发房性早搏时，RMSSD会升高；深度呼吸引起的窦性心律不齐会让SDNN变大；而房颤伴缓慢心室率的时间窗内，pNN50会因为长间期片段增多而下降。所以我在实际项目里从来不拿单个时域特征做规则判定，而是把这一组指标作为特征向量送入后续的分类器。

3.2 把时域特征和Shannon熵写成一个函数

Shannon熵对分布形状非常敏感，宽而扁平的RR间期分布会得到更高熵值，恰好对应房颤的随机性。这一节把上述统计量封装成一个独立函数，输入毫秒单位的RR间期向量，输出一个结构体：

function feat = afTimeFeatures(rrMs) rrDiff = diff(rrMs); feat.meanRR = mean(rrMs); feat.sdnn = std(rrMs); feat.rmssd = sqrt(mean(rrDiff .^ 2)); feat.pnn50 = 100 * sum(abs(rrDiff) > 50) / numel(rrDiff); [counts, ~] = histcounts(rrMs, 200:10:1600); prob = counts / sum(counts); prob(prob == 0) = []; feat.shannonEnt = -sum(prob .* log2(prob)); end

两点说明。第一，直方图的bin宽度是10 ms，范围是200到1600 ms；bin宽度决定熵的数值量级，过窄会把噪声当成信息，过宽则丢失分布细节，10 ms是我在多组数据上试出来的稳定选择。第二，概率为0的bin在熵公式里没有意义，需要先剔除再求和，否则log2(0)会直接产生NaN，后面整条特征管线都会断掉。

3.3 时域特征对照表与窗口长度选择

不同时域特征对分析窗口长度的敏感度不同，只用30秒数据算频域特征没有意义，但算RMSSD已经足够。实际处理中，固定长度窗口比变长窗口更便于不同记录之间比较特征值。下表是我常用的参数组合，供搭建实验时直接参考：

表里的窗口和阈值参数，我在fs=250 Hz数据上调过，换成500 Hz采样率时，movmax窗口要翻倍到3秒，MinPeakDistance也要从0.2秒换算成对应采样点数。另一个容易忽略的点是滑窗策略必须全流程一致。有人把前1分钟用30秒窗口、后1分钟用60秒窗口，最后合并到一个特征文件里，模型性能自然不稳定。跨记录对比时，固定窗口长度比纠结特征本身的绝对值更重要。

4. 频域与非线性特征：从功率谱到样本熵

4.1 对RR间期序列做重采样：这一步不能省

对RR间期序列做频谱分析前，必须先认识到RR间期在时间轴上不是等间距的。直接对序列下标做FFT，得到的是伪谱，峰值位置与实际频率完全对不上。标准做法是把RR间期按时间戳插值到均匀采样网格上，再做Welch法功率谱估计：

tt = cumsum(rrMs) / 1000; % 相对时间秒 tt = tt - tt(1); fsRR = 1; % 插值目标采样率1 Hz ttResamp = 0:1/fsRR:tt(end); rrResamp = interp1(tt, rrMs, ttResamp, 'linear'); [pxx, freq] = pwelch(detrend(rrResamp, 'linear'), ... hann(512), 256, 512, fsRR);

代码逻辑分三段：先把每个RR间期累加得到心拍发生的时间点，然后以1 Hz为目标采样率做线性插值，最后用512点汉宁窗做功率谱估计。1 Hz在RR间期分析中足够，心率波动的主要能量集中在0.4 Hz以下。detrend去掉线性趋势，是为了避免整段信号缓慢升高的趋势在极低频区造成虚假功率。插值时记得到tt(end)为止，末尾不足一个插值步长的部分直接舍去，不影响主频带估计。下面是我常用的谱估计参数：

窗长增大能压低频谱方差，但也会消耗更多有效数据来做平均。512点窗在五分钟片段里是比较好的平衡点，过短的窗会让频带估计毛刺非常多。

4.2 LF/HF比值在房颤场景要慎用

传统HRV分析把功率谱划分为LF频带和HF频带。窦性心律下，HF峰与呼吸节律同步，LF峰与压力反射和血管运动张力相关。但房颤时RR间期的波动主要来自房室结的随机传导，呼吸性窦性心律不齐机制失效，LF和HF的生理学意义不复存在，直接套用LF/HF比值得到的是一个缺乏解释力的数字。我在房颤项目中一般更多关注0.01到0.4 Hz全频段的积分功率，以及频谱形状的宽带化程度，而不是纠结LF/HF的绝对值。

如果一定要和窦性心律做同一套特征对比，可以把LF和HF的绝对功率都保留下来，让分类器自己去学习权重，不要预先做比值假设。房颤数据里HF带能量普遍被宽带随机成分稀释，LF/HF比值经常出现异常大的波动，这种波动不是神经调节信息，而是房室结随机传导的伪影。

4.3 心房主导频率：f波主频的实用提取流程

如果想把房颤检测从心室反应不规则深入到心房活动异常，需要提取f波主导频率。困难在于体表心电图上QRS-T波能量远大于f波，直接用频谱分析会被QRS波主导。最常见的处理是平均心拍消减：把所有心拍按R波对齐叠加求平均，得到QRS-T模板，再从原始信号逐拍减去该模板，剩下的残差以f波为主。核心代码如下：

win = round(fs * 0.3); seg = zeros(2 * win + 1, numel(locs)); usable = locs > win & locs < numel(ecgFilt) - win; for i = find(usable)' seg(:, i) = ecgFilt(locs(i) - win : locs(i) + win); end qrsTemplate = mean(seg(:, usable), 2); residual = ecgFilt; for i = find(usable)' idx = locs(i) - win : locs(i) + win; residual(idx) = residual(idx) - qrsTemplate; end

注意R波靠近信号首尾时窗口会越界，所以先用usable掩码剔除这些不完整心拍。模板窗口取0.6秒，足够覆盖QRS波和T波的大部分能量。消减完成后再对residual做时频分析，在4-10 Hz频带内找每个时间窗的峰值频率，该频率记为这一段的主导频率。这个特征对电极位置和呼吸干扰比较敏感，单导联数据算出来的主导频率在不同记录之间可能偏差0.5 Hz以上。

4.4 样本熵：参数设置与MATLAB实现

样本熵是房颤特征里最依赖的非线性指标之一。它衡量的是序列中出现新模式的概率，值越大说明数据越随机；房颤时RR间期序列的模式重复性差，样本熵显著升高。常用的参数组合是模板长度m=2、相似容差r=0.15*SD，r跟随信号自身的标准差自动缩放，不同振幅的数据之间具有可比性。

function seVal = sampEn(y, m, r) N = numel(y); countB = 0; countA = 0; for i = 1:N-m for j = i+1:N-m if max(abs(y(i:i+m-1) - y(j:j+m-1))) <= r countB = countB + 1; if abs(y(i+m) - y(j+m)) <= r countA = countA + 1; end end end end if countB == 0 seVal = NaN; else seVal = -log(countA / countB); end end

这段代码采用双重循环加Chebyshev距离判断，逻辑直观但时间复杂度为O(N^2)，RR间期数量超过几千个点时会明显变慢。批量实验时可以考虑改成基于分箱加速的近似算法，或者用MEX重写内层循环。countB为零意味着序列完全缺乏模板匹配，这种情况下熵值无定义，返回NaN即可，上游特征表里要对NaN做专门处理，避免训练时整行样本被丢弃。

5. 特征有效性验证与排错：先怀疑数据，再怀疑算法

5.1 两个五分钟内的自检手段

特征提取流程跑通之后，我不会直接算完整批特征，而是先做两个廉价的检查。第一个是RR间期分布直方图：窦性心律是单峰分布，房颤是宽而扁甚至多峰的分布；如果直方图出现一个明显次级峰，大概率是QRS漏检形成的两倍RR间期。第二个是R波标注抽查，把检测结果画在原始信号上：

plot((1:fs*20)/fs, ecgFilt(1:fs*20)); hold on; plot(locs(locs < fs*20)/fs, ... ecgFilt(locs(locs < fs*20)), 'r^');

红色三角标记的位置如果有规律地落在T波波峰上，说明MinPeakDistance或MinPeakProminence设置不当。这两个检查都不涉及复杂统计学方法，但能避免大批错误特征进入后续训练，是成本最低的排错手段。

5.2 每个单一特征先过一个AUC检验

在把特征交给分类器之前，我习惯先用带标签的训练集对每个特征单独做一次ROC分析。MATLAB里perfcurve可以直接返回AUC，数值越高说明该特征单独区分房颤与窦性心律的能力越强：

[~, ~, ~, auc] = perfcurve(trainLabels, featValues, 1);

假设trainLabels是0/1标签向量，featValues是当前特征的列向量。我的经验是AUC大于0.8的特征可以单独作为规则判定的候选；0.6到0.8的特征保留下来做多特征融合；低于0.55的特征先检查是不是计算错误，而不是直接丢弃。某个特征在当前数据集上AUC低，可能是窗口长度不合适，换一个滑窗策略往往能救回来。

5.3 特征文件必须携带元数据

最后一个排错要点：特征值的可重复性依赖整条预处理链路的每个环节。滤波器阶数从4改成2、R波阈值从0.3改成0.5、异常值替换窗口从21改成31，都会让SDNN和样本熵产生可观测的变化。我在保存特征时会把关键参数一并写入struct：

outFeat.meta.fs = fs; outFeat.meta.filterOrder = 4; outFeat.meta.filtRangeHz = [0.5 40]; outFeat.meta.windowLenSec = 300; outFeat.meta.outlierMethod = 'filloutliers movmedian 21'; outFeat.features = featTable; save('af_feat.mat', '-struct', 'outFeat');

save的时候把版本号和时间戳也拼进文件名，比如af_feat_20260607_v3.mat，回退和对比实验时也更好定位。

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