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

조회 수: 16 (최근 30일)
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
Daniel M
Daniel M 2019년 11월 4일
편집: Daniel M 2019년 11월 4일
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?

댓글을 달려면 로그인하십시오.

채택된 답변

Daniel M
Daniel M 2019년 11월 4일
편집: Daniel M 2019년 11월 4일
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
Stephan Lallinger
Stephan Lallinger 2019년 11월 6일
편집: Stephan Lallinger 2019년 11월 6일
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 2019년 11월 6일
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!

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Single-Rate Filters에 대해 자세히 알아보기

제품


릴리스

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by