基于ECG信号的HRV时域与频域分析Matlab代码实现

一、核心流程

graph TD A[ECG数据导入] --> B[预处理] B --> C[R波检测] C --> D[RR间期提取] D --> E[时域分析] D --> F[频域分析] E --> G[可视化] F --> G

---

二、完整Matlab代码实现

1. 数据预处理（去噪与基线校正）

function clean_ecg = preprocess_ecg(ecg, fs)
    % 带通滤波（0.5-40Hz）去除工频干扰和基线漂移
    [b,a] = butter(3, [0.5 40]/(fs/2));
    filtered = filtfilt(b,a,ecg);
    
    % 小波去噪（db4小波，5层分解）
    [c,l] = wavedec(filtered,5,'db4');
    denoised = waverec(wthresh(c, 's', 3*std(c)/2), l, 'db4');
    
    % 移动平均平滑（窗口宽度=50ms）
    window = round(0.05*fs);
    clean_ecg = movmean(denoised, window);
end

2. R波检测（改进Pan-Tompkins算法）

function [r_peaks, locs] = detect_r_peaks(ecg, fs)
    N = length(ecg);
    window = round(0.15*fs);  % 积分窗口
    
    % 平方与积分
    squared = ecg.^2;
    integrated = filter(ones(1,window)/window, 1, squared);
    
    % 动态阈值
    threshold = 0.7*max(integrated(1:1000)) + 0.3*std(integrated(1:1000));
    
    % 峰值检测
    [peaks, locs] = findpeaks(ecg, 'MinPeakHeight', threshold, ...
        'MinPeakDistance', round(0.2*fs), 'Prominence', 0.5*threshold);
end

3. RR间期提取

function rr_intervals = get_rr_intervals(locs, fs)
    rr_samples = diff(locs);
    rr_intervals = rr_samples / fs * 1000;  % 转换为毫秒
    
    % 异常值处理（3σ原则）
    mean_rr = mean(rr_intervals);
    std_rr = std(rr_intervals);
    rr_intervals(rr_intervals < mean_rr-3*std_rr | rr_intervals > mean_rr+3*std_rr) = [];
end

4. 时域分析

function [sdnn, rmssd, pnn50] = time_domain(rr_intervals)
    N = length(rr_intervals);
    
    % 基本指标
    sdnn = std(rr_intervals);
    
    % 差分序列
    diffs = diff(rr_intervals);
    
    % RMSSD
    rmssd = sqrt(mean(diffs.^2));
    
    % pNN50
    nn50 = sum(diffs > 50);
    pnn50 = nn50 / N;
end

5. 频域分析（Welch方法）

function [lf_power, hf_power, lf_hf_ratio] = freq_domain(rr_intervals, fs)
    N = length(rr_intervals);
    window = hamming(256);
    noverlap = 128;
    
    % 功率谱估计
    [Pxx,f] = pwelch(rr_intervals, window, noverlap, [], fs);
    
    % 频段划分
    lf_band = [0.04, 0.15];
    hf_band = [0.15, 0.4];
    
    % 频段积分
    lf_idx = f >= lf_band(1) & f <= lf_band(2);
    hf_idx = f > hf_band(1) & f <= hf_band(2);
    
    lf_power = trapz(f(lf_idx), Pxx(lf_idx));
    hf_power = trapz(f(hf_idx), Pxx(hf_idx));
    
    % LF/HF比值
    lf_hf_ratio = lf_power / hf_power;
end

6. 可视化模块

function plot_hrv_results(ecg, rr_intervals, locs, sdnn, rmssd, lf_hf_ratio)
    figure;
    
    % ECG波形与R峰标记
    subplot(3,1,1);
    plot(ecg);
    hold on;
    plot(locs, ecg(locs), 'ro');
    title('ECG波形与R峰检测');
    xlabel('时间 (s)');
    ylabel('幅值 (mV)');
    
    % RR间期序列
    subplot(3,1,2);
    stem(rr_intervals, 'b', 'LineWidth', 1.5);
    title('RR间期序列');
    xlabel('心跳序号');
    ylabel('间隔 (ms)');
    
    % 时频域指标
    subplot(3,1,3);
    bar([sdnn, rmssd, lf_hf_ratio]);
    set(gca, 'XTickLabel', {'SDNN (ms)', 'RMSSD (ms)', 'LF/HF'});
    ylabel('数值');
    title('HRV时频域指标');
end

---

三、完整调用示例

%% 参数设置
fs = 250; % 采样率 (Hz)
ecg = load('ecg_signal.mat'); % 加载ECG数据（列向量）

%% 预处理
clean_ecg = preprocess_ecg(ecg, fs);

%% R波检测
[r_peaks, locs] = detect_r_peaks(clean_ecg, fs);

%% RR间期提取
rr_intervals = get_rr_intervals(locs, fs);

%% 时域分析
[sdnn, rmssd, pnn50] = time_domain(rr_intervals);

%% 频域分析
[lf_power, hf_power, lf_hf_ratio] = freq_domain(rr_intervals, fs);

%% 可视化
plot_hrv_results(clean_ecg, rr_intervals, locs, sdnn, rmssd, lf_hf_ratio);

---

四、结果示例

指标	正常范围	临床意义
SDNN	50-100 ms	总体自主神经活性
RMSSD	20-80 ms	副交感神经张力
LF/HF	0.5-2.0	交感-副交感平衡

参考代码 采用ECG信号进行HRV的时域分析、频域分析的Matlab代码 www.youwenfan.com/contentcnq/55064.html

五、扩展功能

1. 呼吸耦合分析

   % 呼吸信号同步处理
   [resp_rate, phase] = respiratory_analysis(ecg, fs);
   

2. 非线性分析

   % Poincaré图
   [sd1, sd2] = poincare_plot(rr_intervals);
   

3. 机器学习集成

   % 使用HRV特征训练分类模型
   features = [sdnn, rmssd, lf_hf_ratio];
   model = train_classifier(features, labels);
   

---

六、注意事项

1. 数据质量验证：需检查QRS波群形态（使用ECG诊断工具箱）

2. 采样率要求：建议≥250Hz（MIT-BIH标准）

3. 临床验证：需与专业医疗设备（如Holter）对比

---

参考文献：

ECG信号预处理与R波检测算法优化

HRV时域频域特征计算方法

基于Welch方法的功率谱估计

运动伪影消除技术