基于多模态信号处理与机器学习的心电信号分析及心率变异性评估（MATLAB）

开始
│
├─ 参数设置
│   ├─ 设置采样率、持续时间、基础心率
│   ├─ 配置心率变异性参数
│   └─ 设置噪声参数
│
├─ 心电信号合成
│   ├─ 生成RR间期序列
│   ├─ 构建心搏模板
│   ├─ 合成干净ECG信号
│   └─ 添加噪声干扰
│
├─ 信号预处理
│   ├─ 高通滤波去除基线漂移
│   ├─ 陷波滤波去除电源干扰
│   └─ 带通滤波增强QRS波
│
├─ QRS波检测
│   ├─ 差分运算增强高频成分
│   ├─ 平方运算强化QRS波
│   ├─ 移动积分平滑信号
│   ├─ 自适应阈值检测候选点
│   └─ 局部最大值精炼R峰位置
│
├─ 心率变异性分析
│   ├─ 计算时域指标
│   ├─ 计算频域指标
│   ├─ 计算非线性指标
│   └─ 房颤样不规则性筛查
│
├─ 结果可视化
│   ├─ 绘制原始和滤波后信号
│   ├─ 标记R峰位置
│   ├─ 显示QRS增强信号
│   ├─ 绘制RR间期序列
│   ├─ 展示功率谱密度
│   └─ 呈现庞加莱图
│
└─ 结果报告
    ├─ 输出数值指标
    └─ 提供房颤筛查结果

function ecg_hrv_pipeline()
% ECG + HRV 高级处理流程（生物电工程项目）
% ---------------------------------------------------------------
% 本脚本功能：
% 1) 合成包含心率变异性和常见噪声的逼真心电图信号
% 2) 去噪（基线和电源干扰），增强QRS波，检测R峰
% 3) 计算心率变异的时域、频域和非线性指标
% 4) 筛查房颤样不规则性（简单启发式方法）
% 5) 可视化每个处理步骤
%
% 运行方法：
%   - 保存为 ecg_hrv_pipeline.m
%   - 在MATLAB中按运行（F5）
%
% 如果信号处理工具箱可用，我们使用butter/filtfilt进行零相位IIR滤波
% 否则，我们回退到简单的卷积方法


%% -------------------- 参数设置 --------------------
fs      = 360;              % 采样率（赫兹）
dur_s   = 180;              % 持续时间（秒）
hr_bpm  = 70;               % 标称心率（每分钟心跳次数）
af_mode = false;            % 设置为true可生成不规则的房颤样心律


% 心率变异性（RR间期变异性）设置
sdnn_ms     = 60;           % 目标SDNN（毫秒），用于变异性真实性
lf_hz       = 0.1;          % 低频调制（约梅耶波）
hf_hz       = 0.25;         % 高频调制（呼吸相关）


% 噪声设置
add_baseline_wander = true; % 是否添加基线漂移
baseline_amp_mV     = 0.15; % 基线漂移幅度（毫伏）
baseline_hz         = 0.35; % 基线漂移频率（赫兹）


add_powerline = true;       % 是否添加电源干扰
powerline_freq = 50;        % 根据地区选择50或60赫兹
powerline_amp_mV = 0.05;    % 电源干扰幅度（毫伏）


rng(42); % 设置随机种子以确保结果可重复


%% -------------------- 合成ECG信号 --------------------
N = fs*dur_s; % 计算总样本数
t = (0:N-1)/fs; % 生成时间向量


% 生成具有LF/HF变异性或房颤样抖动的RR间期序列
[rr_s, r_times] = generate_rr_series(dur_s, hr_bpm, sdnn_ms, lf_hz, hf_hz, af_mode);


% 通过累加具有形态学特征的心搏（高斯小波）构建合成ECG
ecg_clean = synth_ecg_from_rr(t, r_times, fs);


% 添加逼真噪声
ecg_noisy = ecg_clean; % 初始化含噪ECG
if add_baseline_wander % 添加基线漂移
    ecg_noisy = ecg_noisy + baseline_amp_mV * sin(2*pi*baseline_hz*t);
end
if add_powerline % 添加电源干扰
    ecg_noisy = ecg_noisy + powerline_amp_mV * sin(2*pi*powerline_freq*t);
end
% 添加白色传感器噪声
ecg_noisy = ecg_noisy + 0.02*randn(size(t));


%% -------------------- 去噪处理 --------------------
% 1) 去除基线漂移（高通滤波～0.5 Hz）
ecg_hp = hp_filter(ecg_noisy, fs, 0.5);


% 2) 电源干扰陷波滤波（50/60 Hz）
ecg_notch = notch_filter(ecg_hp, fs, powerline_freq);


% 3) 带通滤波聚焦QRS波（经典5-15 Hz；我们稍宽一些）
ecg_bp = bp_filter(ecg_notch, fs, 5, 25);


%% -------------------- QRS波增强 + R峰检测 --------------------
% Pan-Tompkins风格：差分 -> 平方 -> 移动积分
diff_sig = [0 diff(ecg_bp)]; % 一阶差分
sq_sig   = diff_sig.^2; % 平方运算增强高频成分
win_ms   = 120; % 积分窗口约120毫秒
win_n    = max(1, round(fs*win_ms/1000)); % 计算窗口样本数
int_sig  = movmean(sq_sig, win_n); % 移动平均积分


% 带不应期的自适应阈值检测
[locs_R, thr] = detect_r_peaks(int_sig, fs);


% 在原始滤波ECG上精炼R峰位置：在检测到的索引周围搜索局部最大值
search_radius = round(0.05*fs); % 搜索半径（样本数）
locs_R = refine_r_peaks(ecg_bp, locs_R, search_radius);


%% -------------------- 心率变异性指标计算 --------------------
RR  = diff(locs_R)/fs; % RR间期（秒）
NN  = RR(~isoutlier(RR, 'median')); % 使用中位数法剔除异常值，获取正常窦性心律间期
if numel(NN) < 3 % 检查是否有足够的数据点
    warning('NN间期数量太少——请尝试更长的持续时间或调整参数');
end


% 时域指标
metrics.meanNN_ms = mean(NN)*1000; % 平均NN间期（毫秒）
metrics.sdnn_ms   = std(NN)*1000; % NN间期标准差（毫秒）
metrics.rmssd_ms  = sqrt(mean(diff(NN).^2))*1000; % 相邻NN间期差值的均方根（毫秒）
metrics.pnn50     = 100*mean(abs(diff(NN)) > 0.050); % 相邻NN间期差值大于50毫秒的百分比


% 频域指标：通过对插值后的心率信号进行Welch谱分析
[lfnu, hfnu, lf_hz_est, hf_hz_est, lf_power, hf_power, total_power] = hrv_freq(NN, locs_R(1:end-1)/fs);


metrics.lfnu = lfnu; metrics.hfnu = hfnu; % 归一化低频和高频功率
metrics.lf_power = lf_power; metrics.hf_power = hf_power; metrics.total_power = total_power; % 绝对功率
metrics.lf_center_hz = lf_hz_est; metrics.hf_center_hz = hf_hz_est; % 频带中心频率


% 非线性指标：庞加莱图（SD1/SD2）
[SD1, SD2] = poincare_SD1_SD2(NN);
metrics.SD1_ms = SD1*1000; metrics.SD2_ms = SD2*1000; % 庞加莱图指标


% 简单房颤样不规则性筛查（启发式方法）：
%   - RR间期的高变异系数，低自相关性，升高的pNN50，降低的SD1/SD2比率稳定性
cv_rr   = std(RR)/mean(RR); % RR间期变异系数
ac_rr   = autocorr_stat(RR, 1); % 一阶自相关系数
metrics.cv_rr = cv_rr; metrics.ac_lag1 = ac_rr; % 存储指标


% 房颤筛查条件判断
af_flag = (cv_rr > 0.15) && (metrics.pnn50 > 15) && (ac_rr < 0.35);
metrics.af_screen_positive = af_flag; % 房颤筛查结果


%% -------------------- 结果报告与可视化 --------------------
figure('Color','w','Name','ECG & HRV处理流程','Position',[80 80 1200 800]); % 创建图形窗口


% 子图1：原始和含噪ECG对比
subplot(4,1,1);
plot(t, ecg_noisy, 'Color', [0.2 0.2 0.2]); hold on;
plot(t, ecg_clean, 'r');
legend('含噪ECG','干净ECG（真实值）'); xlabel('时间（秒）'); ylabel('毫伏');
title('含噪声的合成ECG信号');


% 子图2：滤波后ECG和检测到的R峰
subplot(4,1,2);
plot(t, ecg_notch); hold on;
stem(locs_R/fs, ecg_notch(locs_R), 'g','filled','MarkerSize',3);
xlabel('时间（秒）'); ylabel('毫伏');
title('滤波后ECG + 检测到的R峰');


% 子图3：QRS增强信号和阈值
subplot(4,1,3);
plot(t, int_sig); hold on; yline(thr, '--r');
xlabel('时间（秒）'); ylabel('任意单位');
title('QRS增强（积分信号）+ 检测阈值');


% 子图4：RR间期序列（心率变异性）
subplot(4,1,4);
plot(locs_R(1:end-1)/fs, RR*1000, '.-'); grid on;
xlabel('时间（秒）'); ylabel('RR间期（毫秒）');
title('RR间期序列（心率变异性）');


% 在命令窗口打印指标
disp('--- 心率变异性指标 ---');
disp(metrics);


% 第二个图形窗口：频域分析和非线性分析
figure('Color','w','Name','心率变异性频域和非线性分析','Position',[120 120 1100 400]);
% 子图1：功率谱密度
subplot(1,3,1);
[pxx,f] = hrv_psd(NN, locs_R(1:end-1)/fs);
plot(f, pxx); xlim([0 0.5]); grid on;
xlabel('频率（赫兹）'); ylabel('功率（毫秒²/赫兹）'); title('心率变异性功率谱密度（Welch方法）');


% 子图2：庞加莱图
subplot(1,3,2);
NN1 = NN(1:end-1); NN2 = NN(2:end);
plot(NN1*1000, NN2*1000, '.'); grid on; axis equal;
xlabel('NN_n（毫秒）'); ylabel('NN_{n+1}（毫秒）');
title(sprintf('庞加莱图（SD1=%.1f毫秒, SD2=%.1f毫秒）', metrics.SD1_ms, metrics.SD2_ms));


% 子图3：归一化LF/HF功率
subplot(1,3,3);
bar([metrics.lfnu, metrics.hfnu]); set(gca,'XTickLabel',{'LFnu','HFnu'}); ylim([0 100]); grid on;
title('归一化LF/HF功率（%）');


% 简单文本面板输出
fprintf('\n房颤样不规则性筛查：%s\n', ternary(metrics.af_screen_positive, '阳性（启发式）', '阴性'));
fprintf('平均NN间期 = %.0f毫秒 | SDNN = %.0f毫秒 | RMSSD = %.0f毫秒 | pNN50 = %.1f%%\n', ...
    metrics.meanNN_ms, metrics.sdnn_ms, metrics.rmssd_ms, metrics.pnn50);
fprintf('LFnu = %.1f%% | HFnu = %.1f%% | RR间期变异系数 = %.2f | 滞后1自相关 = %.2f\n\n', ...
    metrics.lfnu, metrics.hfnu, metrics.cv_rr, metrics.ac_lag1);


end


%% ================== 辅助函数 ==================


function [rr_s, r_times] = generate_rr_series(dur_s, hr_bpm, sdnn_ms, lf_hz, hf_hz, af_mode)
% 生成具有LF/HF调制或房颤样抖动的RR间期序列
mean_rr = 60/hr_bpm; % 平均RR间期（秒）
tstep = 0.25; % 调制时间步长
tmod = 0:tstep:dur_s; % 调制时间向量


if af_mode % 房颤模式
    % 房颤样：更大、更随机的变异性，带有1/f^α噪声
    base = mean_rr + 0.30*randn(size(tmod));
    base = max(0.3, base); % 确保RR间期不小于0.3秒
else % 正常模式
    % 准正弦LF/HF调制
    lf = 0.05*sin(2*pi*lf_hz*tmod + 2*pi*rand);
    hf = 0.03*sin(2*pi*hf_hz*tmod + 2*pi*rand);
    base = mean_rr + lf + hf + 0.01*randn(size(tmod));
end


% 缩放到目标SDNN
rr_interp = interp1(tmod, base, 0:mean_rr:dur_s+10, 'linear', 'extrap');
rr_s = rr_interp(:);
scale = (sdnn_ms/1000) / std(rr_s);
rr_s = (rr_s - mean(rr_s))*scale + mean_rr;
rr_s(rr_s < 0.3) = 0.3; % 确保RR间期不小于0.3秒


r_times = cumsum(rr_s); % 累积RR间期得到R峰时间
r_times = r_times(r_times < dur_s); % 截断到持续时间范围内
rr_s = diff([0; r_times]); % 从R峰时间重新计算RR间期
end


function ecg = synth_ecg_from_rr(t, r_times, fs)
% 简单心搏模板：使用高斯波形叠加模拟P、Q、R、S、T波
% 幅度（毫伏）和宽度调整以符合真实情况
ecg = zeros(size(t)); % 初始化ECG信号
% 模板参数相对于R峰时间（秒）
comp = [ ... %   延迟(秒)   幅度(毫伏)  宽度(秒)
          -0.20       0.10     0.05;   % P波
          -0.04      -0.15     0.015;  % Q波
           0.00       1.00     0.012;  % R波
           0.02      -0.25     0.015;  % S波
           0.25       0.35     0.08];  % T波
for k=1:numel(r_times) % 对每个R峰时间
    for c=1:size(comp,1) % 对每个波形成分
        mu = r_times(k) + comp(c,1); % 成分中心时间
        A  = comp(c,2); % 成分幅度
        w  = comp(c,3); % 成分宽度
        ecg = ecg + A * exp(-0.5*((t-mu)/w).^2); % 添加高斯波形
    end
end
% 添加轻微形态变异性
ecg = ecg + 0.005*randn(size(ecg));
end


function y = hp_filter(x, fs, fc)
% 高通滤波～去除基线漂移
if has_spt() % 检查是否有信号处理工具箱
    [b,a] = butter(2, fc/(fs/2), 'high'); % 设计2阶高通巴特沃斯滤波器
    y = filtfilt(b,a,x); % 零相位滤波
else % 无信号处理工具箱的回退方法
    % 通过移动中位数和移动平均进行简单去趋势
    w = max(1, round(fs*0.6)); % 窗口大小
    x_med = movmedian(x, w); % 移动中位数
    y = x - movmean(x_med, w); % 减去移动平均
end
end


function y = notch_filter(x, fs, f0)
% 电源干扰陷波滤波
Q = 30; % 品质因数
if has_spt() % 检查是否有信号处理工具箱
    w0 = f0/(fs/2); % 归一化中心频率
    bw = w0/Q; % 带宽
    [b,a] = iirnotch(w0, bw); % 设计陷波滤波器
    y = filtfilt(b,a,x); % 零相位滤波
else % 无信号处理工具箱的回退方法
    % 粗略方法：减去窄带正弦拟合
    t = (0:numel(x)-1)/fs; % 时间向量
    ref = sin(2*pi*f0*t); % 参考正弦波
    alpha = (ref*x.')/(ref*ref.'); % 计算最佳幅度
    y = x - alpha*ref; % 减去估计的干扰
end
end


function y = bp_filter(x, fs, f1, f2)
% 带通滤波聚焦QRS波能量
if has_spt() % 检查是否有信号处理工具箱
    [b,a] = butter(3, [f1 f2]/(fs/2), 'bandpass'); % 设计3阶带通巴特沃斯滤波器
    y = filtfilt(b,a,x); % 零相位滤波
else % 无信号处理工具箱的回退方法
    % 移动平均高通强调
    y = x - movmean(x, round(fs*0.15)); % 简单高通滤波
end
end


function tf = has_spt()
% 检测信号处理工具箱是否存在（butter/filtfilt）
tf = exist('butter','file')==2 && exist('filtfilt','file')==2;
end


function [locs, thr] = detect_r_peaks(int_sig, fs)
% 在积分后的QRS增强信号上使用自适应阈值检测
int_sig = int_sig(:); % 确保列向量
% 使用百分位数估计噪声和信号水平
noise_est = prctile(int_sig, 60); % 噪声估计（第60百分位数）
sig_est   = prctile(int_sig, 98); % 信号估计（第98百分位数）
thr = 0.3*noise_est + 0.7*sig_est; % 自适应阈值


% 不应期约200毫秒
ref = round(0.2*fs); % 不应期样本数
cand = find(int_sig > thr); % 找到超过阈值的候选点


if isempty(cand) % 如果没有候选点
    locs = [];
    return;
end


locs = []; % 初始化R峰位置
k = 1; % 初始化索引
while k <= numel(cand) % 遍历所有候选点
    idx = cand(k); % 当前候选点索引
    % 形成一个不应期窗口
    w_end = min(numel(int_sig), idx+ref); % 窗口结束索引
    [~, localmax] = max(int_sig(idx:w_end)); % 在窗口内找到局部最大值
    locs(end+1) = idx + localmax - 1; % 存储局部最大值位置
    % 跳过不应期
    k = find(cand > idx+ref, 1, 'first'); % 找到下一个不应期外的候选点
    if isempty(k), break; end % 如果没有更多候选点，退出循环
end
locs = unique(locs); % 确保位置唯一
end


function locs_ref = refine_r_peaks(ecg, locs, rad)
% 在带通滤波后的ECG上精炼R峰位置为局部最大值
locs_ref = locs; % 初始化精炼后的位置
N = numel(ecg); % ECG信号长度
for i=1:numel(locs) % 对每个检测到的位置
    a = max(1, locs(i)-rad); % 搜索窗口开始
    b = min(N, locs(i)+rad); % 搜索窗口结束
    [~, mx] = max(ecg(a:b)); % 在窗口内找到最大值
    locs_ref(i) = a + mx - 1; % 更新位置
end
% 移除过于接近的峰（<200毫秒）
minsep = round(0.2 * (length(ecg)/max(1,ceil(max(locs)/ (length(ecg)))))); % 最小间隔（未使用）
% 确保严格递增且间隔合适
locs_ref = sort(locs_ref); % 排序
if numel(locs_ref) > 1 % 如果有多个峰
    RR = diff(locs_ref); % 计算RR间期（样本数）
    keep = [true, RR > round(0.2* (length(ecg)/ (ecg_time_length(ecg))))]; % 保留合适间隔的峰
    locs_ref = locs_ref(keep); % 应用筛选
end
end


function T = ecg_time_length(ecg)
% 辅助函数：用于精炼间隔
T = numel(ecg); % 仅用于避免额外参数
end


function [lfnu, hfnu, lf_center, hf_center, lf_power, hf_power, total_power] = hrv_freq(NN, tNN)
% 将NN(t)插值到均匀网格并计算功率谱密度；积分LF/HF频带
% 频带（成人）：VLF 0.003-0.04 Hz（此处忽略），LF 0.04-0.15，HF 0.15-0.40
if numel(NN) < 4 % 检查是否有足够的数据点
    lfnu=NaN; hfnu=NaN; lf_center=NaN; hf_center=NaN;
    lf_power=NaN; hf_power=NaN; total_power=NaN; return;
end
fs_hrv = 4; % 4 Hz插值频率
tu = tNN(1):1/fs_hrv:tNN(end); % 均匀时间网格
x  = interp1(tNN, NN, tu, 'pchip'); % 三次Hermite插值


[pxx, f] = pwelch(x - mean(x), 256, 128, 1024, fs_hrv); % Welch功率谱估计


lf_idx = f>=0.04 & f<0.15; % LF频带索引
hf_idx = f>=0.15 & f<=0.40; % HF频带索引


lf_power = trapz(f(lf_idx), pxx(lf_idx)); % LF绝对功率
hf_power = trapz(f(hf_idx), pxx(hf_idx)); % HF绝对功率
total_power = trapz(f(f>=0.04 & f<=0.40), pxx(f>=0.04 & f<=0.40)); % 总功率


lfnu = 100 * lf_power / (lf_power + hf_power); % 归一化LF功率
hfnu = 100 * hf_power / (lf_power + hf_power); % 归一化HF功率


lf_center = centroid(f(lf_idx), pxx(lf_idx)); % LF中心频率
hf_center = centroid(f(hf_idx), pxx(hf_idx)); % HF中心频率
end


function [pxx,f] = hrv_psd(NN, tNN)
% 辅助函数：显示功率谱密度
if numel(NN) < 4 % 检查是否有足够的数据点
    pxx = zeros(1,512); f = linspace(0,0.5,512); return;
end
fs_hrv = 4; % 4 Hz插值频率
tu = tNN(1):1/fs_hrv:tNN(end); % 均匀时间网格
x  = interp1(tNN, NN, tu, 'pchip'); % 三次Hermite插值
[pxx,f] = pwelch(x - mean(x), 256, 128, 1024, fs_hrv); % Welch功率谱估计
end


function c = centroid(f, p)
% 计算频谱中心频率
if isempty(f) || isempty(p) || sum(p)<=0 % 检查输入有效性
    c = NaN; return;
end
c = sum(f(:).*p(:))/sum(p(:)); % 加权平均频率
end


function [SD1, SD2] = poincare_SD1_SD2(NN)
% SD1：短期变异性，SD2：长期变异性
NN1 = NN(1:end-1); % 前一个NN间期
NN2 = NN(2:end); % 后一个NN间期
diffs = (NN2 - NN1)/sqrt(2); % 差值分量
sums  = (NN2 + NN1)/sqrt(2); % 和分量
SD1 = std(diffs); % SD1：差值分量的标准差
SD2 = std(sums); % SD2：和分量的标准差
end


function ac1 = autocorr_stat(x, lag)
% 计算滞后自相关系数
x = x(:) - mean(x); % 去均值
if numel(x) < lag+1 % 检查是否有足够的数据点
    ac1 = NaN; return;
end
ac1 = sum( x(1:end-lag).*x(1+lag:end) ) / sum( x.^2 ); % 自相关计算
end


function out = ternary(cond, a, b)
% 三元操作符函数
if cond, out=a; else, out=b; end
end