How to use wavelet filters for peak finding ECG data

I have a 170 datasets of ECG data that im looking to process for RR peaks and average Waveform shaping. Each dataset is approximatley 2 minutes long, and it all is very noisy and varied. I am interested in finding the peaks but using "findpeaks()" alone is proving insufficient since there are significant noise artificats and human subject variability. I've heard others have used wavelet filters, but im a little stuck with how to use them. Is there an ECG example I can follow to learn which functions are useful for this type of denoising?

As an example when I plot a specific case I see several siginifant large amplitude noise artifacts. The usable ECG data is between ~70 to 75 sec and from 90 seconds to about 120 seconds. Several good peaks have been missed, and several noise peaks have been incorrectly labeled. Furthermore the findpeaks function requires hard coded .6 and .3 values that I wouldnt know apriori. My hope is that wavelet filtering can reduce the number of false positive and false negatives.

Any advice would be appreciated.

ss = 6; ee = 10; 
t = dataset{ss,ee}(:,1);        %Time Vector
v = dataset{ss,ee}(:,2);        %Voltage Vector
[rpk,tpk] = findpeaks(v,t,'MinPeakDistance',.6,'MinPeakHeight',.3);

plot(t,v)
hold on; grid on;
xlabel('Time (sec)');
ylabel('Amplitude (mV)');

[Mathieu NOE]

hello

some suggestions :

Wavelet Denoising - MATLAB & Simulink - MathWorks France

Wavelets Based Denoising - File Exchange - MATLAB Central (mathworks.com)

[Mathieu NOE]

Wavelet signal denoising - MATLAB wdenoise - MathWorks France

Answer 1

[Star Strider]

There are several options depending on what the noise is. The best way to determine that is to take the Fourier transform of the signal using fft or pspectrum. If the noise is band-limited, use a lowpass filter to remove it (choosing 'ImpulseResponse','iir' for best results). If it is broadband, either use wavelet denoising or the sgolayfilt function (I usually choose a 3-degree polynomial and then adjust the framelen value to get the result I want).

[Peter Brehm]

Okay,

This seems like a method to filter the out of band noise (High Pass?). I'm unsure how to properly select order and framelen. I am also still struggling with in band noise. The 60Hz noise obviously got cleared up which is good. But it still requires that I have apriori knowledge of peak height in order to correctly identify them. Im also worried about the distortion of the existing waveforms. Even with this method, I still have a high number of false positives and false negatives.

load('ECGdata_sub6_ee10.mat')   % Load in raw data
order = 4;                      % How to pick this?
framelen = 201;                 % How to pick this?
vsg = sgolayfilt(v,order,framelen); %Filter dataset with Sgolay filter

mph = .1;                       % min Peak Height (How do I Pick this?)
mpd = .6;                       % min peak distance (RR interval dependent)
[rpk,tpk] = findpeaks(vsg,t,'MinPeakDistance',mpd,'MinPeakHeight',mph);

[f,P1] = MakeSpectrum(t,v);     %Generate frequency content for plot
[f,Psg] = MakeSpectrum(t,vsg);  %Generate Frequency content for plot

% Create Plots
figure(1)
subplot(2,1,1)
plot(t,v)
hold on; grid on;
plot(t,vsg)
plot(tpk,rpk,'ok')
xlabel('Time (sec)'); ylabel('Amplitude (mV)');
subplot(2,1,2)
plot(t,v)
hold on; grid on;
plot(t,vsg)
plot(tpk,rpk,'ok')
xlabel('Time (sec)'); ylabel('Amplitude (mV)');
xlim([80,90])

figure(2)
plot(f,P1)
hold on; grid on;
plot(f,Psg);
xlabel('Frequency (Hz)');ylabel('Amplitude (mV)');
xlim([0,100])

function [f,P1] = MakeSpectrum(t,X)
NumFSamp = length(X);           %(#) Number of Freq Samples
dt = mean(diff(t));             %(sec/#) Average Time btwn samples
fs = 1/dt;                      %(#/sec) Sampleing Frequency
Y = fft(X);                     %() Fourier Transform
Y(abs(Y)<1e-10)=0;              % Set anything smaller than 1e-10 to 0
P2 = abs(Y/NumFSamp);           % Magnitude of 2 sided spectrum
P1 = 2*P2(1:floor(NumFSamp/2)+1);   % (V) Single Sided Spectrum 
f = fs*(0:(NumFSamp/2))/NumFSamp;
end

[Star Strider]

Unfortunately, I doubt that there is any way to salvage most of this record, although some of it can be recovered —

LD = load('ECGdata_sub6_ee10.mat');
t = LD.t;
v = LD.v;

Fs = 1/mean(diff(t));                                           % Sampling Frequency

figure
plot(t,v)
grid

figure
plot(t,v)
grid
axis([60 80  -1 1])

v = bandstop(v, [58 62], Fs, 'ImpulseResponse','iir');          % Notch Out 60 Hz Line Noise
vf = bandpass(v, [2.5 75], Fs, 'ImpulseResponse','iir');        % Remove Some Baseline Drift & High-Frequency Noise
framelen = 51;                                                  % Experiment To Get The Best Result
vf = sgolayfilt(vf, 3, framelen);                               % Suppress  As Much Broadband Noise As Possible

figure
plot(t,vf)
grid
axis([60 80  -1 1])

figure
plot(t,vf)
grid

[pks1,locs1] = findpeaks(vf, 'MinPeakProminence',0.6);
[pks2,locs2] = findpeaks(vf, 'MinPeakProminence',0.4);
locsd = setdiff(locs2,locs1,'stable');

figure
plot(t,vf)
hold on
% plot(t(locs1),vf(locs1), '+r')
% plot(t(locs2),vf(locs2), '+r')
plot(t(locsd),vf(locsd), '+r')
hold off
grid

This is my best effort. There is still a fair amount of noise, and if I was reading this record, I would send it back to be done again, since there are too many problems with it for it to be clinically meaningful. The ‘framelen’ value must be long enough to suppress as much noise as possible, yet short enough to not hide any meaningful information in the EKG trace. That is not a straightforward decision, and there is no way I know of to automate it using a Fourier transform or other analytic results. Much of the remaining analysis will necessarily have to be done manually. since there is no reliable way to separate the ‘good’ parts of the record from the ‘corrupted’ parts of the record using automated techniques. The two findpeaks calls and the setdiff call are the best I can do with respect to isolaitng the R-deflections.

It is always best to use a right leg common reference (sometimes erroneously termed a ‘ground’) in order to subtract whatever that picks up from the other leads. Also, all leads must have coaxial cables with the shield grounded, and always use a earth-chassis ground as well, being mindful of ground faults with other instrumentation attached to the patient.

I did this offlilne because it is stilll not possible to work with .mat files with the online Run feature.

.

[Peter Brehm]

Thank you so much for this suggestion. I hadnt seen the setdiff() function before. Im sure I will make use of it in the future. For some more context, I am a grad student investigating new types of electrodes. I have 16 sets of electrodes, plus the traditional Adhesive gel electrode. For each electrode type (17), i have about 2 minute of data from 10 human subjects. So in total its 170 datasets. This is too many for me to manually step through so I am trying to write a script to automate the peak finding process.

Some are very "well behaved" like the traditional gel electrode, but others are not and have significant noise/motion artefacts. Im trying to locate the Rpeaks so that I can compare the various electrodes against each other, and determine their skin-electrode impedance parameters. I learned today that a past grad student used Daubechie wavelet filters so I will go investigate those next.

And thanks again!

[Star Strider]

My pleasure!

Apparently I have not actually solved the problem, so I will delete my answer.

[Peter Brehm]

haha, please dont, I accepted the answer! I havent submitted many matlab questions, so i just forgot accept!

[Star Strider]

Thank you!

I usually delete my answers if they don’t solve the problem (or if they aren’t accepted), so I am pleased that this one appears to have helped.

Wavelets can be good at denoising, however they are not always successful in my experience, and for that reason I usually do either frequuency-based filtering or use sgolayfilt for broadband noise.

There are many wavelet families that would be worth experimenting with. For a discrete signal, Haar wavelets could be worth trying as well. See Discrete Multiresolution Analysis for a relevant discussion.

Another option might be independent component analysis using the rica function or relevant File Exchange contributions. I have used ICA for EKG analysis, however I do not have any recent experience with it.