Matlab implementation illustrates the application of Pan-Tompkins Algorithm to a short ECG data. This Matlab code is previously published at http://farukuysal.net/QRSDetection.aspx and can be found many different web page as example code. The original real time decision algorithm substituted with an alternative peak detection algorithm. Alternative peak detection algorithm can not support real time but very effective when using overlapping windows. Second version of this algorithm (not included this module) can detect P and T wave depending on the QRS detection.
% QRS Detection Example % shows the effect of each filter according to Pan-Tompkins algorithm. % Note that, the decision algorithm is different then the mentioned algorithm. % by Faruk UYSAL clear all close all x1 = load('ecg3.dat'); % load the ECG signal from the file fs = 200; % Sampling rate N = length (x1); % Signal length t = [0:N-1]/fs; % time index figure(1) subplot(2,1,1) plot(t,x1) xlabel('second');ylabel('Volts');title('Input ECG Signal') subplot(2,1,2) plot(t(200:600),x1(200:600)) xlabel('second');ylabel('Volts');title('Input ECG Signal 1-3 second') xlim([1 3])
x1 = x1 - mean (x1 ); % cancel DC conponents x1 = x1/ max( abs(x1 )); % normalize to one figure(2) subplot(2,1,1) plot(t,x1) xlabel('second');ylabel('Volts');title(' ECG Signal after cancellation DC drift and normalization') subplot(2,1,2) plot(t(200:600),x1(200:600)) xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second') xlim([1 3])
% LPF (1-z^-6)^2/(1-z^-1)^2 b=[1 0 0 0 0 0 -2 0 0 0 0 0 1]; a=[1 -2 1]; h_LP=filter(b,a,[1 zeros(1,12)]); % transfer function of LPF x2 = conv (x1 ,h_LP); %x2 = x2 (6+[1: N]); %cancle delay x2 = x2/ max( abs(x2 )); % normalize , for convenience . figure(3) subplot(2,1,1) plot([0:length(x2)-1]/fs,x2) xlabel('second');ylabel('Volts');title(' ECG Signal after LPF') xlim([0 max(t)]) subplot(2,1,2) plot(t(200:600),x2(200:600)) xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second') xlim([1 3])
% HPF = Allpass-(Lowpass) = z^-16-[(1-z^-32)/(1-z^-1)] b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]; a = [1 -1]; h_HP=filter(b,a,[1 zeros(1,32)]); % impulse response iof HPF x3 = conv (x2 ,h_HP); %x3 = x3 (16+[1: N]); %cancle delay x3 = x3/ max( abs(x3 )); figure(4) subplot(2,1,1) plot([0:length(x3)-1]/fs,x3) xlabel('second');ylabel('Volts');title(' ECG Signal after HPF') xlim([0 max(t)]) subplot(2,1,2) plot(t(200:600),x3(200:600)) xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second') xlim([1 3])
% Make impulse response h = [-1 -2 0 2 1]/8; % Apply filter x4 = conv (x3 ,h); x4 = x4 (2+[1: N]); x4 = x4/ max( abs(x4 )); figure(5) subplot(2,1,1) plot([0:length(x4)-1]/fs,x4) xlabel('second');ylabel('Volts');title(' ECG Signal after Derivative') subplot(2,1,2) plot(t(200:600),x4(200:600)) xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second') xlim([1 3])
x5 = x4 .^2; x5 = x5/ max( abs(x5 )); figure(6) subplot(2,1,1) plot([0:length(x5)-1]/fs,x5) xlabel('second');ylabel('Volts');title(' ECG Signal Squarting') subplot(2,1,2) plot(t(200:600),x5(200:600)) xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second') xlim([1 3])
% Make impulse response h = ones (1 ,31)/31; Delay = 15; % Delay in samples % Apply filter x6 = conv (x5 ,h); x6 = x6 (15+[1: N]); x6 = x6/ max( abs(x6 )); figure(7) subplot(2,1,1) plot([0:length(x6)-1]/fs,x6) xlabel('second');ylabel('Volts');title(' ECG Signal after Averaging') subplot(2,1,2) plot(t(200:600),x6(200:600)) xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second') xlim([1 3])
figure(7) subplot(2,1,1) max_h = max(x6); thresh = mean (x6 ); poss_reg =(x6>thresh*max_h)'; figure (8) subplot(2,1,1) hold on plot (t(200:600),x1(200:600)/max(x1)) box on xlabel('second');ylabel('Integrated') xlim([1 3]) subplot(2,1,2) plot (t(200:600),x6(200:600)/max(x6)) xlabel('second');ylabel('Integrated') xlim([1 3]) left = find(diff([0 poss_reg])==1); right = find(diff([poss_reg 0])==-1); left=left-(6+16); % cancle delay because of LP and HP right=right-(6+16);% cancle delay because of LP and HP for i=1:length(left) [R_value(i) R_loc(i)] = max( x1(left(i):right(i)) ); R_loc(i) = R_loc(i)-1+left(i); % add offset [Q_value(i) Q_loc(i)] = min( x1(left(i):R_loc(i)) ); Q_loc(i) = Q_loc(i)-1+left(i); % add offset [S_value(i) S_loc(i)] = min( x1(left(i):right(i)) ); S_loc(i) = S_loc(i)-1+left(i); % add offset end % there is no selective wave Q_loc=Q_loc(find(Q_loc~=0)); R_loc=R_loc(find(R_loc~=0)); S_loc=S_loc(find(S_loc~=0)); figure subplot(2,1,1) title('ECG Signal with R points'); plot (t,x1/max(x1) , t(R_loc) ,R_value , 'r^', t(S_loc) ,S_value, '*',t(Q_loc) , Q_value, 'o'); legend('ECG','R','S','Q'); subplot(2,1,2) plot (t,x1/max(x1) , t(R_loc) ,R_value , 'r^', t(S_loc) ,S_value, '*',t(Q_loc) , Q_value, 'o'); xlim([1 3])
Pan, Jiapu; Tompkins, Willis J., "A Real-Time QRS Detection Algorithm," Biomedical Engineering, IEEE Transactions on , vol.BME-32, no.3, pp.230,236, March 1985 doi: 10.1109/TBME.1985.325532