心电疲劳检测：从QRS波微变提取生理可解释特征

简介：本资源是一套基于MATLAB与机器学习算法实现的心电信号分析与驾驶/作业疲劳状态识别的完整开发方案，面向本科毕业设计、课程设计及工程实践初学者，解决生理信号建模、特征提取与分类判别等典型问题。压缩包共24个文件（4.17MB），含15个核心MATLAB函数（如Pan_Tompkins、LSTM_train、features_extract等）、3个交互式脚本（mlx）、2个Excel格式数据集（原始ECG与时域特征）、2个.mat模型文件（含训练好的net.mat与tm.mat）、1个README.md说明文档及1个注释文本，覆盖信号预处理、R波检测、特征工程、LSTM分类全流程。目前已有44人学习下载，提供开箱即用的可运行代码、结构清晰的模块化实现、配套数据集与详细注释，支持直接调试、参数调优与功能扩展，是理解生物医学信号智能分析落地的优质实践范例。

1. 心电疲劳检测不是“测心率”，而是从QRS波群微变里揪出大脑皮层供能不足的信号：这套 MATLAB+机器学习方案，专治毕业设计里“数据有、模型跑不通、结果像玄学”的三连翻车

你是不是也试过：用现成的心电信号（ECG）数据，套个 SVM 或者 LSTM，训练完准确率看着还行，但一换人、一换采集设备、甚至同一个人隔天再测，模型就突然“失忆”？不是代码写错了，是漏掉了心电疲劳检测最核心的物理前提——疲劳不是让心跳变快或变慢，而是让心电信号的时域形态发生亚毫秒级畸变：R波上升支斜率降低、T波不对称性增强、RR间期变异性（HRV）低频段能量衰减。这套资源不走“端到端深度学习”那种黑匣子路线，而是用 MATLAB 把 Pan-Tompkins 算法做稳、把 ARMA 模型提取的时域特征做实、再用 Birge-Massart 准则做特征筛选，最后喂给可解释性强的分类器。它适合真正要交源码、要讲清每一步物理意义、要答辩时被导师问“这个特征为什么能表征疲劳”的场景——比如某高校生物医学工程方向的课程设计，或者某实验室为模拟项目X搭建的基线系统。所有脚本都按模块切分，.mlx是交互式笔记，.m是可批量调用函数，.xlsx里存的是真实采集的多状态（清醒/轻度疲劳/重度疲劳）ECG 片段及其人工标注标签，不是合成噪声。

2. 从原始 ECG 到可用特征：Pan-Tompkins + 波形阈值分割 + 时域特征工程的三层过滤链

心电疲劳检测的成败，70% 取决于前处理是否鲁棒。这套资源没用现成的 WFDB 工具箱直接读.dat，而是把信号预处理拆成三个可验证、可调试、可替换的环节：先用 Pan-Tompkins 算法定位 R 波，再用wave_thresh_v1.m做自适应波形分割，最后用features_extract.m提取 18 维时域特征。这不是为了炫技，而是因为真实场景中，导联接触不良、工频干扰、基线漂移会让端到端模型直接崩溃。下面带你走通这条链路。

2.1 Pan-Tompkins 定位 R 波：为什么不用 findpeaks？因为真实 ECG 的 R 波不是“尖峰”，而是“宽峰”

Pan_Tompkins.m是整套流程的基石。它不依赖findpeaks这种通用峰值检测器，而是复现了经典论文里的五步滤波+微分+平方+移动窗积分流程。关键参数全暴露在函数头部：

function [r_peaks, r_amplitudes] = Pan_Tompkins(ecg_signal, fs) % fs: 采样率（Hz），必须精确！本数据集为 250 Hz % ecg_signal: 列向量，原始单导联 ECG（单位：mV） % 输出：r_peaks - R 波位置索引；r_amplitudes - 对应幅度（用于后续质量评估） % 【参数说明】 % 1. 带通滤波：[0.5, 40] Hz —— 滤掉基线漂移和高频肌电，保留 QRS 主频带 % 2. 微分算子：[-1, -2, 0, 2, 1] —— 增强 R 波上升沿，抑制 P/T 波 % 3. 平方运算：把负值归零，突出能量集中区 % 4. 移动窗积分：窗长 = round(0.15*fs) —— 0.15s 是典型 QRS 宽度，确保包络完整 % 5. 自适应阈值：初始阈值 = 0.7 * max(integrated_signal)，后续用双阈值动态更新

提示：r_amplitudes不是拿来分类的，而是用来筛掉低信噪比片段。我在某跨平台系统里加了一行逻辑：if r_amplitudes(i) < 0.3*mean(r_amplitudes), continue; end，直接剔除 12% 的伪迹片段，测试集准确率提升 4.2%。

2.2wave_thresh_v1.m：用局部标准差替代固定阈值，解决不同人 R 波幅度差异大的问题

Pan_Tompkins输出的是 R 波位置，但疲劳分析需要的是单个心跳周期的完整波形（P-Q-R-S-T）。wave_thresh_v1.m就是干这个的。它不简单地以 R 为中心截取 ±200ms，而是动态找 Q 和 T 的边界：

function [q_points, t_points] = wave_thresh_v1(ecg_signal, r_peaks, fs) % 输入：ecg_signal（同上）、r_peaks（来自 Pan_Tompkins）、fs % 输出：q_points（Q 波起始索引）、t_points（T 波结束索引） % 【核心逻辑】 % 1. 在每个 R 前 150ms 内，计算滑动窗（窗长=20ms）的标准差 % 2. Q 点定义为：标准差首次超过均值+1.5倍标准差的位置（即噪声突增点） % 3. 在每个 R 后 300ms 内，找 T 波终点：信号回归基线后，连续 50ms 的标准差 < 0.05 mV % 4. 若某次分割失败（如 T 波平缓），则回退到固定偏移：t_point = r_peak + round(0.35*fs) % 【为什么 v1 而不是 v0？】 % v0 用固定阈值 0.15 mV 找 Q，但在某同学采集的瘦弱受试者数据上，Q 波幅度仅 0.08 mV，全军覆没； % v1 改用局部统计量，对个体差异鲁棒得多。

这段代码的妙处在于：它把“找 Q/T 边界”这个主观任务，转化成了可量化的统计判据。你不需要调参，只需要确认fs是否正确——而这个值，wfdb_install.m会自动从.hea文件里读取（如果用 WFDB 数据），或由你手动填入（如心电采集数据.xlsx是 Excel 导出，需查原始设备设置）。

2.3features_extract.m：18 维时域特征不是随便堆的，每一维都对应一个生理可解释指标

features_extract.m是疲劳检测的“特征心脏”。它输入单个心跳波形（Q-T 区间），输出 18 个数值。这些不是 PCA 降维后的黑盒特征，而是明确定义的时域指标：

function features = features_extract(qrs_wave, rr_intervals, fs) % qrs_wave: 单个心跳波形（Q 到 T），列向量 % rr_intervals: 前 10 个 RR 间期（ms），行向量，用于计算 HRV 特征 % fs: 采样率 % 【关键步骤注释】 % 步骤1：标准化波形（消除个体基线差异） qrs_norm = (qrs_wave - mean(qrs_wave)) / std(qrs_wave); % 步骤2：计算 ARMA 残差（调用 ARMA_coeffs.m） [~, ~, resid] = ARMA_coeffs(qrs_norm, 2, 1); % p=2, q=1 features(18) = var(resid); % F18：残差方差，越大说明波形越“毛糙” % 步骤3：HRV 特征（只用最近 10 个 RR，避免长时记忆干扰） features(4) = std(rr_intervals); % F4：SDNN features(5) = mean(diff(rr_intervals).^2)^0.5; % F5：RMSSD，更敏感于迷走神经张力

注意：rr_intervals必须是同一受试者连续心跳的间隔。如果你的数据是分段采集的（如心电采集数据.xlsx里每行是独立 30s 片段），那么rr_intervals应取该片段内前 10 个 RR，不能跨片段拼接——否则 F4/F5 会严重失真。

3. 特征筛选与模型训练：Birge-Massart 准则砍掉冗余特征，LSTM 与 SVM 双路验证

有了 18 维特征，下一步不是直接扔进模型。真实项目里，特征维度高 ≠ 效果好，反而容易过拟合小样本。这套资源用Birge_Massart.m做特征筛选，再用LSTM_train.m和Copy_of_main.mlx分别训练时序模型与传统模型，形成交叉验证闭环。

3.1Birge_Massart.m：用统计显著性代替经验阈值，精准剔除“看起来有用、实际噪声”的特征

Birge_Massart.m实现的是基于最小描述长度（MDL）的特征选择。它不看相关系数，而是问：“保留这个特征，能让模型在测试集上的预测误差下降多少？这个下降值是否显著大于它引入的复杂度代价？” 具体到代码：

function [selected_features, scores] = Birge_Massart(X, y, alpha) % X: 特征矩阵（n_samples x n_features），已标准化 % y: 标签向量（n_samples x 1），[0,1] 表示非疲劳/疲劳 % alpha: 显著性水平，默认 0.05 % 【算法逻辑】 % 1. 对每个特征 f_i，计算其与 y 的互信息 I(f_i; y) % 2. 计算该特征的“描述长度代价”：DL_i = log2(n_samples) + 0.5*log2(var(f_i)) % 3. 计算“收益-代价比”：ratio_i = I(f_i; y) / DL_i % 4. 保留 ratio_i > alpha * max(ratio) 的特征（alpha=0.05 意味着只留 top 5%） % 【实测效果】 % 在本数据集上，18 维 → 选中 9 维：F1(R波幅)、F2(QRS宽)、F4(SDNN)、F5(RMSSD)、F7(T波不对称度)、 % F11(P波/R波比)、F13(QT间期)、F16(RR变异性熵)、F18(ARMA残差方差) % 剔除了 F6(F波幅度)、F9(ST段斜率)等在疲劳状态下变化不显著的维度

血泪经验：我第一次用全部 18 维训 LSTM，在 5 折交叉验证里，有一折的准确率只有 63%（其他四折 85%+）。用Birge_Massart筛一遍后，五折全部稳定在 82%~86%。波动变小，才是模型真正学到规律的标志。

3.2LSTM_train.m：不是堆层数，而是用“心跳序列”建模疲劳的累积效应

LSTM_train.m的输入不是单个心跳的 9 维特征，而是连续 N 个心跳的特征矩阵（N×9）。这抓住了疲劳的本质——它是生理状态的缓慢漂移，不是单次心跳的突变。代码关键部分：

function net = LSTM_train(X_seq, y_seq, options) % X_seq: 三维数组（N_timesteps x N_features x N_sequences） % y_seq: 标签（N_sequences x 1），每个序列一个标签（疲劳/非疲劳） % options: 结构体，含 'MaxEpochs', 'MiniBatchSize', 'InitialLearnRate' % 【网络结构精要】 layers = [ sequenceInputLayer(9, 'Normalization','zscore') % 输入9维特征，z-score标准化 lstmLayer(64, 'OutputMode','last') % 只取最后一个LSTM单元输出（代表整个序列的“状态总结”） dropoutLayer(0.3) % 防止过拟合，因数据量小（<500序列） fullyConnectedLayer(2) % 二分类 softmaxLayer classificationLayer]; % 【训练参数真相】 % 'MiniBatchSize' = 8：不是越大越好！本数据集序列总数仅 427，batch=8 保证每 epoch 覆盖足够多样本 % 'InitialLearnRate' = 0.01：太大易震荡，太小收敛慢；用 trainNetwork 自带的 'adam' 优化器 % 'Plots' = 'training-progress'：必须开着，观察 validation loss 是否早停

提示：X_seq的构造是关键。data_split.m里做了滑动窗口：取连续 15 个心跳（约 15s），标签取这 15s 内的主导状态（如 15 个心跳里 12 个标记为疲劳，则整个序列标为疲劳）。这比单心跳分类更能抵抗瞬时伪迹。

3.3Copy_of_main.mlx：用 SVM 做“快速验证锚点”，5 行代码跑通 baseline

Copy_of_main.mlx是交互式笔记本，专为课程设计答辩准备。它用fitcsvm训练一个线性 SVM，5 行核心代码搞定：

% mlx 笔记本中的可执行块 X_train = X_selected(train_idx,:); % 已筛选的9维特征 y_train = y(train_idx); X_test = X_selected(test_idx,:); y_test = y(test_idx); svmModel = fitcsvm(X_train, y_train, 'KernelFunction','linear', 'Standardize',true); y_pred = predict(svmModel, X_test); acc_svm = sum(y_pred == y_test) / length(y_test); % 通常 78%~82%

为什么必须有 SVM？
因为 LSTM 需要调参、需要 GPU、需要解释“为什么这个序列被判疲劳”。而 SVM 是你的可信锚点：当 LSTM 准确率是 85%，SVM 是 80%，你知道提升是真实的；如果 LSTM 是 95% 而 SVM 只有 60%，那大概率是数据泄露或标签错误。某同学在答辩时被问“你怎么证明不是过拟合？”，他当场打开Copy_of_main.mlx，把 SVM 和 LSTM 的混淆矩阵并排一放，评委立刻点头——这就是可解释性的力量。

4. 避坑：那些让心电疲劳检测从 85% 掉到 50% 的真实陷阱（附现象、原因、解法）

这套资源经过严格测试，但“严格测试”不等于“无坑”。我在帮某实验室部署模拟项目X时，踩过以下五个典型坑，每一个都曾让准确率腰斩。这里不讲理论，只说现象、原因、解法，照着做就能绕开。

4.1 现象：Pan_Tompkins.m在某段数据上完全找不到 R 波，r_peaks返回空数组

原因：该段 ECG 存在严重基线漂移（drift > 1 mV），而Pan_Tompkins的带通滤波下限是 0.5 Hz，无法抑制这种超低频漂移。
解法：在调用Pan_Tompkins前，加一级高通滤波：

% 加在 main.mlx 或你自己的主脚本里 ecg_filtered = highpass(ecg_signal, 0.3, fs, 'ImpulseResponse','iir'); % 截止频率0.3Hz [r_peaks, ~] = Pan_Tompkins(ecg_filtered, fs);

注意：highpass是 MATLAB R2018a+ 内置函数。若用旧版，用designfilt('highpassiir','FilterOrder',4,'HalfPowerFrequency',0.3,'SampleRate',fs)替代。

4.2 现象：wave_thresh_v1.m找到的 T 波终点总滞后 100ms，导致 QRS 波形被截断

原因：wave_thresh_v1.m默认假设 T 波在 R 后 300ms 内结束，但某些受试者（尤其老年人）QT 间期可达 450ms。
解法：修改wave_thresh_v1.m中的搜索上限：

% 找到原代码第 47 行左右： search_window_end = r_peak + round(0.3 * fs); % 原为 0.3s % 改为： search_window_end = r_peak + round(0.45 * fs); % 支持最长 QT 间期 450ms

提示：改完后，务必用plot_confusion.m查看 T 波分割效果——正常应看到 T 波终点紧贴波形回落基线处，而非悬在半空。

4.3 现象：features_extract.m计算的 F18（ARMA 残差方差）在所有样本上都是 NaN

原因：ARMA_coeffs.m内部用了arima函数拟合，而arima要求输入波形长度 ≥ 50 点。若wave_thresh_v1截取的波形太短（如 < 40 点），arima报错返回 NaN。
解法：在features_extract.m开头加长度检查：

if length(qrs_wave) < 50 features(18) = 0.1; % 设一个合理默认值，避免传播 NaN warning('QRS wave too short for ARMA fitting, using default F18=0.1'); return; end

4.4 现象：LSTM_train.m训练时 validation loss 持续上升，但 training loss 下降

原因：data_split.m默认随机打乱数据，但 ECG 数据有强烈时间相关性——把同一受试者的前后片段分到 train/test，造成数据泄露。
解法：改用“受试者分层划分”：

% 在 data_split.m 中，找到随机划分部分，替换为： subject_ids = [1,1,1,2,2,2,3,3,3,...]; % 你的数据中每个样本对应的受试者ID cv = cvpartition(subject_ids, 'HoldOut', 0.2); % 按受试者分层 train_idx = training(cv); test_idx = test(cv);

关键：subject_ids必须是你自己根据数据来源手动标注的，心电采集数据.xlsx里没有这一列，需你补充。

4.5 现象：plot_confusion.m画出的混淆矩阵里，“疲劳”类召回率（Recall）只有 40%，但“非疲劳”类高达 95%

原因：数据集严重不平衡——数据集(时域特征).xlsx中，非疲劳样本 320 个，疲劳样本仅 107 个。模型学会“永远预测非疲劳”来刷准确率。
解法：三管齐下：

1. 重采样：在data_split.m后加smote（需 Statistics and Machine Learning Toolbox）：

   [X_balanced, y_balanced] = smote(X_train, y_train, 'NumNeighbors',3);

2. 损失加权：LSTM_train.m中，trainingOptions加'ClassWeights','balanced'；
3. 评估改用 F1-score：classificationReport = classificationReport(y_test, y_pred)，盯住F1-Score列，不是 Accuracy。

5. 模型部署与实时检测：把.mat网络转成函数，用main.mlx一键加载新数据跑预测

毕业设计和课程设计的终点，不是“模型训出来”，而是“别人能用你的代码跑通”。这套资源的终极价值，在于net.mat和tm.mat这两个文件——它们是训练好的模型和特征标准化参数，让你跳过耗时的训练过程，直接进入预测环节。下面教你如何用main.mlx把新 ECG 数据喂进去，5 分钟得到疲劳判断。

5.1net.mat与tm.mat：不是黑盒，而是可追溯的部署资产

net.mat里存的是trainNetwork输出的SeriesNetwork对象，tm.mat里存的是fitcsvm或trainNetwork用到的Standardize参数（均值、标准差）。它们的关系是：

重要提醒：tm.mat的mu/sigma是对筛选后特征（9 维）计算的，不是对原始 18 维！所以你不能直接用X_new（18 维）去除tm.mat里的mu——必须先用Birge_Massart.m筛出相同 9 维，再标准化。

5.2main.mlx：三步走通预测流，支持 Excel/CSV/实时串口数据

main.mlx是交互式主入口，它把整个 pipeline 封装成三个清晰步骤。打开它，按顺序执行即可：

步骤 1：加载并预处理新数据

% 【单元格1】加载你的新ECG数据（支持Excel/CSV） new_data = readmatrix('my_new_ecg.csv'); % 或 readtable('my_new_ecg.xlsx') ecg_signal = new_data(:,1); % 假设第一列是ECG fs = 250; % 务必确认你的采样率！ % 【单元格2】运行全流程（自动调用Pan_Tompkins→wave_thresh_v1→features_extract→Birge_Massart） [r_peaks, features_matrix] = rdata(ecg_signal, fs); % rdata.m 是整合脚本，内部已硬编码调用所有必要函数

步骤 2：加载模型与标准化参数

% 【单元格3】 load('net.mat'); % 加载LSTM网络 load('tm.mat'); % 加载标准化参数 % 对 features_matrix 做标准化（注意：只取Birge-Massart选中的9维） selected_idx = [1,2,4,5,7,11,13,16,18]; % 这是本资源默认选中的9维索引 X_new = features_matrix(:, selected_idx); X_new_norm = (X_new - tm.mu) ./ tm.sigma;

步骤 3：预测与可视化

% 【单元格4】 % LSTM预测（输出概率） [~, scores] = classify(trainedNet, X_new_norm); prob_fatigue = scores(:,2); % 第二列是"疲劳"类概率 % 生成决策（阈值0.5） pred_label = scores(:,2) > 0.5; % 【单元格5】画图：原始ECG + R波标记 + 疲劳概率曲线 figure; subplot(2,1,1); plot(ecg_signal); hold on; plot(r_peaks, ecg_signal(r_peaks), 'ro'); title('ECG with R-peaks'); subplot(2,1,2); plot(prob_fatigue); ylabel('Fatigue Probability'); xlabel('Heartbeat Index'); yline(0.5, '--r', 'Decision Threshold');

实战技巧：main.mlx里有个隐藏开关——把rdata.m中的use_lstm = true改成false，它就会自动切换到 SVM 预测路径（用tm.mat里的svmModel）。这样你能在同一界面对比两种模型，答辩时随时切换演示。

5.3 实时检测扩展：用wfdb.mlx接硬件，把串口数据流喂进 pipeline

wfdb.mlx不是读.dat文件，而是演示如何接入真实硬件。它用serialport对象监听 COM 端口，每收到 1024 点数据就触发一次处理：

% wfdb.mlx 中的关键循环 s = serialport('COM3', 115200); % 根据你的设备改端口号和波特率 configureCallback(s, 'terminator', @onDataReceived); function onDataReceived(src, event) data = read(s, 1024, 'double'); % 读1024点浮点数 % 后续调用 rdata(data, fs) 流程... % 注意：此处 fs 必须与硬件设置严格一致！ end

血泪教训：某同学用 Arduino 采集 ECG，波特率设 9600，但read函数默认超时 1 秒，导致每次只读到 960 点，rdata报错。解决方案：read(s, 1024, 'double', 'Timeout', 0.05)，把超时设为 50ms，并在 Arduino 端用Serial.flush()清空缓冲区。从那以后我每次接串口，都强制走一遍read+timeout+flush三件套。

希望帮到你。

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