How can i filter ECG signals with high motion artifact ?

I am working on a mini project where we need to analyse ECG signals for people while running on treadmill, all the codes I've found do not eliminate the high peaks resulted from motion of the person so i cant detect the R peaks. Could anyone help me please because this is the first time for me dealing with Matlab >< am really lost

Answer 1

[Star Strider]

First, do a fft on your data, to see what the frequency of the baseline drift is. It should be much lower than your EKG frequencies. Then, if you have the Signal Processing Toolbox, design a bandpass filter with the low frequency cutoff high enough to eliminate your baseline drift (usually 1 to 5 Hz), and a high frequency cutoff of between about 45 to 100 Hz, depending on your signal. (My filter design procedure is here: How to design a lowpass filter for ocean wave data in Matlab?) You may also need a notch filter to eliminate mains frequency noise. That’s in the MATLAB documentation: Remove the 60 Hz Hum from a Signal. Filter design usually requires some experimentation, so you will have to determine what frequencies are best for your signals. Be sure to sample them with a sampling frequency of at least 250 Hz so that your filters will work most efficiently. If you sample at a lower frequency, you may not be able to design the filter you want.

[Ru'a Khaled]

thank you for answering, but how can i see the freq of baseline drift after taking the fft ? i mean the code that should be written. Another thing, is the sampling freq the same as Fs on the fdatool window ? do i have to assign it as the sampling rate of my data aquisition kit i have used to obtain my samples?

[Star Strider]

My pleasure.

You will have to look at the fft to see the baseline drift signal. It should be very low frequency. See the code between the first two plot figures in the fft link I provided in my Answer. That is the code you will use to calculate and plot your Fourier-transformed data.

The sampling frequency, ‘Fs’ in my code and fdatool, is the sampling frequency of your data acquisition kit. It is the sampling frequency that you design your filters with.

Note that the highest resolvable frequency in a sampled signal is the Nyquist frequency ‘Fn’, and that Fn=Fs/2. Your data acquisition kit most likely has a number of pre-set sampling frequencies. Choose one that is at least 250 Hz, and design your filters using the frequency you are sampling with.

[Ru'a Khaled]

thank you so much for explaining. I did some filtration and this is the best result i could achieve, how can i detect the R peaks while my signal is still like this ?

[Star Strider]

Please attach your signal to your original post, preferably as a .mat file (attach it using the ‘paperclip’ icon). I’ll see what I can come up with as a filter design. I need the sampling times vector as well as the data vector.

[Ru'a Khaled]

ok I've attached it ^^ the sampling rate of acquisition kit I used was 200 :( I am done with samples so I can't change it and re-collect samples

[Star Strider]

It took a while to design a filter I’m happy with, and that eliminates as much noise as possible and still retains the important aspects of the EKG morphology.

Your signal contains motion artifact that is relatively easy to filter out, but also contains considerable muscle (EMG) artifact that is impossible to completely eliminate from the EKG because it has broadband characteristics. There is also some 50 Hz mains frequency noise, but the filter cutoff is low enough to eliminate it.

This is my filter design and signal analysis:

D = load('Ru''a Khaled test1.mat');

t = double(D.b1(:,1));
s = double(D.b1(:,2));

figure(1)
plot(t, s)
grid

Ts = mean(diff(t));                             % Sampling Interval (sec)
Fs = 1/Ts;                                      % Sampling Frequency (Hz)
Fn = Fs/2;                                      % Nyquist Frequency (Hz)
L = length(t);

fts = fft(s)/L;                                 % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;             % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector

figure(2)
plot(Fv, abs(fts(Iv))*2)                        % Plot FFT
grid
xlabel('Frequency')
ylabel('Amplitude')

Wp =  [1.5  20]/Fn;                             % Passband (Normalised)
Ws =  [0.5  45]/Fn;                             % Stopband (Normalised)
Rp =  1;                                        % Passband Ripple (dB)
Rs = 50;                                        % Stopband Ripple (dB)
[n,Wn] = buttord(Wp, Ws, Rp, Rs);               % Filter Order
[b,a]  = butter(n,Wn);                          % Transfer Function
[sos,g] = tf2sos(b,a);                          % Convert To Second-Order-Section For Stability

figure(3)
freqz(sos, 4096, Fs)                            % Filter Bode Plot

sf = filtfilt(sos, g, s);                       % Filter Signal

xlim5 = [ 85   100];                            % Window (Time)
figure(5)                                       % Compare Unfiltered & Filtered Signal
subplot(2,1,1)
plot(t, s)
axis([xlim5    ylim])
title('Unfiltered EKG')
grid
subplot(2,1,2)
plot(t, sf)
axis([xlim5    ylim])
grid
title('Bandpass-Filtered EKG')

I don’t know your experimental setup. You may be able to eliminate at least some of the EMG artifact by using a ‘ground’ electrode, and subtracting the ground signal from your EKG signal in hardware (differential amplifier). Biomedical instrumentation is not trivial to design and implement, and more data are always better than fewer data. Later, you can always ignore data you don’t need, but you can’t go back and create data you do need.

[Ru'a Khaled]

thank you so much this was better than my filtering . Concerning the experimental setup, unfortunately I don't have much time to collect data again :( Can you help me with detecting R waves from this signal in order to calculate instantaneous heart rate ? I know am asking a lot but It's really my first time on Matlab and i don't have much time to finish this so learning from scratch won't help me ><

[Star Strider]

My pleasure.

It took me more than a hour of experimenting this morning before I was finally happy with that filter design.

Detecting R-waves from your signal is going to be a real challenge considering the broadband EMG noise. If you have access to the Wavelet Toolbox, I would rely heavily on that, experimenting first using the Daubechies wavelet family, most probably db2, db3, or db4. It would be worth your time to do a PubMed search on wavelet detection of QRS complexes in noisy EKGs, since yours is quite definitely not the first time this problem has arisen! You will have to experiment, because the skeletal muscle motor unit action potentials (MUAPs) your EKG picked up from the skeletal muscles are similar in magnitude and morphology to the QRS complex (which should come a no surprise since skeletal muscle and the myocardium are all striated muscle and have similar action potentials). This is the reason to always use a ground electrode with an EKG so you can use a simple op-amp differential amplifier to subtract the MUAPs in real time from your EKG signals.) Wavelets are your only real hope of detecting QRS complexes in your signal, and that is going to require some significant experimentation as well. Nevertheless, you should be able to do it once you decide on the appropriate wavelet. Fortunately, the Wavelet Toolbox has some really neat GUIs that can make your task much easier.

[Ru'a Khaled]

oh!! am sorry that you had to spend all that time to help me ..and thank you for your advice ^^

[Star Strider]

My pleasure!

That is typical for filter design. I wanted to be sure I gave you a filter that worked correctly and was as robust as I could design it.

Also, if my Answer solved your problem, please Accept it.