Butterworth filter for ecg signal in healthy people

Question

Christian Otte 2021-4-27

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/814710-butterworth-filter-for-ecg-signal-in-healthy-people

回答： Mathieu NOE 2021-4-29

Skærmbillede 2021-04-27 kl. 16.15.47.png

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?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2021-4-28

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 2021-4-28

ECG_sample.mat

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 2021-4-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/814710-butterworth-filter-for-ecg-signal-in-healthy-people#answer_686215

在 MATLAB Online 中打开

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).

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2021-4-28

编辑：Star Strider 2021-4-28

在 MATLAB Online 中打开

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 2021-4-28

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 2021-4-28

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 2021-4-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/814710-butterworth-filter-for-ecg-signal-in-healthy-people#answer_687630

在 MATLAB Online 中打开

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