QRS Detection: Pan Tompkins Algorithm

One of the most straightforward aspects of ECG analysis is to detect QRS complexes, which are the super obvious "peaks" that occur regularly in the signal. Each represents one heartbeat, and therefore extracting them can yield valuable information for downstream analysis. It's also generally the first step of all automated ECG analysis.

The most widely used algorithm for QRS detection is the Pan Tompkins algorithm.

I'm going to walk through the algorithm on the motion artifact data since it's simpler, and then apply it to the segment of arrhythmia data.

In [1]:
using CSV, Plots, Images #, TimeSeries
using StatsBase
plotlyjs()

Plotly javascript loaded.

To load again call 

init_notebook(true)

WARNING: Method definition ==(Base.Nullable{S}, Base.Nullable{T}) in module Base at nullable.jl:238 overwritten in module NullableArrays at C:\Users\sandy\AppData\Local\JuliaPro-0.6.0.1\pkgs-0.6.0.1\v0.6\NullableArrays\src\operators.jl:128.
WARNING: Compat.AsyncCondition is deprecated, use
Out[1]:
Plots.PlotlyJSBackend()
Base.AsyncCondition instead.
  likely near In[1]:1
WARNING:
In [2]:
fs_motion = 500;
sample_motion = CSV.read("test22_00j.csv", header = ["Time", "ECG1", "ECG2", "ECG3", "ECG4"], datarow=1, nullable=false)
plot(sample_motion[:,1],sample_motion[:,4],title = "Sample ECG during Jumping with Motion Artifact",xlabel="Time (sec)")

Compat.AsyncCondition is deprecated, use Base.AsyncCondition instead.
  likely near In[1]:1
Out[2]:
In [3]:
fs = fs_motion
T = 1/fs
x = sample_motion[:,4];

#Mean shift the signal
x .= x .- mean(x);

Low-pass Filter

In [4]:
x_LP = zeros(size(x));
x_LP[1:12] = x[1:12];
for n = 13:length(x_LP)
    x_LP[n] = 2x_LP[n-1] - x_LP[n-2] + x[n] - 2x[n-6] + x[n-12]
end
plot(sample_motion[:,1],x_LP,title = "Sample ECG after Low-Pass Filter",xlabel="Time (sec)")
Out[4]:

High-pass Filter

In [5]:
x_HP = zeros(size(x))
for n = 33:length(x_HP)
    x_HP[n] = 32x_LP[n-16] - (x_HP[n-1] + x_LP[n] - x_LP[n-32])
end
plot(sample_motion[:,1],x_HP,title = "Sample ECG after Band-Pass Filter",xlabel="Time (sec)")
Out[5]:

Derivative

In [6]:
x_D = zeros(size(x))
for n = 3:length(x_D)-2
    x_D[n] = (1/(8*T)) * (-x_HP[n-2] - 2x_HP[n-1] + 2x_HP[n+1] + x_HP[n+2])
end
plot(sample_motion[:,1],x_D,title = "Sample ECG after Derivative",xlabel="Time (sec)")
Out[6]:

Squaring

In [7]:
x_S = x_D.^2
plot(sample_motion[:,1],x_S,title = "Sample ECG after Derivative Squaring",xlabel="Time (sec)")
Out[7]:

Moving-Window Integration (Moving Average Filter)

We want the width of the window to be such that it encompasses an entire QRS complex, but not too much more. Normally, QRS complexes are about 100 ms wide, but in the case of some pathologies, can be up to 150 ms. We set the window to be 150 ms wide.

In [8]:
w = Int(fld(0.150,T))
x_MAF = zeros(length(x)-w+1,)
for n = 1:length(x_MAF)
    x_MAF[n] = mean(x_S[n:n+w-1])
end
x_MAF = [x_MAF[1]*ones(w-1,1); x_MAF]
plot(sample_motion[:,1],x_MAF,title = "Sample ECG after Derivative Squaring and Moving Average Filter",xlabel="Time (sec)")
Out[8]:

Now superimposing that with the delayed original signal (total delay = 24 samples from first three steps = 6+16+2) (Amplitude adjusted to be visible on the plot)

In [9]:
x_delay = x
x_delay[1:24] = x[1]
x_delay[25:end] = x[1:end-24]
plot!(sample_motion[:,1],100000000x_delay)
Out[9]:

The QRS complexes correspond to the rising edge of each of the "bumps" in moving average filtered signal. After shifting the original signal by the total delay (above), you can see that the QRS complexes line up with the rising edges.

Fiducial Marks

Find local maxima and make sure they are at least 200 ms apart using the existing Julia function for finding local maxima.

In [10]:
locs_cart = findlocalmaxima(x_MAF)
locs_test = [idx[1] for idx in locs_cart]
pks_test = x_MAF[locs_test];
In [11]:
plot(sample_motion[:,1],x_MAF,title = "MAF Signal",xlabel="Time (sec)")
plot!(sample_motion[locs_test,1],pks_test,line = :scatter)
Out[11]:

Clearly wayyy too many peaks... implement the findpeaks function from MATLAB that allows for specification of more parameters to define peak such as min peak distance.

In [12]:
include("findpeaks.jl")
Out[12]:
findpeaks (generic function with 1 method)
In [13]:
#Adjust filtered signal
x_HP .= x_HP .- sample_motion[:,1].*((x_HP[end]-x_HP[1])/(sample_motion[end,1]-sample_motion[1,1]));
x_HP .= x_HP .- mean(x_HP);
In [14]:
pks,locs = findpeaks(x_MAF,[],Int(0.2/T),[],[])
plot(sample_motion[:,1],x_MAF,title = "MAF Signal with Peaks Identified",xlabel="Time (sec)")
plot!(sample_motion[locs,1],pks,line = :scatter)
Out[14]:

Yay! Success! Though the timepoints are not actually perfectly at the beginning of the plateaus.... (Zoomed in and checked)

Use the first 2 seconds of signal to determine thresholds for signal and noise for the INTEGRATION WAVEFORM

These just set the initial thresholds for signal and noise that will then be adjusted as the peaks are found.

In [15]:
THR_SIG = maximum(x_MAF[1:2*fs])*0.25; # 0.25 of the max amplitude 
THR_NOISE = mean(x_MAF[1:2*fs])*0.5; # 0.5 of the mean signal is considered to be noise
SPKI = THR_SIG;
NPKI = THR_NOISE;

Use the first 2 seconds of signal to determine thresholds for signal and noise for the FILTERED WAVEFORM

In [16]:
THR_SIG1 = maximum(x_HP[1:2*fs])*0.25; # 0.25 of the max amplitude 
THR_NOISE1 = mean(x_HP[1:2*fs])*0.5; 
SPKF = THR_SIG1; # Signal level in Bandpassed filter
NPKF = THR_NOISE1; # Noise level in Bandpassed filter

Thresholding and Online Decision Rule

This just goes through and checks each peak against the previous 8 to see if it should still be considered a peak or not.

In [17]:
# Initialize lots of stuff
qrs_c =zeros(0) #amplitude of R
qrs_i = zeros(Int64,0) #index
y_i = zeros(0)
x_i = zeros(Int64,0)
SIG_LEV = 0; 
nois_c =zeros(0)
nois_i = zeros(Int64,0);
Tskip = 0; # becomes one when a T wave is detected
not_nois = 0; # it is not noise when not_nois = 1
selected_RR = zeros(0) # Selected RR intervals
m_selected_RR = 0;
mean_RR = 0;
qrs_i_raw = zeros(Int64,0);
qrs_amp_raw = zeros(0)
ser_back = 0; 
test_m = 0;
SIGL_buf = zeros(0)
NOISL_buf = zeros(0)
THRS_buf = zeros(0)
SIGL_buf1 = zeros(0)
NOISL_buf1 = zeros(0)
THRS_buf1 = zeros(0)
ax = zeros(1,6);
locs_temp = zeros(Int64,0)
x_i_t = zeros(Int64,0);
In [18]:
for i = 1:length(pks)
    
    #Locate the corresponding peak in the filtered signal 
    if (locs[i]-round(0.150*fs) >= 1) && (locs[i] <= length(x_HP))
        y_i = maximum(x_HP[Int(locs[i]-round(0.150*fs)):locs[i]]) #Since I pad to make the same length, should line up exactly
        x_i = indmax(x_HP[Int(locs[i]-round(0.150*fs)):locs[i]])
    else
        if i == 1
            y_i = maximum(x_HP[1:locs[i]])
            x_i = indmax(x_HP[1:locs[i]])
            ser_back = 1
        elseif locs[i] >= length(x_HP)
            y_i = maximum(x_HP[Int(locs[i]-round(0.150*fs)):length(x_HP)])
            x_i = indmax(x_HP[Int(locs[i]-round(0.150*fs)):length(x_HP)])
        end
    end
    
    
    #Update the heart_rate (last 8 and selected last 8)
    if length(qrs_c) > 8
        
        diffRR = diff(qrs_i[length(qrs_i)-8:length(qrs_i)]); #calculate length of last 8 RR intervals
        mean_RR = mean(diffRR); #calculate the mean of last 8 RR intervals
        comp = qrs_i[length(qrs_i)]-qrs_i[length(qrs_i)-1]; #most recent RR
        
        if comp <= 0.92*mean_RR || comp >= 1.16*mean_RR
            # lower down thresholds to detect better in Integrated Waveform
            THR_SIG *= 0.5;
            #THR_NOISE = 0.5*(THR_SIG); 
                
            # lower down thresholds to detect better in Filtered Waveform 
            THR_SIG1 *= 0.5;
            #THR_NOISE1 = 0.5*(THR_SIG1);          
        else
            m_selected_RR = mean_RR; #Shouldn't this be if ALL 8 are within those bounds, not just the last one?
        end 
    end
    
    #Calculate the mean of the last 8 R waves to make sure that QRS is not missing (If no R detected, trigger a search back)
    if m_selected_RR > 0
        test_m = m_selected_RR; 
    elseif mean_RR > 0 && m_selected_RR == 0 #else if the regular RR available use it   
        test_m = mean_RR;   
    else
        test_m = 0;
    end
        
    if test_m > 0
        if (locs[i] - qrs_i[length(qrs_i)]) >= round(1.66*test_m) #A QRS is missed 
            pks_temp = maximum(x_MAF[Int(qrs_i[length(qrs_i)] + round(0.200*fs)):Int(locs[i] - round(0.200*fs))]) #search back and locate the max in this interval
            locs_temp = indmax(x_MAF[Int(qrs_i[length(qrs_i)] + round(0.200*fs)):Int(locs[i] - round(0.200*fs))])
            locs_temp = Int(qrs_i[length(qrs_i)] + round(0.200*fs) + locs_temp - 1); #Compute actual index of location in overall signal 
             
            if pks_temp > THR_NOISE #As long as it passes the noise barrier
                qrs_c = append!(qrs_c, pks_temp)
                qrs_i = append!(qrs_i, locs_temp);
              
                #Find the location in filtered sig
                if locs_temp <= length(x_HP)
                    y_i_t = maximum(x_HP[Int(locs_temp-round(0.150*fs)):locs_temp])
                    x_i_t = indmax(x_HP[Int(locs_temp-round(0.150*fs)):locs_temp])
                else
                    y_i_t = maximum(x_HP[Int(locs_temp-round(0.150*fs)):length(x_HP)])
                    x_i_t = indmax(x_HP[Int(locs_temp-round(0.150*fs)):length(x_HP)])
                end

                #Take care of bandpass signal threshold
                if y_i_t > THR_NOISE1  
                    qrs_i_raw = append!(qrs_i_raw, Int(locs_temp-round(0.150*fs) + (x_i_t - 1))) #Save actual index on filtered signal
                    qrs_amp_raw = append!(qrs_amp_raw, y_i_t) #Save amplitude on filtered signal 
                    SPKF = 0.25*y_i_t + 0.75*SPKF; #From second thresh - shouldn't it be overall peak for the first part
                end
                not_nois = 1;
                SPKI = 0.25*pks_temp + 0.75*SPKI #When found with the second threshold 
            end
        else
              not_nois = 0;
        end
    end
    
    
    #Find noise and QRS peaks
    if pks[i] >= THR_SIG
        
        #If current QRS candidate occurs within 360 ms of the previous QRS,the algorithm determines if it's T wave or QRS
        if length(qrs_c) >= 3 #More than three previous ones detected already
            if (locs[i]-qrs_i[length(qrs_i)]) <= round(0.360*fs) #Within 360ms of last peak 
                Slope1 = mean(diff(x_MAF[Int(locs[i]-round(0.075*fs)):locs[i]])); #Mean slope of the waveform at that position based on 75 ms period before
                Slope2 = mean(diff(x_MAF[Int(qrs_i[length(qrs_i)]-round(0.075*fs)):qrs_i[length(qrs_i)]])); #Mean slope of previous peak
                
                if abs(Slope1) <= abs(0.5*(Slope2))  #Slope less then 0.5 of previous peak
                    nois_c = append!(nois_c, pks[i]); #Add it to noise
                    nois_i = append!(nois_i, locs[i]);
                    Tskip = 1; #T wave identification FOR LATER - KEEP TRACK OF LOCATIONS OF T WAVES AS WELL
                    
                    #Adjust noise level in both filtered and integrated signal
                    NPKF = 0.125*y_i + 0.875*NPKF;
                    NPKI = 0.125*pks[i] + 0.875*NPKI; 
                else
                    Tskip = 0;
                end
            end
        end
        
        #If not a T wave
        if Tskip == 0  #Tskip is 1 when a T wave is detected - so NOT a T wave    
            qrs_c = append!(qrs_c, pks[i]) #Save the peak
            qrs_i = append!(qrs_i, locs[i])
        
            #Filtered signal - check threshold
            if y_i >= THR_SIG1
                if ser_back == 1 #Peak at the very beginning
                   qrs_i_raw = append!(qrs_i_raw, x_i) #Save index on filtered signal 
                else
                   qrs_i_raw = append!(qrs_i_raw, Int(locs[i]-round(0.150*fs) + (x_i - 1))) #Save index on filtered signal 
                end
                qrs_amp_raw = append!(qrs_amp_raw, y_i) #Save amplitude on filtered signal 
                SPKF = 0.125*y_i + 0.875*SPKF #Adjust threshold for bandpass filtered sig
            end
         
            #Adjust Signal level
            SPKI = 0.125*pks[i] + 0.875*SPKI 
        end
        
    
    elseif (THR_NOISE <= pks[i]) && (pks[i] < THR_SIG) #Higher than noise, but not high enough for peak
        #Adjust Noise level in filtered signal
        NPKF = 0.125*y_i + 0.875*NPKF
        #Adjust Noise level in integrated signal
        NPKI = 0.125*pks[i] + 0.875*NPKI 
        
      
    elseif pks[i] < THR_NOISE #Under noise threshold
        nois_c = append!(nois_c, pks[i]); #Add it to noise
        nois_i = append!(nois_i, locs[i]);
        
        #Adjust noise level in filtered signal
        NPKF = 0.125*y_i + 0.875*NPKF;    
        #Adjust noise level in integrated signal
        NPKI = 0.125*pks[i] + 0.875*NPKI;             
    end

                                            
    #Adjust the threshold with SNR
    if NPKI != 0 || SPKI != 0
        THR_SIG = NPKI + 0.25*(abs(SPKI - NPKI));
        THR_NOISE = 0.5*(THR_SIG);
    end
    
    #Adjust the threshold with SNR for filtered signal
    if NPKF != 0 || SPKF != 0
        THR_SIG1 = NPKF + 0.25*(abs(SPKF - NPKF));
        THR_NOISE1 = 0.5*(THR_SIG1);
    end
    
    
    #Keep track of thresholds of integrated signal
    SIGL_buf = append!(SIGL_buf, SPKI);
    NOISL_buf = append!(NOISL_buf, NPKI);
    THRS_buf = append!(THRS_buf, THR_SIG);

    #Keep track of thresholds of filtered signal
    SIGL_buf1 = append!(SIGL_buf1, SPKF);
    NOISL_buf1 = append!(NOISL_buf1, NPKF);
    THRS_buf1 = append!(THRS_buf1, THR_SIG1);

    Tskip = 0; #reset parameters
    not_nois = 0; #reset parameters
    ser_back = 0;  #reset bandpass param   
end

Plotting the Final Results

In [19]:
plot(sample_motion[:,1],x_MAF,title = "MAF Signal after Algorithm",xlabel="Time (sec)")
plot!(sample_motion[qrs_i,1],qrs_c,line = :scatter)
Out[19]:
In [20]:
plot(sample_motion[:,1],x_HP,title = "Filtered Signal after Algorithm",xlabel="Time (sec)")
plot!(sample_motion[qrs_i_raw,1],qrs_amp_raw,line = :scatter)
Out[20]:

Yay!!!!

So what the heck do all the other parameters that we were keeping track of contain?

Looking at the signal and noise levels and thresholds

In [21]:
plot(sample_motion[:,1],x_MAF,title = "MAF Signal after Algorithm",xlabel="Time (sec)",label = "Integrated Signal")
plot!(sample_motion[qrs_i,1],qrs_c,line = :scatter, label = "Identified Peaks")
plot!(sample_motion[qrs_i,1],SIGL_buf, label = "Signal Level")
plot!(sample_motion[qrs_i,1],NOISL_buf, label = "Noise Level")
plot!(sample_motion[qrs_i,1],THRS_buf, label = "Threshold")
Out[21]:
In [22]:
plot(sample_motion[:,1],x_HP,title = "Filtered Signal after Algorithm",xlabel="Time (sec)",label = "Filtered Signal")
plot!(sample_motion[qrs_i_raw,1],qrs_amp_raw,line = :scatter,label = "Identified Peaks") 
plot!(sample_motion[qrs_i_raw,1],SIGL_buf1,label = "Signal Level")
plot!(sample_motion[qrs_i_raw,1],NOISL_buf1, label = "Noise Level")
plot!(sample_motion[qrs_i_raw,1],THRS_buf1, label = "Threshold")
Out[22]:

So I combined everything so far into a single Pan Tompkins function...

In [23]:
include("pan_tompkins.jl")
include("run_pan_tompkins.jl")
Out[23]:
run_pan_tompkins (generic function with 1 method)

Testing on an Arrhythmia Sample

In [24]:
fs_arrhyth = 360;
sample_arrhyth = CSV.read("207.csv",header = ["Time", "ECG1", "ECG2"], datarow=1, nullable=false)
plot(sample_arrhyth[360*515:360*535,1],sample_arrhyth[360*515:360*535,3],title = "Sample ECG during Arrhythmia",xlabel="Time (sec)")
Out[24]:
In [25]:
#Test on arrhythmia sample
fs = fs_arrhyth;
x = sample_arrhyth[360*1700:360*1730,3];
t = sample_arrhyth[360*1700:360*1730,1];

x_HP, x_MAF, qrs_i, qrs_c, SIGL_buf, NOISL_buf, THRS_buf, qrs_i_raw, qrs_amp_raw, SIGL_buf1, NOISL_buf1, THRS_buf1 = run_pan_tompkins(x,t,fs);
In [26]:
#Plot results
#MAF Signal
plot(t,x_MAF,title = "MAF Signal after Algorithm",xlabel="Time (sec)",label = "Integrated Signal")
plot!(t[qrs_i],qrs_c,line = :scatter, label = "Identified Peaks")
plot!(t[qrs_i],SIGL_buf, label = "Signal Level")
plot!(t[qrs_i],NOISL_buf, label = "Noise Level")
plot!(t[qrs_i],THRS_buf, label = "Threshold")
Out[26]:
In [27]:
#Filtered signal
plot(t,x_HP,title = "Filtered Signal after Algorithm",xlabel="Time (sec)",label = "Filtered Signal")
plot!(t[qrs_i_raw],qrs_amp_raw,line = :scatter,label = "Identified Peaks") 
plot!(t[qrs_i_raw],SIGL_buf1,label = "Signal Level")
plot!(t[qrs_i_raw],NOISL_buf1, label = "Noise Level")
plot!(t[qrs_i_raw],THRS_buf1, label = "Threshold")
Out[27]:

In this case, the signal and noise thresholds aren't particularly useful, but one could imagine a case where knowing the rough thresholds could help with calibration of equipment, or choosing settings etc.

In [ ]: