心电信号降噪方法（滤波器/移动平均/小波等，MATLAB环境）

news2024/9/17 3:23:06

对于一个正常的、完整的心动周期，对应的心电图波形如下图所示，各个波形都对应着心脏兴奋活动的生理过程，包含P波，PR段，QRS波群，ST段，T波，U波。

（1）P波心电图中，起始的P波对应心房的兴奋过程，一个完整的心动周期起始于窦房结，P波的起始，代表窦房结将心脏兴奋传递至心房使得心房全部受到激励的过程。P波的方向及形态与兴奋在心房内的传播途径密切相关，在各个导联中，P波形态不尽一致，一般来说，P波顶部圆钝，其幅值不高于0.25mV，时限宽度在0.11s以内。

（2）QRS复合波在心电图中，QRS复合波也称为QRS波群，表征着心室除极的心脏电位变化状态，由于心室各部位受到激励顺序有先后，因此QRS复合波是代表心室几个部分兴奋过程的综合波形，是由方向向下的Q波，波形宽大且幅值较高的R波，以及波形方向向下的S波组合形成的。QRS复合波在心电图中波形最为明显，有波形宽，斜率大的特点，因此最容易被识别，常作为心电信号研究的基础，应用广泛。一般正常情况下，QRS波的时限宽度在0.06-0.10s。而在不同的导联中，Q波，R波，S波的振幅、方向不尽一致，致使QRS波形态各异。一般来说Q波振幅小于同一导联中R波的1/4~1/2。在双极肢体导联I、II、III中，R波的振幅一般不超过1.5mV，2.5mV，2.0mV。在标准导联中，S波的振幅的变化量不超过0.6mV。

（3）S-T段在QRS波群的终点到与之相继出现的T波起点中间的水平连线是S-T段，表征着从心室除极结束至缓慢复极的心电活动状态。一般来说，S-T段正常情况下其时限宽度在0.10-0.14s，接近于基线，偏移量较小，一般下移量在0.05mV以内，ST段的抬高量在不同导联中有不同范围，例如在V3导联中抬高量在0.5mV内，ST段的偏移量超过正常范围针对心肌类疾病诊断意义重大。

（4）T波在QRS波群之后出现的顶部圆钝的波形是T波，表征心室复极的状态，通常T波与其同一心动周期的QRS主波方向统一，但也存在T波倒置，双向等情况，其所占时限较长，一般为0.05-0.250s，在不同导联中其幅值也不尽相同，一般应大于与其相邻的R波幅值的1/10。

（5）Q-T间期Q-T间期是从Q波的起始点到T波的终点之间的距离，表征的是心室除极开始到复极完成的心脏活动过程。此段间期与心搏速率相关联，心率越慢，间期越长，相反则越短。一般情况下，若心率在70次/分左右，那么Q-T间期在0.32-0.44s左右，超过正常Q-T间期最高值0.03s被称为Q-T间期显著延长。

（6）U波U波出现在T波之后，是因心室复极化而形成，波形宽而低平，振幅较小，一般不超过0.05mV，时限宽度在0.02~0.04s。U波一般与T波方向相同，其振幅约在0.03mV左右，U波不一定出现在所有导联中，在胸导联中其波形较为明显。

鉴于此，采用滤波器/移动平均/小波等方法进行心电信号降噪，运行环境为MATLAB，运行代码如下：

clc
clear all;
close all;

folderType = 'AVRT';
fileName = strcat('JS10493');
% folderType = 'AF';
% fileName = strcat('JS10361');
matFileName = strcat(fileName, '.mat');
heaFileName = strcat(fileName, '.hea');

data_path = fullfile('data/raw/', folderType, matFileName);

data = load(data_path);

ecg_orig = data.val(1, :);
ecg_orig(isinf(ecg_orig)|isnan(ecg_orig)) = 0;


fs = 125;
N = length(ecg_orig);
t = (0:N-1) / fs;
f = (0:N-1) * (fs / N);
frequencies = f(1:N/2+1);

%RESAMPLE SIGNAL To 125Hz
% [ecg, t] = resample_signal_125(ecg_orig, t_orig);

%Assigning new parameter
% fs = 125;
% N = length(ecg);
% f = (0:N-1) * (fs / N);


%Applying high pass and low pass filters
hl_filter = hpf_lpf(ecg_orig, fs);
pl_filter = powerline_removal(hl_filter, fs);
movSignal = moving_avg(pl_filter, 5);
waveletSignal = wavelet_denoise(movSignal, 5);

% Create a single figure with subplots for each signal
figure;

% Original ECG
subplot(5, 1, 1);
plot(t, ecg_orig);
xlabel('Time (s)');
ylabel('Amplitude (mV)');
title('Original Signal');

% HPF + LPF Filtered
subplot(5, 1, 2);
plot(t, hl_filter);
xlabel('Time (s)');
ylabel('Amplitude (mV)');
title('HPF + LPF Signal');

% Powerline Removal
subplot(5, 1, 3);
plot(t, pl_filter);
xlabel('Time (s)');
ylabel('Amplitude (mV)');
title('Powerline Removal');

% Moving Average applied
subplot(5, 1, 4);
plot(t, movSignal);
xlabel('Time (s)');
ylabel('Amplitude (mV)');
title('Moving Average Applied');

% Wavelet Denoised
subplot(5, 1, 5);
plot(t, waveletSignal);
xlabel('Time (s)');
ylabel('Amplitude (mV)');
title('Wavelet Denoised');


out = waveletSignal;

figure
subplot(4, 1, 1);

[Pxx,F] = periodogram(ecg_orig,[],length(ecg_orig),fs);
plot(F,10*log10(Pxx))

title('Original Signal Spectrum');
xlabel('Frequency');

subplot(4, 1, 2);

[Pxx_out,F_out] = periodogram(out,[],length(out),fs);   
plot(F,10*log10(Pxx_out))

title('Output Signal Spectrum');
xlabel('Frequency');

subplot(4, 1, 3)
input_fft = fft(ecg_orig)/N;
input_fft = 2 * abs(input_fft(1:N/2+1));
plot(frequencies, input_fft)

title('Original Signal FFT');
xlabel('Frequency');

subplot(4, 1, 4)
output_fft = fft(out)/N;
output_fft = 2 * abs(output_fft(1:N/2+1));
plot(frequencies, output_fft)

title('Output Signal FFT');
xlabel('Frequency');

output_folder_path = fullfile('data/processed/', folderType);
if ~exist(output_folder_path, 'dir')
  mkdir(output_folder_path);
end

output_path = fullfile('data/processed/', folderType, matFileName);
save(output_path, 'out');

source_file_path_hea = fullfile('data/raw/', folderType, heaFileName);
destination_file_path_hea = fullfile("data/processed/", folderType, heaFileName);

copyfile(source_file_path_hea, destination_file_path_hea);

%完整代码：https://mbd.pub/o/bread/mbd-ZZiUmJhw

%assess_signals(ecg, out, fs)