How to get rid of ringing/overshoot in an ecg signal while using notch filter?

4 vues (au cours des 30 derniers jours)
Hallo all,
i have an ecg signal which contains 50 Hz noise. I want to get rid off the noise.
First of all i do a pre processing step with a highpass and lowpass and. After this i use the notch filter.
Notch filter code:
f_noise = 50; %noise to filter
fs = 500;
f_notch = f_noise/fs; %normalized noise
N =[exp(j*2*pi*f_notch);exp(-j*2*pi*f_notch)]; %place zeros
P = 0.5*N; %place additional poles with same angular as zeros
b = poly(N); a = poly(P); %get filter coeff.
Problem here i got an "overshoot" in the QRS complex (see Picture).
A other member suggested this notch filter to implement:
%Notchfilter by Daniel. M
Fs = 500;
Fn = fs/2; % Nyquist frequency
numHarmonics = 0; % Let's do 50, 100, 150, 200, 250; 0 -> filter only 50 Hz
lineFreq = 50; % Hz
for fq = ((0:numHarmonics)+1) * lineFreq
Fl = fq + [-1, 1]; % notch around Fl. Could try [-2, 2] if too tight
[z,p,k] = butter(1, Fl/Fn, 'stop');
sos = zp2sos(z,p,k);
ecg_signal_005_150_notch = filtfilt(sos, 1, ecg_signal_005_150); % assumes data is [time x ... dimensions]
% overwrites data, and filters sequentially for each notch
end
This solution removes the "overshoot" perfectly but it adds "ringing" (see picture).
Does somebody know a good trade-off between ringing and overshoot?
Best regards
  1 commentaire
Daniel M
Daniel M le 4 Nov 2019
Modifié(e) : Daniel M le 4 Nov 2019
I took another crack at this and couldn't come up with anything satisfactory (within my available time). Can you not just do a band pass filter from 0.5 to 40?

Connectez-vous pour commenter.

Réponse acceptée

Daniel M
Daniel M le 4 Nov 2019
Modifié(e) : Daniel M le 4 Nov 2019
clearvars
clc
close all
fs = 500;
ecg = xlsread('Noisy_50HZ_ECG.xlsx');
ecg = detrend(ecg.');
ecg = sgolayfilt(ecg,0,5);
T = 1/fs;
L = length(ecg);
t = (0:L-1)*T;
filtfreq = [0.5 40];
xlim = [49 51];
%%%% Apply my own filters
[B,A] = butter(4, filtfreq./(fs/2));
f_ecg = filtfilt(B,A,ecg);
%%% Notch filter
data = f_ecg;
Fn = fs/2; % Nyquist frequency
numHarmonics = 3; % Let's do 50, 100, 150, 200
lineFreq = 50; % Hz
for fq = ((0:numHarmonics)+1) * lineFreq
Fl = fq + [-1, 1];
[z,p,k] = butter(6, Fl/Fn, 'stop');
sos = zp2sos(z,p,k);
data = filtfilt(sos, 1, data);
end
figure
plot(t,ecg,'b.-',t,data,'r-')
set(gca,'XLim',[48 52])
figure
plot(t,data)
set(gca,'XLim',[49 51])
I also tried an FIR notch, using the fieldtrip toolbox, different bandpasses, etc. But this is what I think worked best.
  4 commentaires
Stephan Lallinger
Stephan Lallinger le 6 Nov 2019
Modifié(e) : Stephan Lallinger le 6 Nov 2019
Maybe a notch filter is not the best solution for an ecg with the bandwidth of 0.05 - 150 Hz.
I tried the "broad FIR" solution before the notch and it wasn't that good.
I tried also serverals windows like Gaussian, Kaiser and Hann but maybe i could improve here.
Iam also working on adaptive filter solutions with reference noise and without.
Best solutions i got is with your code and this settings:
for fq = ((0:numHarmonics)+1) * lineFreq
Fl = fq + [-0.1, 0.1]; % notch around Fl. Could try [-2, 2] if too tight
[a,b] = butter(1, Fl/Fn, 'stop');
%sos = zp2sos(z,p,k);
ecg_bandpassed = filtfilt(a,b, ecg_bandpassed); % assumes data is [time x ... dimensions]
% overwrites data, and filters sequentially for each notch
end
I needs roughly 5 seconds to settle down but then it is quite good.
To be honest i was suprised because the corners are very sharp and the order is just 1.
It also seems zero phase filtering is much better then using just "filter" or convolution.
Normally i use convolution for filtering i have no experience with "filtfilt" function.
For clearifaction i used not the 50 Hz noisy signal. Instead i used a clear ecg and added the noise via Matlab.
Daniel M
Daniel M le 6 Nov 2019
It makes sense that it would work well with a "perfect" 50 Hz noise signal. But real world applications are not so ideal, as you have learned.
If you eventually find a solution that works well for you, please post here for myself and others to learn from. Good luck!

Connectez-vous pour commenter.

Plus de réponses (0)

Produits


Version

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by