基于MATLAB的心电信号去噪
心电信号(ECG)是诊断心脏疾病的重要工具，但在采集过程中常受到各种噪声干扰。本文将介绍基于MATLAB的心电信号去噪方法，并提供完整的代码实现。
噪声类型分析
- 基线漂移：低频噪声(0.5Hz以下)，由呼吸运动和电极移动引起
- 工频干扰：50/60Hz交流电源干扰
- 肌电噪声：高频噪声(5-2000Hz)，由肌肉活动产生
- 运动伪影：由患者移动引起的不规则干扰
MATLAB去噪方法
1. 基本滤波方法
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
2. 小波变换去噪
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
3. 自适应滤波
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
4. 经验模态分解(EMD)
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
5. 综合去噪流程
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
方法对比分析
| 去噪方法 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| 基本滤波 | 实现简单，计算量小 | 对非平稳噪声效果差，可能扭曲QRS波 | 实时处理，硬件资源有限 |
| 小波变换 | 有效处理非平稳信号，保留特征点 | 小波基和阈值选择影响效果 | 临床分析，研究应用 |
| 自适应滤波 | 对工频干扰去除效果好，实时性强 | 需要参考信号，对其它噪声效果有限 | 工频干扰严重环境 |
| EMD | 完全自适应，无需预设参数 | 计算量大，端点效应问题 | 复杂噪声环境 |
| 综合方法 | 全面处理各类噪声，效果最优 | 实现复杂，计算量大 | 高精度医疗诊断 |
性能评估指标
信噪比(SNR)：衡量信号与噪声的比例
$$SNR = 10 \log_{10}\left(\frac{P_{\text{signal}}}{P_{\text{noise}}}\right)$$均方根误差(RMSE)：评估去噪信号与原始信号的差异
$$RMSE = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(x_i - \hat{x}_i)^2}$$相关系数：评估波形相似度
$$r = \frac{\sum_{i=1}^{N}(x_i - \bar{x})(\hat{x}_i - \bar{\hat{x}})}{\sqrt{\sum_{i=1}^{N}(x_i - \bar{x})^2}\sqrt{\sum_{i=1}^{N}(\hat{x}_i - \bar{\hat{x}})^2}}$$
实际应用建议
预处理：
- 使用高通滤波(0.5Hz)去除基线漂移
- 应用工频陷波滤波器去除50/60Hz干扰
主去噪：
- 首选小波去噪(db6小波，5-6层分解)
- 对肌电噪声严重信号，结合EMD方法
后处理：
- 使用移动平均平滑波形
- 应用幅度归一化
QRS检测：
- 使用Pan-Tompkins等算法检测R峰
- 基于RR间期进行心律分析
结果示例
去噪性能评估:
信噪比(SNR): 18.76 dB
均方根误差(RMSE): 0.0234
相关系数: 0.9826
数据资源
MIT-BIH心律失常数据库：
% 下载并加载MIT-BIH数据 [signal, Fs, tm] = rdsamp('mitdb/100'); ecg = signal(:,1); % 第一导联PTB诊断数据库：
% 访问PTB数据库 ptb = load('ptbdb.mat'); ecg = ptb.data(1).signal(:,1); fs = ptb.data(1).fs;合成ECG生成：
% 生成模拟ECG信号 fs = 360; t = 0:1/fs:10; ecg = ecg_simulation(t);
本文介绍的MATLAB心电信号去噪方法已在临床和研究中广泛应用，实际应用中需根据具体噪声特点和系统要求选择合适的方法组合。
