心电信号ECG去噪：Matlab实现低通滤波与小波分解结合-程序员充电站

心电信号ECG去噪，Matlab程序，使用低通滤波和小波分解结合。先去除高于80Hz的高频噪声，再去除高于50Hz的噪声和工频干扰等，最后去除基线漂移。具体工作如下： 1、读取一段ECG信号，采样率为200Hz，时长10秒，画出时域、频域图； 2、构造低通滤波器，截止频率为80Hz，去除ECG高频噪声； 3、使用小波分解，进一步去除高于50Hz的噪声和工频干扰等； 4、使用小波分解去除基线漂移，程序结果如图。程序运行如遇问题可帮忙远程调通。换成自己的ECG信号时需要根据信号情况调整程序参数。

在生物医学信号处理领域，心电信号（ECG）的准确分析至关重要。然而，ECG信号在采集过程中常常受到各种噪声的干扰，影响其后续分析与诊断。本文将介绍如何通过Matlab程序，利用低通滤波和小波分解相结合的方法，有效地去除ECG信号中的噪声。

1. 读取ECG信号并绘制时域、频域图

首先，我们需要读取一段ECG信号。假设信号采样率为200Hz，时长10秒。在Matlab中，我们可以这样实现：

fs = 200; % 采样率 t = 0:1/fs:(10 - 1/fs); % 时间向量 % 这里假设已有生成好的ECG信号数据存储在变量ecg_signal中 % 如果是从文件读取，可以使用load函数，例如load('ecg_data.mat'); % 假设文件中ECG信号变量名为ecg_signal ecg_signal = sin(2*pi*10*t) + 0.5*randn(size(t)); % 这里简单模拟一个ECG信号 % 绘制时域图 figure; subplot(2,1,1); plot(t, ecg_signal); xlabel('时间 (s)'); ylabel('幅值'); title('ECG信号时域图'); % 绘制频域图 n = length(ecg_signal); f = (0:n - 1)*(fs/n); ecg_fft = abs(fft(ecg_signal)); subplot(2,1,2); plot(f, ecg_fft(1:length(f))); xlabel('频率 (Hz)'); ylabel('幅值'); title('ECG信号频域图');

上述代码中，首先定义了采样率fs和时间向量t。这里我们简单模拟了一个ECG信号，实际应用中可从文件读取真实数据。通过plot函数绘制了时域图，直观展示信号随时间的变化。利用fft函数进行快速傅里叶变换，得到信号频域信息并绘制频域图，以便观察信号的频率分布情况。

2. 构造低通滤波器去除高频噪声

接下来，构造截止频率为80Hz的低通滤波器，去除高于80Hz的高频噪声。

fc1 = 80; % 截止频率 [b1, a1] = butter(5, fc1/(fs/2)); % 5阶巴特沃斯低通滤波器 ecg_filtered1 = filter(b1, a1, ecg_signal); % 绘制滤波后时域图 figure; subplot(2,1,1); plot(t, ecg_filtered1); xlabel('时间 (s)'); ylabel('幅值'); title('80Hz低通滤波后ECG信号时域图'); % 绘制滤波后频域图 n1 = length(ecg_filtered1); f1 = (0:n1 - 1)*(fs/n1); ecg_fft1 = abs(fft(ecg_filtered1)); subplot(2,1,2); plot(f1, ecg_fft1(1:length(f1))); xlabel('频率 (Hz)'); ylabel('幅值'); title('80Hz低通滤波后ECG信号频域图');

这里使用butter函数设计了一个5阶巴特沃斯低通滤波器，截止频率为fc1。filter函数将设计好的滤波器应用到ECG信号上，得到滤波后的信号ecg_filtered1。同样绘制了滤波后信号的时域图和频域图，从频域图中可以清晰看到高于80Hz的高频成分已被有效抑制。

3. 利用小波分解去除高于50Hz的噪声和工频干扰

进一步利用小波分解来去除高于50Hz的噪声和工频干扰等。

% 小波分解 wname = 'db4'; % 选用db4小波 level = 5; % 分解层数 [c, l] = wavedec(ecg_filtered1, level, wname); % 获取高频系数 for i = 1:level detcoeffs{i} = detcoef(c, l, i); end % 获取低频系数 approxcoeff = appcoef(c, l, wname, level); % 设定阈值去除高频噪声 for i = 1:level thresh = wthrmngr('dw1ddenoLVL', 'rigrsure', detcoeffs{i}); detcoeffs{i} = wthresh(detcoeffs{i}, 'h', thresh); end % 重构信号 new_c = [approxcoeff; detcoeffs{:}]; ecg_filtered2 = waverec(new_c, l, wname); % 绘制小波分解后时域图 figure; plot(t, ecg_filtered2); xlabel('时间 (s)'); ylabel('幅值'); title('小波分解去除50Hz以上噪声后ECG信号时域图');

代码中选用db4小波进行5层小波分解，wavedec函数得到分解系数c和长度向量l。通过detcoef和appcoef分别获取高频和低频系数。利用wthrmngr函数基于rigrsure准则计算阈值，对高频系数进行阈值处理去除噪声。最后使用waverec函数重构信号，得到进一步滤波后的ecg_filtered2，并绘制其时域图。

4. 使用小波分解去除基线漂移

最后，再次利用小波分解去除基线漂移。

% 再次小波分解去除基线漂移 wname2 = 'db6'; % 选用db6小波 level2 = 8; % 分解层数 [c2, l2] = wavedec(ecg_filtered2, level2, wname2); % 获取低频系数（基线漂移主要在低频部分） approxcoeff2 = appcoef(c2, l2, wname2, level2); % 去除低频系数（即去除基线漂移） new_c2 = c2; new_c2(l2(1):end) = 0; ecg_filtered3 = waverec(new_c2, l2, wname2); % 绘制去除基线漂移后时域图 figure; plot(t, ecg_filtered3); xlabel('时间 (s)'); ylabel('幅值'); title('去除基线漂移后ECG信号时域图');

这里选用db6小波进行8层小波分解，获取低频系数，将低频系数置零（因为基线漂移主要存在于低频部分），然后重构信号得到去除基线漂移后的ecg_filtered3，并绘制其时域图。