简介基于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间期清理参数表参数取值作用下限300 ms过滤检测错误导致的过短间期上限2000 ms过滤长停搏或漏检造成的伪长间期movmedian窗口21点保持局部趋势占整段比例小ThresholdFactor55倍MAD减少生理性长间期被误删如果后续只算时域特征直接删除异常点更省事但要继续做频谱分析就建议保留长度避免时间戳出现空洞。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 msbin宽度决定熵的数值量级过窄会把噪声当成信息过宽则丢失分布细节10 ms是我在多组数据上试出来的稳定选择。第二概率为0的bin在熵公式里没有意义需要先剔除再求和否则log2(0)会直接产生NaN后面整条特征管线都会断掉。3.3 时域特征对照表与窗口长度选择不同时域特征对分析窗口长度的敏感度不同只用30秒数据算频域特征没有意义但算RMSSD已经足够。实际处理中固定长度窗口比变长窗口更便于不同记录之间比较特征值。下表是我常用的参数组合供搭建实验时直接参考特征计算方式房颤时典型变化主要干扰来源meanRR所有RR间期的均值通常缩短但个体差异大基础心率水平SDNN所有RR间期的标准差增大呼吸性窦性心律不齐RMSSD相邻差值平方后取均值的平方根显著增大房性早搏pNN50相邻差值超过50 ms的比例升高长间期集中片段Shannon熵RR间期直方图的信息熵增大bin宽度设置不当表里的窗口和阈值参数我在fs250 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)为止末尾不足一个插值步长的部分直接舍去不影响主频带估计。下面是我常用的谱估计参数参数取值说明插值采样率1 Hz覆盖0.5 Hz以下的心率频带窗函数Hann 512点频率分辨率约0.002 Hz重叠比例50%256点重叠趋势处理linear去除极低频趋势项窗长增大能压低频谱方差但也会消耗更多有效数据来做平均。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间期序列的模式重复性差样本熵显著升高。常用的参数组合是模板长度m2、相似容差r0.15*SDr跟随信号自身的标准差自动缩放不同振幅的数据之间具有可比性。function seVal sampEn(y, m, r) N numel(y); countB 0; countA 0; for i 1:N-m for j i1:N-m if max(abs(y(i:im-1) - y(j:jm-1))) r countB countB 1; if abs(y(im) - y(jm)) 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和样本熵产生可观测的变化。我在保存特征时会把关键参数一并写入structoutFeat.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回退和对比实验时也更好定位。本文还有配套的精品资源点击获取