从心电信号提取呼吸波形并计算HRV频谱参数
从心电信号中提取呼吸波形并计算HRV频谱参数是一个完整的生理信号处理流程。
一、完整处理流程概述
- 输入：原始心电信号（ECG）
- 预处理：滤波、去噪、R波检测
- 呼吸波形提取：从ECG中提取呼吸信号（EDR）
- 呼吸基频计算：从呼吸波形中提取呼吸频率
- HRV分析：从R-R间期计算心率变异性
- 频谱分析：计算HRV的功率谱密度
- 参数计算：提取LF、HF功率及衍生参数
二、MATLAB完整实现代码
%% 主函数：从ECG信号提取呼吸波形并计算HRV频谱参数
function [resp_signal, resp_rate, HRV_params] = extract_resp_from_ecg(ecg_signal, fs, varargin)
% 从心电信号提取呼吸波形并计算HRV频谱参数
% 输入：
% ecg_signal - 原始ECG信号（向量）
% fs - 采样频率（Hz）
% varargin - 可选参数：'plot_flag'（是否绘图）
% 输出：
% resp_signal - 提取的呼吸波形
% resp_rate - 呼吸基频（Hz）
% HRV_params - HRV频谱参数结构体
%% 1. 参数设置与初始化
if nargin > 2
plot_flag = varargin{1};
else
plot_flag = true;
end
% 频带定义（根据Task Force标准）
LF_band = [0.04, 0.15]; % 低频带 (Hz)
HF_band = [0.15, 0.4]; % 高频带 (Hz)
VLF_band = [0.003, 0.04]; % 超低频带 (Hz)
% 呼吸信号提取参数
resp_band = [0.1, 0.5]; % 正常呼吸频率范围 (Hz)
%% 2. ECG信号预处理
disp('步骤1: ECG信号预处理...');
ecg_clean = preprocess_ecg(ecg_signal, fs);
%% 3. R波检测与R-R间期计算
disp('步骤2: R波检测...');
[r_peaks, rr_intervals] = detect_r_peaks(ecg_clean, fs);
% 检查R波检测质量
if length(r_peaks) < 10
error('检测到的R波数量不足，请检查ECG信号质量');
end
%% 4. 从ECG提取呼吸波形（EDR）
disp('步骤3: 提取呼吸波形...');
resp_signal = extract_edr_signal(ecg_clean, r_peaks, fs);
%% 5. 计算呼吸基频
disp('步骤4: 计算呼吸基频...');
resp_rate = compute_respiratory_rate(resp_signal, fs, resp_band);
%% 6. HRV分析：R-R间期插值与频谱计算
disp('步骤5: HRV频谱分析...');
[lf_power, hf_power, lfnu, hfnu, lf_hf_ratio, total_power] = ...
compute_hrv_spectrum(rr_intervals, fs, LF_band, HF_band, VLF_band);
%% 7. 存储结果
HRV_params = struct();
HRV_params.LF_power = lf_power; % 低频绝对功率 (ms²)
HRV_params.HF_power = hf_power; % 高频绝对功率 (ms²)
HRV_params.LF_nu = lfnu; % 低频归一化功率 (nu)
HRV_params.HF_nu = hfnu; % 高频归一化功率 (nu)
HRV_params.LF_HF_ratio = lf_hf_ratio; % LF/HF比值
HRV_params.total_power = total_power; % 总功率 (ms²)
HRV_params.resp_rate_hz = resp_rate; % 呼吸频率 (Hz)
HRV_params.resp_rate_bpm = resp_rate * 60; % 呼吸频率 (次/分钟)
%% 8. 可视化结果
if plot_flag
visualize_results(ecg_clean, r_peaks, resp_signal, rr_intervals, ...
HRV_params, fs, LF_band, HF_band);
end
disp('处理完成！');
fprintf('呼吸频率: %.2f Hz (%.1f 次/分钟)\n', resp_rate, resp_rate*60);
fprintf('LF功率: %.2f ms², HF功率: %.2f ms²\n', lf_power, hf_power);
fprintf('LF/HF比值: %.2f\n', lf_hf_ratio);
end
%% 子函数1: ECG信号预处理
function ecg_clean = preprocess_ecg(ecg_signal, fs)
% ECG信号预处理：去基线漂移、滤波、去噪
% 1. 去基线漂移（高通滤波，截止频率0.5Hz）
[b_hp, a_hp] = butter(2, 0.5/(fs/2), 'high');
ecg_hp = filtfilt(b_hp, a_hp, ecg_signal);
% 2. 带通滤波（0.5-40Hz，保留QRS波信息）
[b_bp, a_bp] = butter(4, [0.5, 40]/(fs/2), 'bandpass');
ecg_bp = filtfilt(b_bp, a_bp, ecg_hp);
% 3. 工频干扰去除（50Hz陷波滤波）
wo = 50/(fs/2); % 归一化频率
bw = wo/35; % 带宽
[b_notch, a_notch] = iirnotch(wo, bw);
ecg_clean = filtfilt(b_notch, a_notch, ecg_bp);
% 4. 平滑处理（可选）
ecg_clean = smooth(ecg_clean, 5);
end
%% 子函数2: R波检测
function [r_peaks, rr_intervals] = detect_r_peaks(ecg_signal, fs)
% 使用Pan-Tompkins算法检测R波
% 1. 微分（突出QRS斜率）
diff_signal = diff(ecg_signal);
diff_signal = [diff_signal(1); diff_signal]; % 保持长度一致
% 2. 平方（增强高频成分）
squared_signal = diff_signal .^ 2;
% 3. 移动平均滤波（平滑）
window_size = round(0.15 * fs); % 150ms窗口
ma_filter = ones(window_size, 1) / window_size;
integrated_signal = conv(squared_signal, ma_filter, 'same');
% 4. 自适应阈值检测
% 寻找局部最大值
[~, locs] = findpeaks(integrated_signal, 'MinPeakHeight', ...
mean(integrated_signal) + 0.5*std(integrated_signal), ...
'MinPeakDistance', round(0.3*fs)); % 最小间隔300ms
% 5. 在原始ECG信号中精确定位R波
r_peaks = zeros(length(locs), 1);
search_window = round(0.05 * fs); % 在检测位置前后50ms搜索
for i = 1:length(locs)
start_idx = max(1, locs(i) - search_window);
end_idx = min(length(ecg_signal), locs(i) + search_window);
[~, max_idx] = max(ecg_signal(start_idx:end_idx));
r_peaks(i) = start_idx + max_idx - 1;
end
% 6. 计算R-R间期（单位：秒）
rr_intervals = diff(r_peaks) / fs;
% 7. 异常R-R间期检测与校正（基于生理范围）
min_rr = 0.3; % 最小R-R间期300ms（对应200bpm）
max_rr = 2.0; % 最大R-R间期2s（对应30bpm）
% 标记异常值
abnormal_idx = (rr_intervals < min_rr) | (rr_intervals > max_rr);
if any(abnormal_idx)
warning('检测到 %d 个异常R-R间期，将进行校正', sum(abnormal_idx));
% 使用中值滤波校正
rr_intervals(abnormal_idx) = median(rr_intervals(~abnormal_idx));
end
end
%% 子函数3: 从ECG提取呼吸波形（EDR）
function resp_signal = extract_edr_signal(ecg_signal, r_peaks, fs)
% 从ECG信号提取呼吸波形（ECG-derived respiration, EDR）
% 方法1: R波振幅调制法（最常用）
% 方法2: QRS面积法
% 方法3: 心电轴旋转法
% 这里使用R波振幅调制法
r_amplitudes = ecg_signal(r_peaks);
% 创建时间向量（R波时间点）
r_times = r_peaks / fs;
% 插值到均匀时间序列
time_vector = (0:length(ecg_signal)-1)' / fs;
resp_signal_raw = interp1(r_times, r_amplitudes, time_vector, 'spline', 'extrap');
% 带通滤波提取呼吸频段（0.1-0.5 Hz，对应6-30次/分钟）
[b_resp, a_resp] = butter(4, [0.1, 0.5]/(fs/2), 'bandpass');
resp_signal = filtfilt(b_resp, a_resp, resp_signal_raw);
% 归一化
resp_signal = (resp_signal - mean(resp_signal)) / std(resp_signal);
end
%% 子函数4: 计算呼吸基频
function resp_rate = compute_respiratory_rate(resp_signal, fs, resp_band)
% 从呼吸波形计算呼吸基频
% 1. 计算功率谱密度
[pxx, f] = pwelch(resp_signal, [], [], [], fs);
% 2. 限制在呼吸频段内
f_idx = (f >= resp_band(1)) & (f <= resp_band(2));
f_resp = f(f_idx);
pxx_resp = pxx(f_idx);
% 3. 寻找主峰（最大功率对应的频率）
[~, max_idx] = max(pxx_resp);
resp_rate = f_resp(max_idx);
% 4. 验证：确保有明确的峰值
% 计算峰值显著性（峰值功率与平均功率之比）
peak_prominence = pxx_resp(max_idx) / mean(pxx_resp);
if peak_prominence < 2
warning('呼吸信号峰值不显著，可能信号质量不佳');
end
end
%% 子函数5: 计算HRV频谱参数
function [lf_power, hf_power, lfnu, hfnu, lf_hf_ratio, total_power] = ...
compute_hrv_spectrum(rr_intervals, fs, LF_band, HF_band, VLF_band)
% 从R-R间期计算HRV频谱参数
% 1. R-R间期插值（转换为均匀时间序列）
% 创建R波时间点
r_times = cumsum([0; rr_intervals]); % 累积时间
% 插值到4Hz采样率（根据Task Force推荐）
interp_fs = 4; % Hz
interp_time = (0:1/interp_fs:r_times(end))';
rr_interp = interp1(r_times, [rr_intervals; rr_intervals(end)], interp_time, 'spline');
% 2. 去趋势（消除超低频趋势）
rr_detrended = detrend(rr_interp);
% 3. 计算功率谱密度（使用Lomb-Scargle周期图处理非均匀采样）
% 对于R-R间期，更推荐使用Lomb-Scargle方法
% 这里使用Welch方法作为近似（已插值为均匀采样）
[pxx, f] = pwelch(rr_detrended, [], [], [], interp_fs);
% 4. 计算各频带功率（单位：ms²）
% 注意：R-R间期单位是秒，转换为毫秒
rr_ms = rr_detrended * 1000; % 转换为毫秒
pxx_ms2 = pxx * (1000^2); % 功率单位转换为ms²
% 计算频带功率（积分）
lf_idx = (f >= LF_band(1)) & (f <= LF_band(2));
hf_idx = (f >= HF_band(1)) & (f <= HF_band(2));
vlf_idx = (f >= VLF_band(1)) & (f <= VLF_band(2));
lf_power = sum(pxx_ms2(lf_idx)) * (f(2)-f(1)); % 积分求功率
hf_power = sum(pxx_ms2(hf_idx)) * (f(2)-f(1));
vlf_power = sum(pxx_ms2(vlf_idx)) * (f(2)-f(1));
total_power = sum(pxx_ms2(f >= 0.003 & f <= 0.4)) * (f(2)-f(1));
% 5. 计算归一化功率和比值
lfnu = lf_power / (total_power - vlf_power) * 100;
hfnu = hf_power / (total_power - vlf_power) * 100;
lf_hf_ratio = lf_power / hf_power;
end
%% 子函数6: 结果可视化
function visualize_results(ecg_clean, r_peaks, resp_signal, rr_intervals, ...
HRV_params, fs, LF_band, HF_band)
% 可视化所有结果
figure('Position', [100, 100, 1200, 800]);
% 子图1: 原始ECG与R波检测
subplot(3, 3, 1);
t_ecg = (0:length(ecg_clean)-1)/fs;
plot(t_ecg, ecg_clean, 'b', 'LineWidth', 1);
hold on;
plot(r_peaks/fs, ecg_clean(r_peaks), 'ro', 'MarkerSize', 8, 'LineWidth', 2);
xlabel('时间 (s)');
ylabel('幅度');
title('ECG信号与R波检测');
legend('ECG信号', 'R波位置', 'Location', 'best');
grid on;
% 子图2: 提取的呼吸波形
subplot(3, 3, 2);
t_resp = (0:length(resp_signal)-1)/fs;
plot(t_resp, resp_signal, 'g', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('幅度');
title(sprintf('提取的呼吸波形 (频率: %.2f Hz)', HRV_params.resp_rate_hz));
grid on;
% 子图3: R-R间期序列
subplot(3, 3, 3);
t_rr = cumsum([0; rr_intervals]);
plot(t_rr(1:end-1), rr_intervals*1000, 'b-o', 'LineWidth', 1.5, 'MarkerSize', 4);
xlabel('时间 (s)');
ylabel('R-R间期 (ms)');
title('R-R间期序列 (HRV)');
grid on;
% 子图4: 呼吸信号频谱
subplot(3, 3, 4);
[pxx_resp, f_resp] = pwelch(resp_signal, [], [], [], fs);
plot(f_resp, 10*log10(pxx_resp), 'b', 'LineWidth', 1.5);
xlabel('频率 (Hz)');
ylabel('功率谱密度 (dB/Hz)');
title('呼吸信号功率谱');
xlim([0, 1]);
grid on;
hold on;
% 标记呼吸频率
plot([HRV_params.resp_rate_hz, HRV_params.resp_rate_hz], ylim, 'r--', 'LineWidth', 1.5);
legend('呼吸谱', sprintf('呼吸频率: %.2f Hz', HRV_params.resp_rate_hz));
% 子图5: HRV频谱（周期图）
subplot(3, 3, 5);
% 重新计算HRV频谱用于绘图
interp_fs = 4;
r_times = cumsum([0; rr_intervals]);
interp_time = (0:1/interp_fs:r_times(end))';
rr_interp = interp1(r_times, [rr_intervals; rr_intervals(end)], interp_time, 'spline');
rr_detrended = detrend(rr_interp);
[pxx_hrv, f_hrv] = pwelch(rr_detrended*1000, [], [], [], interp_fs);
plot(f_hrv, pxx_hrv, 'b', 'LineWidth', 1.5);
xlabel('频率 (Hz)');
ylabel('功率谱密度 (ms²/Hz)');
title('HRV功率谱密度');
xlim([0, 0.5]);
grid on;
hold on;
% 标记LF和HF频带
fill([LF_band(1), LF_band(2), LF_band(2), LF_band(1)], ...
[0, 0, max(pxx_hrv)*1.1, max(pxx_hrv)*1.1], 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
fill([HF_band(1), HF_band(2), HF_band(2), HF_band(1)], ...
[0, 0, max(pxx_hrv)*1.1, max(pxx_hrv)*1.1], 'g', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
legend('HRV谱', 'LF频带', 'HF频带', 'Location', 'best');
% 子图6: LF/HF功率比
subplot(3, 3, 6);
bar_values = [HRV_params.LF_power, HRV_params.HF_power, HRV_params.LF_HF_ratio];
bar_categories = {'LF功率', 'HF功率', 'LF/HF比值'};
bar(bar_values(1:2));
hold on;
plot(3, bar_values(3), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
set(gca, 'XTickLabel', bar_categories(1:2));
ylabel('功率 (ms²) / 比值');
title('HRV频谱参数');
text(3, bar_values(3), sprintf('LF/HF=%.2f', bar_values(3)), ...
'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center');
grid on;
% 子图7: Poincaré散点图（非线性HRV分析）
subplot(3, 3, 7);
rr_n = rr_intervals(1:end-1) * 1000; % 转换为ms
rr_n1 = rr_intervals(2:end) * 1000;
plot(rr_n, rr_n1, 'b.', 'MarkerSize', 10);
hold on;
plot([min(rr_n), max(rr_n)], [min(rr_n), max(rr_n)], 'r--', 'LineWidth', 1.5);
xlabel('RR_n (ms)');
ylabel('RR_{n+1} (ms)');
title('Poincaré散点图');
axis equal;
grid on;
% 计算SD1和SD2
rr_diff = rr_n1 - rr_n;
SD1 = std(rr_diff) / sqrt(2);
SD2 = std(rr_n + rr_n1) / sqrt(2);
text(0.7, 0.9, sprintf('SD1=%.1f ms\nSD2=%.1f ms', SD1, SD2), ...
'Units', 'normalized', 'BackgroundColor', 'white');
% 子图8: 参数汇总表
subplot(3, 3, 8);
axis off;
param_text = {
sprintf('呼吸频率: %.2f Hz (%.1f bpm)', ...
HRV_params.resp_rate_hz, HRV_params.resp_rate_bpm);
sprintf('LF功率: %.2f ms²', HRV_params.LF_power);
sprintf('HF功率: %.2f ms²', HRV_params.HF_power);
sprintf('LFnu: %.1f nu', HRV_params.LF_nu);
sprintf('HFnu: %.1f nu', HRV_params.HF_nu);
sprintf('LF/HF比值: %.2f', HRV_params.LF_HF_ratio);
sprintf('总功率: %.2f ms²', HRV_params.total_power);
};
text(0.1, 0.9, param_text, 'FontSize', 10, 'VerticalAlignment', 'top');
title('参数汇总');
% 子图9: 时频分析（可选）
subplot(3, 3, 9);
% 使用短时傅里叶变换查看HRV时频特性
if length(rr_intervals) > 50
[s, f_stft, t_stft] = spectrogram(rr_intervals*1000, 32, 30, 64, 1/mean(rr_intervals));
imagesc(t_stft, f_stft, 10*log10(abs(s)));
axis xy;
xlabel('时间 (s)');
ylabel('频率 (Hz)');
title('HRV时频分析');
colorbar;
caxis([-20, 20]);
else
text(0.5, 0.5, '数据不足进行时频分析', ...
'HorizontalAlignment', 'center', 'FontSize', 12);
axis off;
end
sgtitle('心电信号呼吸提取与HRV分析结果', 'FontSize', 14, 'FontWeight', 'bold');
end
三、使用示例
%% 示例1: 使用模拟ECG数据测试
% 生成模拟ECG信号（含呼吸调制）
fs = 1000; % 采样率1000Hz
t = 0:1/fs:60; % 60秒数据
% 生成ECG信号（简化模拟）
ecg_sim = generate_simulated_ecg(t, fs);
% 调用主函数
[resp_signal, resp_rate, HRV_params] = extract_resp_from_ecg(ecg_sim, fs, true);
%% 示例2: 加载真实ECG数据
% 假设有ECG数据文件（如MIT-BIH格式）
% [signal, fs] = load_ecg_data('ecg_record.dat');
% 调用主函数
% [resp_signal, resp_rate, HRV_params] = extract_resp_from_ecg(signal, fs, true);
%% 辅助函数：生成模拟ECG信号
function ecg_signal = generate_simulated_ecg(t, fs)
% 生成模拟ECG信号（含呼吸调制）
% 基础ECG波形（使用ECG合成模型）
heart_rate = 70; % 心率70bpm
resp_rate = 0.25; % 呼吸频率0.25Hz (15次/分钟)
resp_amplitude = 0.1; % 呼吸调制幅度
% 生成ECG基本波形
ecg_base = zeros(size(t));
for i = 1:length(t)
% 简化ECG模型：P波、QRS波、T波
phase = mod(t(i) * heart_rate/60, 1) * 2*pi;
% P波
if phase < pi/8
ecg_base(i) = 0.1 * sin(8*phase);
% QRS波
elseif phase < pi/4
ecg_base(i) = sin(16*phase);
% T波
elseif phase < pi/2
ecg_base(i) = 0.3 * sin(4*phase);
end
end
% 添加呼吸调制（R波振幅变化）
r_peak_times = 0:60/heart_rate:max(t);
r_amplitudes = 1 + resp_amplitude * sin(2*pi*resp_rate*r_peak_times);
% 创建调制后的ECG
ecg_signal = ecg_base;
for i = 1:length(r_peak_times)-1
idx = find(t >= r_peak_times(i) & t < r_peak_times(i+1));
if ~isempty(idx)
ecg_signal(idx) = ecg_signal(idx) * r_amplitudes(i);
end
end
% 添加噪声
noise_level = 0.05;
ecg_signal = ecg_signal + noise_level * randn(size(ecg_signal));
end
参考代码 从心电波形提取出呼吸波形，计算呼吸基频，通过HRV曲线来得到该频谱，最后计算LFa，RFa参数值 www.youwenfan.com/contentcst/63232.html
四、关键参数解释
1. 呼吸基频 (Respiratory Rate)
- 正常范围：0.1-0.5 Hz（6-30次/分钟）
- 临床意义：反映呼吸系统的功能状态
2. HRV频谱参数
-
LF (低频功率, 0.04-0.15 Hz)
- 反映交感神经和副交感神经的共同调节
- 与压力、情绪状态相关
- 单位：ms²
-
HF (高频功率, 0.15-0.4 Hz)
- 主要反映副交感神经（迷走神经）活动
- 与呼吸性窦性心律不齐相关
- 单位：ms²
-
LFnu, HFnu (归一化功率)
- 消除总功率变化的影响
- 计算公式：LFnu = LF/(总功率-VLF)×100
- 更稳定地反映自主神经平衡
-
LF/HF比值
- 反映交感-迷走神经平衡
- 正常范围：1.5-2.0
- 比值升高：交感神经活动增强
- 比值降低：副交感神经活动增强
五、注意事项
-
信号质量要求
- ECG信号采样率建议≥250Hz
- 信号应无明显运动伪影
- 记录时间建议≥5分钟以获得稳定频谱
-
算法局限性
- R波检测准确性直接影响所有后续分析
- 呼吸提取在心律失常患者中可能不准确
- 运动伪影会严重影响结果
-
临床应用
- 可用于睡眠呼吸监测
- 评估自主神经功能
- 监测压力和心理状态
- 评估呼吸-心脏耦合
-
验证方法
- 与呼吸带或胸腹运动传感器对比
- 在不同呼吸模式下验证
- 评估算法的鲁棒性
六、扩展功能
%% 扩展功能1: 批量处理多个ECG文件
function batch_process_ecg_files(file_list, output_dir)
% 批量处理ECG文件
for i = 1:length(file_list)
% 加载数据
[ecg_signal, fs] = load_ecg_file(file_list{i});
% 处理
[resp_signal, resp_rate, HRV_params] = extract_resp_from_ecg(ecg_signal, fs, false);
% 保存结果
save(fullfile(output_dir, sprintf('result_%d.mat', i)), ...
'resp_signal', 'resp_rate', 'HRV_params');
end
end
%% 扩展功能2: 实时监测（简化版）
function real_time_monitor(ecg_stream, fs, window_size)
% 实时监测呼吸和HRV参数
% ecg_stream: 实时ECG数据流
persistent buffer;
if isempty(buffer)
buffer = [];
end
% 添加新数据到缓冲区
buffer = [buffer; ecg_stream];
% 当缓冲区足够大时处理
if length(buffer) >= window_size * fs
% 提取窗口数据
window_data = buffer(end-window_size*fs+1:end);
% 处理
[resp_signal, resp_rate, HRV_params] = extract_resp_from_ecg(window_data, fs, false);
% 更新显示（此处简化）
fprintf('实时监测 - 呼吸: %.1f bpm, LF/HF: %.2f\n', ...
resp_rate*60, HRV_params.LF_HF_ratio);
% 清空旧数据（保持滑动窗口）
buffer = buffer(end-round(window_size*fs/2):end);
end
end
七、总结
本代码提供了从心电信号提取呼吸波形并计算HRV频谱参数的完整解决方案，包括：
- ECG预处理：滤波、去噪
- R波检测：Pan-Tompkins算法
- 呼吸提取：R波振幅调制法
- 频谱分析：功率谱密度计算
- 参数提取：LF、HF功率及衍生参数
- 可视化：多图展示处理结果
