Pan-Tompkins算法实现ECG信号R波检测

Pan-Tompkins算法实现ECG信号R波检测

算法概述
Pan-Tompkins算法是实时QRS波群检测的经典方法，通过信号处理链增强QRS特征并抑制噪声。核心步骤包括微分、平方、积分处理和自适应阈值检测。

实现步骤详解

信号预处理
输入ECG信号首先进行微分处理，突出QRS波群的斜率特征：

diff_ecg = [0;diff(ecgData)];  % 一阶差分保持长度一致

特征增强
平方运算放大高频分量，强化R波与其他波形的差异：

squared_ecg = diff_ecg.^2;

能量整合
移动平均窗口合并QRS区域能量，窗口宽度与心率范围相关：

window_width = round(0.15 * Fs);    % 150ms典型窗口
integrated_ecg = movmean(squared_ecg, window_width);

峰值检测
采用动态阈值策略，避免固定阈值对信号幅值变化的敏感性：

threshold = 0.05 * max(integrated_ecg);  % 自适应阈值比例
min_peak_dist = round(0.4 * Fs);         % 200BPM对应最小间隔

精确定位
在积分信号检测的候选位置附近搜索原始信号最大点，实现亚采样精度：

search_window = max(1,Peaks_locs(i)-round(0.05*Fs)) : ... 
                min(length(ecgData),Peaks_locs(i)+round(0.05*Fs));
[~, max_idx] = max(ecgData(search_window));

参数调优建议

• 窗口宽度：0.12-0.20秒对应正常心率范围
• 阈值比例：0.05-0.25根据信号质量调整
• 搜索半径：±50ms适应不同QRS波宽度

临床应用注意事项

1. 对运动伪影敏感，建议结合加速度计数据
2. 极低振幅QRS波需调整阈值策略
3. 房颤等心律失常情况需特殊处理

典型输出包含原始信号、处理过程信号和R峰标记位置，可用于后续RR间期分析和心率变异性研究。完整实现需考虑实时处理延迟、边界效应处理等工程细节。

function [qrsLocs, integrated_ecg, integrated_ecg_Peaks_locs] = panTompkins(ecg, Fs)
%% Pan-Tompkins算法提取ECG波形R波
% 输入：
%   ecgFiltered - 已滤波的 ECG 信号 (1D 数组)
%   Fs          - 采样率 (Hz)
%
% 输出：
%   qrsPeaks - 检测到的 R 峰幅值
%   qrsLocs  - 检测到的 R 峰位置（索引） 
%   integrated_ecg - 过程积分结果
%   qrsLocs  - 检测到的 R 峰位置（索引）

    ecgData = ecg;
    %% 微分处理（增强QRS斜率信息）
    diff_ecg = [0;diff(ecgData)];  % 保持长度一致
    
    %% 平方处理（放大高频分量）
    squared_ecg = diff_ecg.^2;
    
    %% 移动窗口积分（合并QRS能量）
    window_width = round(0.15 * Fs);    % 150ms窗口
    integrated_ecg = movmean(squared_ecg, window_width);
    
    %% 峰值检测
    % 动态阈值设置
    % threshold = 0.22 * max(integrated_ecg);
    threshold = 0.05 * max(integrated_ecg);
    min_peak_dist = round(0.4 * Fs);    % 最小峰间距（对应200BPM）
    
    % 寻找积分信号的峰值
    [integrated_ecg_Peaks, integrated_ecg_Peaks_locs] = findpeaks(integrated_ecg,...
        'MinPeakHeight', threshold,...
        'MinPeakDistance', min_peak_dist);
    
    % 在原始信号中精确定位R波
    qrsLocs = zeros(size(integrated_ecg_Peaks_locs));
    for i = 1:length(integrated_ecg_Peaks_locs)
        % 在积分峰附近±50ms范围内搜索最大点
        search_window = max(1,integrated_ecg_Peaks_locs(i)-round(0.05*Fs)) : min(length(ecgData),integrated_ecg_Peaks_locs(i)+round(0.05*Fs));
        [~, max_idx] = max(ecgData(search_window));
        qrsLocs(i) = search_window(max_idx);
    end
end