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

6 次查看（过去 30 天）
Stephan Lallinger
Stephan Lallinger 2019-11-4
评论： Daniel M 2019-11-6
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.
(Reading the data from https://www.mathworks.com/matlabcentral/answers/487640-ecg-spectrum-wrong#comment_762871)
  4 个评论
显示 2更早的评论
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 Center 和 File Exchange 中查找有关 Digital Filtering 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by