Butterworth filter for ecg signal in healthy people

Rp = 1;             % passband ripple
Rs = 40;            % stopband ripple
Fs = 5000;          % sampling frequency
Fn = Fs/2;          % nyquist frequency
Fpass = [0.05 150]; % passband
Fstop = [0.01 200]; % stopband
Wp = Fpass/Fn;      % normalized
Ws = Fstop/Fn;

% Create a butterworth filter 
[N,Wn] = buttord(Wp,Ws,Rp,Rs);  % Returns the lowest order
[B,A] = butter(N,Wn);           % Butterworth filter design
[sos1,g1] = tf2sos(B,A);        % Convert digital filter transfer to second-order section form

ecg = filtfilt(sos1, g1, ecg_raw);
time = 1:length(ecg);

subplot(2,1,1)
plot(time,ecg_raw);
title('Raw ecg')
subplot(2,1,2)
plot(time,ecg)
title('Ecg filtered')

Hi.

Im trying to filter a ecg signal using a butterworth filter. Been looking a lot of postes trying to figure it out. Im new in signal in signal processing so it's possible that i dont fully understand the processing and the parametres which is required doing it.

can you point me the right direction?

[Mathieu NOE]

hello

can you share the ecg data as well ? otherwise it's difficult to test your code

tx

[Christian Otte]

Hi

Thanks for answering

My data file is 2 large sharing, also after zipping it. Or else i would be happy to share it?

[Mathieu NOE]

hello

yes you can zip it

if it's very big maybe we can start to work on the problem on a smaller extract of the data

also at what sampling rate are you working ?

cheers

[Christian Otte]

I would love your input the signal processing.

It's is a ecg signal at a healthy person. It's sampled at 5000Hz across the heart.

Answer 1

[Star Strider]

The filter may be unstable, and if so, may not produce any output. Another possibility is that one or more of the elements of the data are NaN (although unlikely).

The best option is probably to use the bandpass function, specifically:

Fs = 5000;
Fpass = [1 100];
ecg = bandpass(ecg_raw, Fpass, Fs);

Although if there is no significant baselline wander, using a lowpass filter would work.

The signal appears to have a significant noise component, however it is not possible to determine what that is. The way to determine that is to calculate and plot the Fourier transform of the signal —

ecg_raw = data(:,3);
L = numel(ecg_raw);
time = linspace(0, L, L)/Fs;
N = 2^nextpow2(L);
FTs = fft(ecg_raw,N)/L;
Fv = linspace(0, 1, N/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTs(Iv))*2)
grid

If this is 50-60 Hz supply frequency noise, it can be eliminated with a notch filter (the bandstop function). If it is broadband noise, the only effective ways to deal with it are wavelet denoising or the Savitzky-Golay filter (the sgolayfilt function).

[Christian Otte]

Hi Star Strider.

Thanks for taking your time helping me.

I searched my data for any NaN elements, which it didn't have.

I also tried the bandpass function, with a result which visual seem to be the same as the raw data?

It's possible there is a lot os noise, because we had to have the the power cable for the tilting table pluged in during the recordings. So would my next move be to amply a notch filter to the raw data and afterwards filter the data with the butterworth filter?

I'm not comfortable with the Fourier transform, so i'm not quiet sure what to espext from the plot?

[Star Strider]

My pleasure!

The Fourier transform will tell you much about your signal, so please run the code I posted with at least a representative sample of your data to understand what is actually in your data. I can help you interpret it either with a sample of your data (that I can use with my code) or an image of the Fourier transform.

The EKG has a bandwidth of 0 to 100 Hz, although a ‘normal’ EKG (no arrhythmias or other abnormalities) can be filtered with a lowpass filter with a passband of 0 to 45 Hz without losing any significant amount of information. That would also eliminate power line noise, however it would not eliminate broadband noise. A bandpass filter is necessary if the EKG has significant baseline variation. I do not see that in the signal you posted, however there is not enough detail in that image for me to be certain of that.

[Christian Otte]

If it's not too much asking, I would love your interpretation of the ECG signal.

I've made a sample, from the raw data, of 10 seconds with a sampling frequency at 5000Hz. It should be representative thus it's ECG baseline recording from a supposedly healthy person.

Thanks again

[Star Strider]

Not too much at all!

After seeing the details of the signal and the Fourier transform plot, I decided to design a completely new FIR filter that eliminates the baselline drift, the power frequency signal at 50 Hz, and the harmonic at 100 Hz. It is a long filter, however extrremely efficient, and it is only necessary to run it once for each signal.

Try this —

LD = load('ECG_sample.mat');
time_index = LD.time_index;
ecg_raw = time_index(:,2);
Fs = 5000;
Fn = Fs/2;
time = time_index(:,1)/Fs;

L = numel(ecg_raw);
N = 2^nextpow2(L);
FTs = fft(ecg_raw,N)/L;
Fv = linspace(0, 1, N/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTs(Iv))*2)
grid
title('Fourier Transform')
xlim([0  170])

fcomb = [[0.5 1.0], [45 48 52 55], [45 48 52 55]+50]
mags = [[0 1], [0 1], [0 1]];
dev = [[0.5 0.1], [0.1 0.5], [0.1 0.5]];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

ecg_filt = filtfilt(hh, 1, ecg_raw);

figure
freqz(hh, 1, 2^16, Fs)

figure
subplot(2,1,1)
plot(time, ecg_raw)
grid
title('Original Signal')

subplot(2,1,2)
plot(time, ecg_filt)
grid
title('Bandstop-Filtered Signal')

figure
subplot(2,1,1)
plot(time, ecg_raw)
grid
xlim([20 22])
title('Original Signal (Detail)')

subplot(2,1,2)
plot(time, ecg_filt)
grid
xlim([20 22])
title('Bandstop-Filtered Signal (Detail)')

The online application does not allow me to load .mat files (I mentioned this, so perhaps it will be fixed shortly), so I did all this offline, producing for the last plot —

That eliminates the baselline variation and all the noise possible with a frequency-selective filter. The quality of the recoved trace is (in my clinical perspective) the best that can be achieved. It is important that none of the details (S-T segment, isoelectric reference of the period between the end of the P-deflection and the beginning of the QRS complex, and others) be preserved. They are, here.

EDIT — (28 Apr 2021 at 18:34)

Added the previously omitted ‘Fn’ assignment.

.

[Christian Otte]

Thanks a lot. I works perfectly, also with another dataset. Besides, I can't see the variable Fn. But I guess it's based on Nyquist theory and therefore: Fn = Fs/2?

Also, figure 4, which look like it has no plot in it?

My only challenge now is that you delivered coding material at a whole other level than the university teachers prepared me fore :) I only have one more question. Is it possible that I can persuade you commenting each line of the code so I have a chance to do some research regarding how this filter works?

[Star Strider]

I was in the middle of responding when Windows crashed. Again. (I really hate the current version of Windows!)

Anyway, since it is not possible to load .mat files in the online application that I now routinely use to run my code in my Answers and Comments, I downloaded the .mat file and ran it on this computer, offline. I initially forgot to copy the ‘Fn’ assignment, however it is now added where it should be.

The last ‘(Detail)’ figure is specifically for the provided data. Note that it has a xlim call at the end, and this is likely the reason a plot using a different data set is empty, since it may not contain the same x-axis values. (The last figure is not necessary for the code. It simply demonstrates that the FIR filter works even better than I though it would!)

The filter itself is straightforward. A detailed discussion of it is in the kaiserord documentation. Beyond that, it is simply a standard FIR filter with the advantage that it combines a highpass filter (eliminiating the baseline variation) and two stopband filters to eliminate the mains frequency interference and the first harmonic (the only one that appears to add significant noise). There are several excellent signal processing texts available that go into significant detail on all types of filters, of course including FIR filters. They discuss them in much more detail (and with significantly greater authority) than I could, so I refer you to them. Many are listed in the references section at the end of each documentation page.

Answer 2

[Mathieu NOE]

hello

dear all, an alternative using a 5th order IIR filter to remove the usual 50 Hz and harmonics mains hum

LD = load('ECG_sample.mat');
time_index = LD.time_index;
ecg_raw = time_index(:,2);
Fs = 5000;
Fn = Fs/2;
time = time_index(:,1)/Fs;
L = numel(ecg_raw);
N = 2^nextpow2(L);
FTs = fft(ecg_raw,N)/L;
Fv = linspace(0, 1, N/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTs(Iv))*2)
grid
title('Fourier Transform')
xlim([0  170])

%% two notch filters in series to remove 50 and 100 Hz hum
fc = 50;       % notch freq for 50 Hz tone
Q = 7;         % adjust Q factor for wider (low Q) / narrower (high Q) notch ; f = fc the filter has gain = 0
[B1,A1] = mynotchfilter(fc,Q, Fs);

fc = 2*50;       % notch freq for 100 Hz tone
Q = 7;         % adjust Q factor for wider (low Q) / narrower (high Q) notch ; f = fc the filter has gain = 0
[B2,A2] = mynotchfilter(fc,Q, Fs);

% two filter is series
B = conv(B1,B2);
A = conv(A1,A2);

ecg_filt = filtfilt(B, A, ecg_raw);

figure
freqz(B, A, 2^16, Fs)
figure
subplot(2,1,1)
plot(time, ecg_raw)
grid
title('Original Signal')
subplot(2,1,2)
plot(time, ecg_filt)
grid
title('Bandstop-Filtered Signal')
figure
subplot(2,1,1)
plot(time, ecg_raw)
grid
xlim([20 22])
title('Original Signal (Detail)')
subplot(2,1,2)
plot(time, ecg_filt)
grid
xlim([20 22])
title('Bandstop-Filtered Signal (Detail)')

function [B,A] = mynotchfilter(fc,Q, Fs)
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)

% fc = 50;       % notch freq
wc = 2*pi*fc;
% Q = 7;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0

w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);

            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
B=[b0 b1 b2];
A=[a0 a1 a2];
end