-
Jean-Noël Bazin authoredJean-Noël Bazin authored
pan_tompkin.m 16.73 KiB
function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs,gr)
%% function [qrs_amp_raw,qrs_i_raw,delay]=pan_tompkin(ecg,fs)
% Complete implementation of Pan-Tompkins algorithm
%% Inputs
% ecg : raw ecg vector signal 1d signal
% fs : sampling frequency e.g. 200Hz, 400Hz and etc
% gr : flag to plot or not plot (set it 1 to have a plot or set it zero not
% to see any plots
%% Outputs
% qrs_amp_raw : amplitude of R waves amplitudes
% qrs_i_raw : index of R waves
% delay : number of samples which the signal is delayed due to the
% filtering
%% Method
% See Ref and supporting documents on researchgate.
% https://www.researchgate.net/publication/313673153_Matlab_Implementation_of_Pan_Tompkins_ECG_QRS_detector
%% References :
%[1] Sedghamiz. H, "Matlab Implementation of Pan Tompkins ECG QRS
%detector.",2014. (See researchgate)
%[2] PAN.J, TOMPKINS. W.J,"A Real-Time QRS Detection Algorithm" IEEE
%TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-32, NO. 3, MARCH 1985.
%% ============== Licensce ========================================== %%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
% FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
% OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
% TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% Author :
% Hooman Sedghamiz, Feb, 2018
% MSc. Biomedical Engineering, Linkoping University
% Email : Hooman.sedghamiz@gmail.com
%% ============ Update History ================== %%
% Feb 2018 :
% 1- Cleaned up the code and added more comments
% 2- Added to BioSigKit Toolbox
%% ================= Now Part of BioSigKit ==================== %%
if ~isvector(ecg)
error('ecg must be a row or column vector');
end
if nargin < 3
gr = 1; % on default the function always plots
end
ecg = ecg(:); % vectorize
%% ======================= Initialize =============================== %
delay = 0;
skip = 0; % becomes one when a T wave is detected
m_selected_RR = 0;
mean_RR = 0;
ser_back = 0;
ax = zeros(1,6);
%% ============ Noise cancelation(Filtering)( 5-15 Hz) =============== %%
if fs == 200
% ------------------ remove the mean of Signal -----------------------%
ecg = ecg - mean(ecg);
%% ==== Low Pass Filter H(z) = ((1 - z^(-6))^2)/(1 - z^(-1))^2 ==== %%
%%It has come to my attention the original filter doesnt achieve 12 Hz
% b = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
% a = [1 -2 1];
% ecg_l = filter(b,a,ecg);
% delay = 6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Wn = 12*2/fs;
N = 3; % order of 3 less processing
[a,b] = butter(N,Wn,'low'); % bandpass filtering ecg_l = filtfilt(a,b,ecg);
ecg_l = ecg_l/ max(abs(ecg_l));
%% ======================= start figure ============================= %%
if gr
figure;
ax(1) = subplot(321);plot(ecg);axis tight;title('Raw signal');
ax(2)=subplot(322);plot(ecg_l);axis tight;title('Low pass filtered');
end
%% ==== High Pass filter H(z) = (-1+32z^(-16)+z^(-32))/(1+z^(-1)) ==== %%
%%It has come to my attention the original filter doesn achieve 5 Hz
% b = zeros(1,33);
% b(1) = -1; b(17) = 32; b(33) = 1;
% a = [1 1];
% ecg_h = filter(b,a,ecg_l); % Without Delay
% delay = delay + 16;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Wn = 5*2/fs;
N = 3; % order of 3 less processing
[a,b] = butter(N,Wn,'high'); % bandpass filtering
ecg_h = filtfilt(a,b,ecg_l);
ecg_h = ecg_h/ max(abs(ecg_h));
if gr
ax(3)=subplot(323);plot(ecg_h);axis tight;title('High Pass Filtered');
end
else
%% bandpass filter for Noise cancelation of other sampling frequencies(Filtering)
f1=5; % cuttoff low frequency to get rid of baseline wander
f2=15; % cuttoff frequency to discard high frequency noise
Wn=[f1 f2]*2/fs; % cutt off based on fs
N = 3; % order of 3 less processing
[a,b] = butter(N,Wn); % bandpass filtering
ecg_h = filtfilt(a,b,ecg);
ecg_h = ecg_h/ max( abs(ecg_h));
if gr
ax(1) = subplot(3,2,[1 2]);plot(ecg);axis tight;title('Raw Signal');
ax(3)=subplot(323);plot(ecg_h);axis tight;title('Band Pass Filtered');
end
end
%% ==================== derivative filter ========================== %%
% ------ H(z) = (1/8T)(-z^(-2) - 2z^(-1) + 2z + z^(2)) --------- %
if fs ~= 200
int_c = (5-1)/(fs*1/40);
b = interp1(1:5,[1 2 0 -2 -1].*(1/8)*fs,1:int_c:5);
else
b = [1 2 0 -2 -1].*(1/8)*fs;
end
ecg_d = filtfilt(b,1,ecg_h);
ecg_d = ecg_d/max(ecg_d);
if gr
ax(4)=subplot(324);plot(ecg_d);
axis tight;
title('Filtered with the derivative filter');
end
%% ========== Squaring nonlinearly enhance the dominant peaks ========== %%
ecg_s = ecg_d.^2;
if gr
ax(5)=subplot(325);
plot(ecg_s);
axis tight;
title('Squared');
end
%% ============ Moving average ================== %%
%-------Y(nt) = (1/N)[x(nT-(N - 1)T)+ x(nT - (N - 2)T)+...+x(nT)]---------%
ecg_m = conv(ecg_s ,ones(1 ,round(0.150*fs))/round(0.150*fs));
delay = delay + round(0.150*fs)/2;
if gr
ax(6)=subplot(326);plot(ecg_m);
axis tight;
title('Averaged with 30 samples length,Black noise,Green Adaptive Threshold,RED Sig Level,Red circles QRS adaptive threshold');
axis tight; end
%% ===================== Fiducial Marks ============================== %%
% Note : a minimum distance of 40 samples is considered between each R wave
% since in physiological point of view no RR wave can occur in less than
% 200 msec distance
[pks,locs] = findpeaks(ecg_m,'MINPEAKDISTANCE',round(0.2*fs));
%% =================== Initialize Some Other Parameters =============== %%
LLp = length(pks);
% ---------------- Stores QRS wrt Sig and Filtered Sig ------------------%
qrs_c = zeros(1,LLp); % amplitude of R
qrs_i = zeros(1,LLp); % index
qrs_i_raw = zeros(1,LLp); % amplitude of R
qrs_amp_raw= zeros(1,LLp); % Index
% ------------------- Noise Buffers ---------------------------------%
nois_c = zeros(1,LLp);
nois_i = zeros(1,LLp);
% ------------------- Buffers for Signal and Noise ----------------- %
SIGL_buf = zeros(1,LLp);
NOISL_buf = zeros(1,LLp);
SIGL_buf1 = zeros(1,LLp);
NOISL_buf1 = zeros(1,LLp);
THRS_buf1 = zeros(1,LLp);
THRS_buf = zeros(1,LLp);
%% initialize the training phase (2 seconds of the signal) to determine the THR_SIG and THR_NOISE
THR_SIG = max(ecg_m(1:2*fs))*1/3; % 0.25 of the max amplitude
THR_NOISE = mean(ecg_m(1:2*fs))*1/2; % 0.5 of the mean signal is considered to be noise
SIG_LEV= THR_SIG;
NOISE_LEV = THR_NOISE;
%% Initialize bandpath filter threshold(2 seconds of the bandpass signal)
THR_SIG1 = max(ecg_h(1:2*fs))*1/3; % 0.25 of the max amplitude
THR_NOISE1 = mean(ecg_h(1:2*fs))*1/2;
SIG_LEV1 = THR_SIG1; % Signal level in Bandpassed filter
NOISE_LEV1 = THR_NOISE1; % Noise level in Bandpassed filter
%% ============ Thresholding and desicion rule ============= %%
Beat_C = 0; % Raw Beats
Beat_C1 = 0; % Filtered Beats
Noise_Count = 0; % Noise Counter
for i = 1 : LLp
%% ===== locate the corresponding peak in the filtered signal === %%
if locs(i)-round(0.150*fs)>= 1 && locs(i)<= length(ecg_h)
[y_i,x_i] = max(ecg_h(locs(i)-round(0.150*fs):locs(i)));
else
if i == 1
[y_i,x_i] = max(ecg_h(1:locs(i)));
ser_back = 1;
elseif locs(i)>= length(ecg_h)
[y_i,x_i] = max(ecg_h(locs(i)-round(0.150*fs):end));
end
end
%% ================= update the heart_rate ==================== %%
if Beat_C >= 9
diffRR = diff(qrs_i(Beat_C-8:Beat_C)); % calculate RR interval
mean_RR = mean(diffRR); % calculate the mean of 8 previous R waves interval
comp =qrs_i(Beat_C)-qrs_i(Beat_C-1); % latest RR
if comp <= 0.92*mean_RR || comp >= 1.16*mean_RR
% ------ lower down thresholds to detect better in MVI -------- %
THR_SIG = 0.5*(THR_SIG);
THR_SIG1 = 0.5*(THR_SIG1);
else
m_selected_RR = mean_RR; % The latest regular beats mean
end
end
%% == calculate the mean last 8 R waves to ensure that QRS is not ==== %%
if m_selected_RR
test_m = m_selected_RR; %if the regular RR availabe use it
elseif mean_RR && m_selected_RR == 0
test_m = mean_RR; else
test_m = 0;
end
if test_m
if (locs(i) - qrs_i(Beat_C)) >= round(1.66*test_m) % it shows a QRS is missed
[pks_temp,locs_temp] = max(ecg_m(qrs_i(Beat_C)+ round(0.200*fs):locs(i)-round(0.200*fs))); % search back and locate the max in this interval
locs_temp = qrs_i(Beat_C)+ round(0.200*fs) + locs_temp -1; % location
if pks_temp > THR_NOISE
Beat_C = Beat_C + 1;
qrs_c(Beat_C) = pks_temp;
qrs_i(Beat_C) = locs_temp;
% ------------- Locate in Filtered Sig ------------- %
if locs_temp <= length(ecg_h)
[y_i_t,x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):locs_temp));
else
[y_i_t,x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):end));
end
% ----------- Band pass Sig Threshold ------------------%
if y_i_t > THR_NOISE1
Beat_C1 = Beat_C1 + 1;
qrs_i_raw(Beat_C1) = locs_temp-round(0.150*fs)+ (x_i_t - 1);% save index of bandpass
qrs_amp_raw(Beat_C1) = y_i_t; % save amplitude of bandpass
SIG_LEV1 = 0.25*y_i_t + 0.75*SIG_LEV1; % when found with the second thres
end
not_nois = 1;
SIG_LEV = 0.25*pks_temp + 0.75*SIG_LEV ; % when found with the second threshold
end
else
not_nois = 0;
end
end
%% =================== find noise and QRS peaks ================== %%
if pks(i) >= THR_SIG
% ------ if No QRS in 360ms of the previous QRS See if T wave ------%
if Beat_C >= 3
if (locs(i)-qrs_i(Beat_C)) <= round(0.3600*fs)
Slope1 = mean(diff(ecg_m(locs(i)-round(0.075*fs):locs(i)))); % mean slope of the waveform at that position
Slope2 = mean(diff(ecg_m(qrs_i(Beat_C)-round(0.075*fs):qrs_i(Beat_C)))); % mean slope of previous R wave
if abs(Slope1) <= abs(0.5*(Slope2)) % slope less then 0.5 of previous R
Noise_Count = Noise_Count + 1;
nois_c(Noise_Count) = pks(i);
nois_i(Noise_Count) = locs(i);
skip = 1; % T wave identification
% ----- adjust noise levels ------ %
NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
else
skip = 0;
end
end
end
%---------- skip is 1 when a T wave is detected -------------- %
if skip == 0
Beat_C = Beat_C + 1;
qrs_c(Beat_C) = pks(i);
qrs_i(Beat_C) = locs(i);
%--------------- bandpass filter check threshold --------------- %
if y_i >= THR_SIG1
Beat_C1 = Beat_C1 + 1;
if ser_back
qrs_i_raw(Beat_C1) = x_i; % save index of bandpass
else
qrs_i_raw(Beat_C1)= locs(i)-round(0.150*fs)+ (x_i - 1); % save index of bandpass
end qrs_amp_raw(Beat_C1) = y_i; % save amplitude of bandpass
SIG_LEV1 = 0.125*y_i + 0.875*SIG_LEV1; % adjust threshold for bandpass filtered sig
end
SIG_LEV = 0.125*pks(i) + 0.875*SIG_LEV ; % adjust Signal level
end
elseif (THR_NOISE <= pks(i)) && (pks(i) < THR_SIG)
NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1; % adjust Noise level in filtered sig
NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV; % adjust Noise level in MVI
elseif pks(i) < THR_NOISE
Noise_Count = Noise_Count + 1;
nois_c(Noise_Count) = pks(i);
nois_i(Noise_Count) = locs(i);
NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1; % noise level in filtered signal
NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV; % adjust Noise level in MVI
end
%% ================== adjust the threshold with SNR ============= %%
if NOISE_LEV ~= 0 || SIG_LEV ~= 0
THR_SIG = NOISE_LEV + 0.25*(abs(SIG_LEV - NOISE_LEV));
THR_NOISE = 0.5*(THR_SIG);
end
%------ adjust the threshold with SNR for bandpassed signal -------- %
if NOISE_LEV1 ~= 0 || SIG_LEV1 ~= 0
THR_SIG1 = NOISE_LEV1 + 0.25*(abs(SIG_LEV1 - NOISE_LEV1));
THR_NOISE1 = 0.5*(THR_SIG1);
end
%--------- take a track of thresholds of smoothed signal -------------%
SIGL_buf(i) = SIG_LEV;
NOISL_buf(i) = NOISE_LEV;
THRS_buf(i) = THR_SIG;
%-------- take a track of thresholds of filtered signal ----------- %
SIGL_buf1(i) = SIG_LEV1;
NOISL_buf1(i) = NOISE_LEV1;
THRS_buf1(i) = THR_SIG1;
% ----------------------- reset parameters -------------------------- %
skip = 0;
not_nois = 0;
ser_back = 0;
end
%% ======================= Adjust Lengths ============================ %%
qrs_i_raw = qrs_i_raw(1:Beat_C1);
qrs_amp_raw = qrs_amp_raw(1:Beat_C1);
qrs_c = qrs_c(1:Beat_C);
qrs_i = qrs_i(1:Beat_C);
%% ======================= Plottings ================================= %%
if gr
hold on,scatter(qrs_i,qrs_c,'m');
hold on,plot(locs,NOISL_buf,'--k','LineWidth',2);
hold on,plot(locs,SIGL_buf,'--r','LineWidth',2);
hold on,plot(locs,THRS_buf,'--g','LineWidth',2);
if any(ax)
ax(~ax) = [];
linkaxes(ax,'x');
zoom on;
end
end
%% ================== overlay on the signals ========================= %%
if gr
figure;
az(1)=subplot(311);
plot(ecg_h);
title('QRS on Filtered Signal');
axis tight;
hold on,scatter(qrs_i_raw,qrs_amp_raw,'m');
hold on,plot(locs,NOISL_buf1,'LineWidth',2,'Linestyle','--','color','k');
hold on,plot(locs,SIGL_buf1,'LineWidth',2,'Linestyle','-.','color','r'); hold on,plot(locs,THRS_buf1,'LineWidth',2,'Linestyle','-.','color','g');
az(2)=subplot(312);plot(ecg_m);
title('QRS on MVI signal and Noise level(black),Signal Level (red) and Adaptive Threshold(green)');axis tight;
hold on,scatter(qrs_i,qrs_c,'m');
hold on,plot(locs,NOISL_buf,'LineWidth',2,'Linestyle','--','color','k');
hold on,plot(locs,SIGL_buf,'LineWidth',2,'Linestyle','-.','color','r');
hold on,plot(locs,THRS_buf,'LineWidth',2,'Linestyle','-.','color','g');
az(3)=subplot(313);
plot(ecg-mean(ecg));
title('Pulse train of the found QRS on ECG signal');
axis tight;
line(repmat(qrs_i_raw,[2 1]),...
repmat([min(ecg-mean(ecg))/2; max(ecg-mean(ecg))/2],size(qrs_i_raw)),...
'LineWidth',2.5,'LineStyle','-.','Color','r');
linkaxes(az,'x');
zoom on;
end
end