基于MATLAB的心电信号去噪

function clean_ecg = basic_filtering(ecg, fs)
    % 带通滤波 (0.5-100Hz)
    [b, a] = butter(3, [0.5, 100]/(fs/2), 'bandpass');
    filtered = filtfilt(b, a, ecg);

    % 陷波滤波去除工频干扰
    wo = 50/(fs/2);  % 工频50Hz
    bw = wo/35;      % 带宽
    [b, a] = iirnotch(wo, bw);
    clean_ecg = filtfilt(b, a, filtered);

    % 绘图
    figure('Name', '基本滤波效果');
    subplot(2,1,1); plot(ecg); title('原始ECG信号');
    subplot(2,1,2); plot(clean_ecg); title('基本滤波后信号');
end

function clean_ecg = wavelet_denoise(ecg, wavelet_name, level)
    % 小波分解
    [C, L] = wavedec(ecg, level, wavelet_name);

    % 计算阈值
    thr = thselect(ecg, 'rigrsure');

    % 应用软阈值处理细节系数
    sorh = 's'; % 软阈值
    keepapp = 1; % 保留近似系数
    denoised = wdencmp('gbl', C, L, wavelet_name, level, thr, sorh, keepapp);

    clean_ecg = denoised;

    % 绘图
    figure('Name', '小波去噪效果');
    subplot(level+2,1,1); plot(ecg); title('原始ECG信号');
    for i = 1:level
        subplot(level+2,1,i+1);
        plot(wrcoef('d', C, L, wavelet_name, i));
        title(['第' num2str(i) '层细节系数']);
    end
    subplot(level+2,1,level+2); plot(clean_ecg); title('小波去噪后信号');
end

function clean_ecg = adaptive_filtering(ecg, fs, mu)
    % 生成参考噪声
    t = (0:length(ecg)-1)/fs;
    ref_noise = sin(2*pi*50*t)'; % 50Hz工频干扰

    % LMS自适应滤波
    order = 32; % 滤波器阶数
    lms = dsp.LMSFilter(order, 'StepSize', mu, 'Method', 'Normalized LMS');
    [clean_ecg, ~, ~] = lms(ref_noise, ecg);

    % 绘图
    figure('Name', '自适应滤波效果');
    subplot(2,1,1); plot(ecg); title('原始ECG信号');
    subplot(2,1,2); plot(clean_ecg); title('自适应滤波后信号');
end

function clean_ecg = emd_denoise(ecg)
    % EMD分解
    imf = emd(ecg);

    % 选择重构分量 (去除高频噪声)
    num_imf = size(imf, 2);
    recon_imf = 2:num_imf-2; % 去除前2个和后2个IMF

    % 重构信号
    clean_ecg = sum(imf(:, recon_imf), 2);

    % 绘图
    figure('Name', 'EMD分解');
    for i = 1:size(imf, 2)
        subplot(size(imf,2),1,i);
        plot(imf(:,i));
        title(['IMF ' num2str(i)]);
    end

    figure('Name', 'EMD去噪效果');
    subplot(2,1,1); plot(ecg); title('原始ECG信号');
    subplot(2,1,2); plot(clean_ecg); title('EMD去噪后信号');
end

function clean_ecg = comprehensive_ecg_denoising(ecg, fs)
    % 步骤1: 去除基线漂移
    [b, a] = butter(3, 0.5/(fs/2), 'high');
    ecg1 = filtfilt(b, a, ecg);

    % 步骤2: 小波去噪
    clean_ecg = wavelet_denoise(ecg1, 'db6', 5);

    % 步骤3: 工频干扰去除
    wo = 50/(fs/2);
    bw = wo/35;
    [b, a] = iirnotch(wo, bw);
    clean_ecg = filtfilt(b, a, clean_ecg);

    % 步骤4: 平滑处理
    clean_ecg = movmean(clean_ecg, 5);

    % 性能评估
    evaluate_performance(ecg, clean_ecg);
end

function evaluate_performance(original, denoised)
    % 计算SNR
    noise = original - denoised;
    signal_power = mean(original.^2);
    noise_power = mean(noise.^2);
    snr = 10*log10(signal_power/noise_power);

    % 计算RMSE
    rmse = sqrt(mean((original - denoised).^2));

    % 计算相关系数
    corr_coef = corrcoef(original, denoised);
    r = corr_coef(1,2);

    % 显示结果
    fprintf('去噪性能评估:\n');
    fprintf('信噪比(SNR): %.2f dB\n', snr);
    fprintf('均方根误差(RMSE): %.4f\n', rmse);
    fprintf('相关系数: %.4f\n', r);

    % 绘制对比图
    figure('Name', '去噪效果对比');
    subplot(3,1,1); plot(original); title('原始ECG信号');
    subplot(3,1,2); plot(denoised); title('去噪后ECG信号');
    subplot(3,1,3); plot(original - denoised); 
    title(['噪声信号 (SNR=' num2str(snr, '%.2f') 'dB)']);
end

%% ECG信号去噪MATLAB实现
clear; close all; clc;

% 1. 加载ECG数据
load('ecg_data.mat'); % 替换为您的ECG数据
fs = 360; % 采样率 (Hz)
t = (0:length(ecg)-1)/fs;

% 2. 添加噪声 (如果数据是干净的)
noise_free = ecg; % 保存干净信号
ecg = ecg + 0.1*randn(size(ecg)); % 高斯白噪声
ecg = ecg + 0.05*sin(2*pi*50*t)'; % 50Hz工频干扰
ecg = ecg + 0.02*sin(2*pi*0.3*t)'; % 基线漂移

% 3. 应用去噪方法
method = 5; % 1-基本滤波, 2-小波去噪, 3-自适应滤波, 4-EMD, 5-综合方法

switch method
    case 1
        clean_ecg = basic_filtering(ecg, fs);
    case 2
        clean_ecg = wavelet_denoise(ecg, 'db6', 6);
    case 3
        clean_ecg = adaptive_filtering(ecg, fs, 0.01);
    case 4
        clean_ecg = emd_denoise(ecg);
    case 5
        clean_ecg = comprehensive_ecg_denoising(ecg, fs);
end

% 4. 特征点检测 (可选)
if exist('noise_free', 'var')
    % 使用去噪后的信号进行QRS检测
    [qrs_amp_raw, qrs_i_raw] = pan_tompkin(clean_ecg, fs, 0);

    figure('Name', 'QRS波检测');
    plot(t, clean_ecg); hold on;
    plot(t(qrs_i_raw), clean_ecg(qrs_i_raw), 'ro');
    title('QRS波检测结果');
    xlabel('时间 (s)');
    ylabel('幅度');
    legend('ECG信号', 'R峰位置');
end