How to Remove the baseline wander from ECG?

Hello there! I've noticed a similar question posted in the community, so I kindly request that you refrain from providing a generic response like "Check this." Instead, I'm reaching out for your assistance in effectively eliminating the baseline wander from my signal. Despite my efforts to remove it, the resulting plot remains unchanged, and I haven't observed any improvement. I would greatly appreciate your help with the code to address this issue. Thank you in advance for your support!

load sig_ecg.mat

fs = 300;
t=linspace(0,length(sig)/fs,length(sig));
figure()
plot(t,sig/max(sig));
grid
title('Signal in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');

figure()
sig=sig-mean(sig); % Per visualizzare meglio lo spettro
f=linspace(-fs/2,fs/2,length(sig));
plot(f,fftshift(abs(fft(sig))));
grid
title('Signal in Frequency Domain'); 
xlabel('Frequency (Hz)');
ylabel('Amplitude');
Lo spettro del segnale mostra un consistente picco a 42Hz che è probabilmente correlato ad una "power line interference", dovuta ad artefatti da strumentazione. L'approccio più semplice per rimuovere quest'interferenza è di applicare un filtraggio IIR stop-band centrato a 42Hz. Inoltre notiamo la presenza di una baseline wander.
F = sig;
f = gca;
EKG = findobj(f, 'Type','line');
t = EKG.XData;                                              % Recover Data
s = EKG.YData;
L = numel(t);                                               % Signal Length
Ts = mean(diff(t));
St = std(diff(t));                                          % Test For Uniform Sampling
tr = linspace(min(t), max(t), L);                           % New, Regularly-Sampled Time Vector
Tsr = mean(diff(tr)) * 1E-3;                                % Convert To Seconds
sr = resample(s, tr);                                       % Resample To Regularly-Sampled Signal
Fsr = 1/Tsr;                                                % Sampling Frequency (Hz)
Fnr = Fsr/2;
figure
plot(tr, sr)
grid
FTsr = fft(sr-mean(sr))/L;                                  % Fourier Transform Of Mean-Corrected Signal
Fv = linspace(0, 1, fix(L/2)+1)*Fnr;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTsr(Iv))*2)
grid
Fs = Fsr;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Ws = 0.5/Fn;                                                % Passband Frequency Vector (Normalised)
Wp = 1.5/Fn;                                                % Stopband Frequency Vector (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs = 150;                                                   % Stopband Attenuation (dB)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'high');                         % Default Here Is A Lowpass Filter
[sos,g] = zp2sos(z,p,k);                                    % Use Second-Order-Section Implementation For Stability
sr_filtered = filtfilt(sos,g,sr);                           % Filter Signal (Here: ‘sr’)
figure
freqz(sos, 2^14, Fs)                                        % Bode Plot Of Filter
set(subplot(2,1,1), 'XLim',[0 Fn])                          % Optional, Change Limits As Necessary
set(subplot(2,1,2), 'XLim',[0 Fn])                          % Optional, Change Limits As Necessary
figure
plot(tr, sr_filtered)
grid
title('Signal in Frequency Domain with removal of baseline wander'); 
xlabel('Frequency (Hz)');
ylabel('Amplitude');

Answer 1

[Star Strider]

I recognise my code, however that version was not intended to work with a signal such as yours.

Try this instead —

load sig_ecg.mat

Fs = 300;
L = numel(sig);
t = linspace(0, L-1, L)/Fs;
% t=linspace(0,length(sig)/fs,length(sig));
figure()
plot(t,sig/max(sig));
grid
title('Original Signal in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');

Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
NFFT = 2^nextpow2(L)

NFFT = 65536

FTsig = fft((sig-mean(sig)).*hann(L), NFFT);
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTsig(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')

% xlim([0 5])

sig = bandstop(sig, [48 52], Fs, 'ImpulseResponse','iir');  % Remove 50 Hz Mains Interference

% figure()
% sig=sig-mean(sig); % Per visualizzare meglio lo spettro
% f=linspace(-fs/2,fs/2,length(sig));
% plot(f,fftshift(abs(fft(sig))));
% grid
% title('Signal in Frequency Domain'); 
% xlabel('Frequency (Hz)');
% ylabel('Amplitude');
% % Lo spettro del segnale mostra un consistente picco a 42Hz che è probabilmente correlato ad una "power line interference", dovuta ad artefatti da strumentazione. L'approccio più semplice per rimuovere quest'interferenza è di applicare un filtraggio IIR stop-band centrato a 42Hz. Inoltre notiamo la presenza di una baseline wander.
% F = sig;
% f = gca;
% EKG = findobj(f, 'Type','line');
% t = EKG.XData;                                              % Recover Data
% s = EKG.YData;
% L = numel(t);                                               % Signal Length
% Ts = mean(diff(t));
% St = std(diff(t));                                          % Test For Uniform Sampling
% tr = linspace(min(t), max(t), L);                           % New, Regularly-Sampled Time Vector
% Tsr = mean(diff(tr)) * 1E-3;                                % Convert To Seconds
% sr = resample(s, tr);                                       % Resample To Regularly-Sampled Signal
% Fsr = 1/Tsr;                                                % Sampling Frequency (Hz)
% Fnr = Fsr/2;
% figure
% plot(tr, sr)
% grid
% FTsr = fft(sr-mean(sr))/L;                                  % Fourier Transform Of Mean-Corrected Signal
% Fv = linspace(0, 1, fix(L/2)+1)*Fnr;
% Iv = 1:numel(Fv);
% figure
% plot(Fv, abs(FTsr(Iv))*2)
% grid
% Fs = Fsr;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Ws = 0.5/Fn;                                                % Passband Frequency Vector (Normalised)
Wp = 2.5/Fn;                                                % Stopband Frequency Vector (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs = 150;                                                   % Stopband Attenuation (dB)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'high');                         % Default Here Is A Lowpass Filter
[sos,g] = zp2sos(z,p,k);                                    % Use Second-Order-Section Implementation For Stability
sig_filtered = filtfilt(sos,g,sig);                           % Filter Signal (Here: isrg)
figure
freqz(sos, 2^14, Fs)                                        % Bode Plot Of Filter
set(subplot(2,1,1), 'XLim',[0 50])                          % Optional, Change Limits As Necessary
set(subplot(2,1,2), 'XLim',[0 50])                          % Optional, Change Limits As Necessary

% figure
% plot(tr, sr_filtered)
% grid
% title('Signal in Frequency Domain with removal of baseline wander'); 
% xlabel('Frequency (Hz)');
% ylabel('Amplitude');

figure
plot(t, sig_filtered)
grid
title('Filtered Signal in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');

This removes the baseline drift, and also eliminates the 50 Hz mains interference with a separate filter. The bandstop function was introduced in R2018a, so you should have it. You can also substitute my filter code with a highpass call, similar to the syntax I used in the bandstop call (although only the 0.5 Hz cutoff frequency, not the 2.5 Hz passband frequency, would need to be provided).

I commented-out parts of the code that are either not appropriate or not helpful for your signal, since this was originally intended to recover data from a .fig file with a time vector that did not have constant, uniform sampling intervals. (I am not certain why I used a two-sided Fourier transform, since I usually use one-sided Fourier transforms for signal processing problems, as I have in revising this.)

.

[HelpAStudent]

Thank you very much that's very helpfull.

Why do you need to calculate the NFFT number ?

And do you know from this point if other filtering is needed?

[Star Strider]

As always, my pleasure!

The ‘NFFT’ calculation has to do with the way the Fourier transform works. The lenght of the fft can be anything greater than or equal to the length of the signal vector, however it is most efficient if the length is a power of 2. See the Wikipedia article on the Fast fourier transform for a discussion and illustration of the reason.

Further filtering does not appear to be necessary at this point. The signal appears to be as clean as it needs to be. There is some remaining high-frequency noise, however it is not significant, so the bandstop filter could be converted to a lowpass filter with a single cutoff frequency of 45 Hz if you want to, with perhaps a slightly better result. The EKG deflections are all visible, and though while not absolutely normal (there are some characteristics that concern me) I cannot identify any specific pathology in it.