Need help to removing motion (breathing?) artifact from ECG signal

I have this single-channel ECG signal with motion artifacts (I believe from breathing).

I've tried several filters, but none have given me a good output.

Not being an expert, I followed the instructions from Kher, 2019 with this code, but a compatibility issue (since I am using version 2018b) prevents me from obtaining any results; by changing some commands the result is unsuitable:

y1 = load('EKG.mat');
y2= (y1 (:,1)); % ECG signal data
a1= (y1 (:,1)); % accelerometer x-axis data
a2= (y1 (:,1)); % accelerometer y-axis data
a3= (y1 (:,1)); % accelerometer z-axis data
y2 = y2/max(y2);
Subplot (3, 1, 1), plot (y2), title ('ECG Signal with motion artifacts'), grid on
a = a1+a2+a3;
a = a/max(a);
mu= 0.0008;

%Hd = adaptfilt.lms(32, mu); %original command
Hd = dsp.LMSFilter('Length', 32, 'StepSize', mu);

% [s2, e] = filter(Hd, a, y2); % original command, don't work in 2018b version

[s2, e] = Hd(a, y2); % command adapted 

fig = figure 
subplot (3, 1, 2)
plot (s2)
title ('Noise (motion artifact) estimate')
grid on
subplot (3, 1, 3)
plot (e)
title ('Adaptively filtered/ Noise free ECG signal')
grid on

I also tried filtering in this other way, but the result is very poor.

ecg_signal = load('EKG.mat');
Fs = 256;
t = (0:length(ecg_signal)-1) / Fs;

fc = 45; % cut frequency
[b, a] = butter(4, fc / (Fs / 2), 'low');

% filter
ecg_filtered = filtfilt(b, a, ecg_signal);

With simple low-pass or high-pass filters I wasn't able to obtain at least acceptable results.

If anyone can help me?

thank you in advance

[Mathieu NOE]

hello

you provided only the ECG signal in your data

No one can work on your first code without the accel data

[Simone]

it was just an unsuccessful attempt. I haven't accel data.

how can I remove the artifacts visible in the figure?

[Mathieu NOE]

can you show us a section with artifacts vs a clean section ? tx

[Simone]

the entire trace of this subject shows these artifacts. this one in the photo is without these artifacts

Answer 1

[Star Strider]

There does not appear to be much noise at all, however the small baseline variation and high-frequency noise can be eliminated with an appropriiately-designed bandpass filter.

Try this —

LD = load('EKG.mat')

LD = struct with fields:
    data256hz: [118752x1 table]

EKG = LD.data256hz{:,1};
Fs = 256;                                                               % No Time Vector Supplied, So Create One
t = linspace(0, numel(EKG)-1, numel(EKG)).'/Fs;

figure
plot(t, EKG)
grid
xlim('tight')
xlabel('Time (s)')
ylabel('Amplitude (V)')
title('Unprocessed EKG Signal')

figure
plot(t, EKG)
grid
xlim([0 5])
xlabel('Time (s)')
ylabel('Amplitude (V)')
title('Unprocessed EKG Signal (Detail)')

[FTs1,Fv] = FFT1(EKG,t);

figure
plot(Fv, abs(FTs1)*2)
grid
xlim('tight')
xlabel('Frequency (Hz)')
ylabel('Magnitude (V)')
title('Fourier Transform')

figure
plot(Fv, abs(FTs1)*2)
grid
xlim([0 10])
xlabel('Frequency (Hz)')
ylabel('Magnitude (V)')
title('Fourier Transform (Detail)')

EKG = [ones(100,1)*EKG(1); EKG];                                        % Pre-Appeend 'ones' Veector To Eliiminate Starting Transient
EKG_filt = bandpass(EKG, [1, 45], Fs, 'ImpulseResponse','iir');         % Filter With Elliptic IIR Bandpass Filter
EKG_filt = EKG_filt(101:end);

figure
plot(t, EKG_filt)
grid
xlim([0 5])
xlabel('Time (s)')
ylabel('Amplitude (V)')
title('Bandpass-Filtered EKG Signal (Detail)')

function [FTs1,Fv] = FFT1(s,t)
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

Use the Fourier transform results to guide the filter design.

.

[Simone]

thanks a lot, this seems accurate without cutting out too much of the signal. Even though it wasn't very noisy, those few artifacts still made it hard to calculate the QRS complex

where can I find the FFT1.m command? It's not in the installed Signal Processing Toolbox

[Star Strider]

As always, my pleasure!

The ‘FFT1’ function is one I wrote, and appended at the end of the code I posted:

function [FTs1,Fv] = FFT1(s,t)

% Arguments:

% s: Signal Vector Or Matrix

% t: Associated Time Vector

t = t(:);

L = numel(t);

if size(s,2) == L

s = s.';

end

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));

Fv = Fs*(0:(NFFT/2))/NFFT;

% Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

FTs1 = FTs(Iv,:);

end

Feel free to copy it and use it. If you save it as a separate file on your MATLAB search path, name that file FFT1.m.

.

Answer 2

[Image Analyst]

I'm sure it's been done before. Rather than me just guessing, I suggest you search the medical literature in PubMed to find articles that show how it's done.

https://pubmed.ncbi.nlm.nih.gov/?term=ecg+breathing