{
    ### Input file: two columns, time of the R peak in seconds and beat type
    ### Ex.: 0.000 1
    ###      0.808 1
    ###      1.618 1
    ###      2.458 1
 
    CONVFMT = "%.3lf";
    
    NN_min = .300;  # NN intervals <= 0.3 s will be excluded
    NN_max = 1.800; # NN intervals >= 1.8 s will be excluded
    SF = 250;    # ECG sampling frequency
    t = $1;      # Time of the R peak in seconds
    ann = $2     # Annotation: beat type (1: sinus; 2: PVC; 3: PAC)

    #### The time of each R peak is rounded (up or down) to the closest multiple of 1/SF.                                          
    #### Ex.: Let SF = 250 Hz (1/SF = 0.004 s). The time of each R peak should be multiples of 0.004 but sometimes it is not due to conversion/export problems.
    #### If, for example, the time of the R peak was 0.809, this value would be rounded down to 0.808.
    #### Users not interested in this "correction," should comment the next five lines.
    
    aux = int(sprintf(CONVFMT, t*SF));
    if (t - aux/SF > 0.5)
	t = (aux+1)/SF;
    if (t - aux/SF  <= 0.5 && t - aux/SF > 0)
	t = aux/SF;
   
    #### Calculation of NN and delta NN intervals after reading 4 data points
    NN = t - t_1;     # Current NN interval
    NN_1 = t_1 - t_2; # Previous NN interval
    NN_2 = t_2 - t_3; # NN interval preceeding the "previous NN interval"
    dNN_formatted = sprintf(CONVFMT, t - 2*t_1 + t_2);      # Current delta NN interval
    dNN = dNN_formatted + 0;
    dNN_1_formatted = sprintf(CONVFMT, t_1 - 2*t_2 + t_3);  # Previous delta NN interval
    dNN_1 = dNN_1_formatted + 0;
    
    #### Verification of whether the current beat and the 3 preceeding ones are "normal sinus" beats (ann == 1)                                        
    if (ann == 1  &&  ann_1 == 1  &&  ann_2 == 1  &&  ann_3 == 1)
	flag_1 = 0;
    else{
	flag_1 = 1;
	non_sinus = 1;
    }
    if (NR <= 4 && (ann < 1 || ann_1 < 1 || ann_2 < 1 || ann_3 < 1))
	non_sinus = 0;

    #### Verification of whether current NN interval and the 3 preceeding one are between NN_min (0.3 sec) and NN_max (1.8 sec)
    if (NN < NN_max  &&  NN > NN_min  &&  NN_1 < NN_max  &&  NN_1 > NN_min  &&  NN_2 < NN_max  &&  NN_2 > NN_min)
	flag_2 = 0;
    else{
	flag_2 = 1;
	outlier = 1; 
    }
    if (NR <= 4 && (NN == 0 || NN_1 == 0 || NN_2 == 0 || NN_3 == 0))
	outlier = 0;

    # Counter for the number of inflection points: IP 
    # Counter for the number of delta NN intervals: n_dNN
    # Counter for the number of delta NN intervals in short (long) accelerative and decelerative segments: short (long) 
    # Counter for the number of delta NN intervals between consecutive inflection points: cont

    #### Hard inflection point criterion: dNN_current * dNN_previous < 0. 
    if (NR > 3 && flag_1 == 0 && flag_2 == 0 && dNN * dNN_1 < 0){
	IP++;
	if (non_sinus != 1 && outlier != 1){
	    if (IP > 1)
		n_dNN += cont;   
	    if (cont <= 2 && IP > 1)
		short += cont;
	    if (cont > 2 && IP > 1)
		long += cont;
#	    print "Hard IP: time =", t_1, " \#  Beat =", NR-2, "  NN =", NN_1;
	    if (IP == 1)
		NN_total = 1;
	    else
		NN_total += cont;
#	    print ("NN_total =", NN_total); 
	}
	if(non_sinus == 1 || outlier == 1){
	    non_sinus = 0; 
	    outlier = 0;
	    NN_total += 1
#	    print "Hard IP: time =", t_1, "  \# Beat =", NR-2, "  NN =", NN_1;
#	    print "NN_total =", NN_total;
	}
	cont = 0;
    }

    #### Soft inflection point criterion: dNN * dNN_prev == 0 and dNN != dNN_prev
    if (NR > 3 && flag_1 == 0 && flag_2 == 0 && dNN * dNN_1 == 0  &&  dNN != dNN_1){
	IP++;
	if (non_sinus != 1 && outlier != 1){
	    if (IP > 1)
		n_dNN += cont;
	    if (dNN_1 != 0) 
		if (cont <= 2 && IP > 1)
		    short += cont;
	    if (cont > 2 && IP > 1)
		    long += cont;
#	    print "Soft IP: time =", t_1, "  \# Beat =", NR-2, "  NN =", NN_1;
	    if (IP == 1)
		NN_total = 1;
	    else
		NN_total += cont;
#	    print ("NN_total =", NN_total); 
	}
	if((non_sinus == 1 || outlier == 1)){
	    non_sinus = 0;
	    outlier = 0;
	    NN_total += 1;
#	    print "Soft IP: time =", t_1, "  \# Beat =", NR-2, "  NN =", NN_1;
#	    print "NN_total =", NN_total;
	}
	cont = 0;
    }
    
    if (NR > 1)
	cont++;
    
    ann_3 = ann_2;
    ann_2 = ann_1;
    ann_1 = ann;
    t_3 = t_2;
    t_2 = t_1;
    t_1 = t;
    N_acc_dec = short + long;

    if (n_dNN > 0)
	PIP = (sprintf(CONVFMT,IP/NN_total*100));
    if (N_acc_dec > 0)
	PNNSS = (sprintf(CONVFMT,short/N_acc_dec*100));
    if (n_dNN > 0)
	PNNLS = (sprintf(CONVFMT,long/n_dNN*100));
}
END{print "PIP =", IP"/"NN_total, "=", PIP";", "  PNNSS =", short"/"N_acc_dec"*100 =", PNNSS";", "  PNNLS =", long"/"n_dNN"*100 =", PNNLS}