基于Pan-Tompkins算法的ECG信号HRV提取方案

一、算法原理与流程

1. 核心流程框架

2. 关键步骤详解

(1) 带通滤波设计

采用级联双滤波器消除基线漂移和工频干扰：

% 低通滤波器（截止11Hz）
b_lp = [0.003, 0.014, 0.023, 0.024, 0.014, 0.003];
a_lp = [1, -2.397, 2.846, -1.956, 0.540, -0.075];

% 高通滤波器（截止5Hz）
b_hp = [0.0884, -0.3536, 0.5355, -0.5355, 0.3536, -0.0884];
a_hp = [1, -2.397, 2.846, -1.956, 0.540, -0.075];

% 级联滤波
ecg_filt = filter(b_hp, a_hp, filter(b_lp, a_lp, ecg));

(2) 自适应阈值检测

动态更新信号峰值估计值(SP)和噪声峰值估计值(NP)：

SP = 0.125*max(ecg_filt) + 0.875*SP_prev;  % 信号峰值估计
NP = 0.125*max(noise_filt) + 0.875*NP_prev; % 噪声峰值估计
TH = 0.25*(SP - NP) + 0.5*NP;              % 动态阈值

(3) R波定位与回溯

采用Hamilton改进的回溯机制：

if isempty(R_peaks)
    % 回溯搜索窗口
    backwin = round(1.66*RR_avg);
    [peaks, locs] = findpeaks(ecg_filt(1:max(1,end-backwin)));
    % 阈值判定
    if max(peaks) > 0.5*TH
        R_peaks = locs(end);
    end
end

二、HRV特征提取实现

1. RR间期序列生成

% R波位置提取
[~, R_locs] = findpeaks(ecg_filt, 'MinPeakHeight', TH*0.8);

% RR间期计算
RR = diff(R_locs)/fs;  % 单位转换为秒
RR = RR(R_locs(2:end) - R_locs(1:end-1) > 200/fs);  % 去除异常间隔

2. 时域特征计算

% 时域指标
SDNN = std(RR)*1000;       % 总标准差(ms)
RMSSD = sqrt(mean(diff(RR).^2))*1000;  % 相邻差值均方根(ms)
pNN50 = sum(abs(diff(RR))>50/fs)/length(RR);  % 50ms变异比例

3. 频域特征计算

% 傅里叶变换
N = length(RR);
f = (0:N-1)*(fs/N);
Pxx = pwelch(RR, [], [], [], fs);

% 频段划分
LF_band = [0.04, 0.15];
HF_band = [0.15, 0.4];
LF = bandpower(Pxx, f, LF_band);
HF = bandpower(Pxx, f, HF_band);
LF_HF = LF/HF;

4. 非线性特征计算

% Poincaré分析
RR_mean = mean(RR);
x = RR - RR_mean;
SD1 = sqrt(0.5*var(x(2:end) + x(1:end-1)));
SD2 = sqrt(0.5*var(x(2:end) - x(1:end-1)));
EA = pi*SD1*SD2/2;  % 椭圆面积

三、MATLAB完整代码示例

function [hrv_features] = extract_HRV(ecg, fs)
    % 参数设置
    win_size = 256;  % 分析窗口长度(秒)
    overlap = 128;   % 窗口重叠

    % 1. Pan-Tompkins R波检测
    [R_locs, ~] = pan_tompkins(ecg, fs);

    % 2. RR间期提取
    RR = compute_rr(R_locs, fs);

    % 3. 时域特征
    SDNN = std(RR)*1000;
    RMSSD = sqrt(mean(diff(RR).^2))*1000;
    pNN50 = sum(abs(diff(RR))>50/fs)/length(RR);

    % 4. 频域特征
    [LF, HF, LF_HF] = compute_spectrum(RR, fs);

    % 5. 非线性特征
    [SD1, SD2, EA] = poincare_analysis(RR);

    % 输出结果
    hrv_features = struct(...
        'SDNN', SDNN, ...
        'RMSSD', RMSSD, ...
        'pNN50', pNN50, ...
        'LF', LF, ...
        'HF', HF, ...
        'LF_HF', LF_HF, ...
        'SD1', SD1, ...
        'SD2', SD2, ...
        'EA', EA);
end

function [R_locs] = pan_tompkins(ecg, fs)
    % 实现带通滤波、微分、平方、积分等步骤
    % ...（具体实现参考前文代码框架）
end

function RR = compute_rr(R_locs, fs)
    % 计算RR间期
    RR = diff(R_locs)/fs;
    RR = RR(R_locs(2:end) - R_locs(1:end-1) > 200/fs);  % 去除异常间隔
end

function [LF, HF, LF_HF] = compute_spectrum(RR, fs)
    % 功率谱密度计算
    N = length(RR);
    f = (0:N-1)*(fs/N);
    Pxx = pwelch(RR, [], [], [], fs);

    LF_band = [0.04, 0.15];
    HF_band = [0.15, 0.4];
    LF = bandpower(Pxx, f, LF_band);
    HF = bandpower(Pxx, f, HF_band);
    LF_HF = LF/HF;
end

function [SD1, SD2, EA] = poincare_analysis(RR)
    % Poincaré分析
    RR_mean = mean(RR);
    x = RR - RR_mean;
    SD1 = sqrt(0.5*var(x(2:end) + x(1:end-1)));
    SD2 = sqrt(0.5*var(x(2:end) - x(1:end-1)));
    EA = pi*SD1*SD2/2;
end

参考代码 利用pan_tompkin算法从ECG信号提取HRV算法 www.youwenfan.com/contentalh/98443.html

四、实验结果验证

数据集测试

使用MIT-BIH数据库进行验证：

五、工程应用建议

1. 硬件部署：在STM32H7系列MCU上实现轻量化部署（参考）
2. 临床验证：结合Holter设备进行24小时连续监测
3. 数据安全：采用AES加密传输生理数据（参考）