--- a
+++ b/HRV.src/statnn.c
@@ -0,0 +1,262 @@
+/* statnn.c		Joe Mietus		Oct 7 2008 */
+
+/* 
+statnn :
+Usage: statnn [options]
+Reads stdin: time(sec), interval(sec), annotation
+and calculates nn interval stats
+Outputs nn/rr, avnn, sdnn, sdann,
+sdnnindx, rmssd and pnn on stderr
+  options :
+  [-l len] : window length (default 5:00)
+  [-m] : RR intervals in msec
+  [-p nndif ...] : nn diffenence for pnn
+  [-s] : short term stats
+         nn/rr, avnn, sdnn, rmssd and pnn on stderr
+  [-L] : print ratio avnn sdnn sdann sdnnindx rmssd pnns on one line
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define MAXBUF 2400
+#define MAXP 10
+#define NNDIF 0.05
+#define ABS(A) ((A) < 0 ? -(A) : (A))
+#define ROFFERR 1e-10
+
+main(argc, argv)
+int argc;
+char *argv[];
+{
+
+    char ann, lastann[2];
+    int i, j, len, lflag, mflag, sflag;
+    int n, nb, nrr, totrr, nnn, totnn, totnnn, np;
+    static int nnx[MAXP];
+    double t, lastt, rr, lastrr, start, end, sum, sum2;
+    double ratio, avnn, sdnn, sdann, sdnnindx, rmssd, pnnx;
+    static double ratbuf[MAXBUF], avbuf[MAXBUF], sdbuf[MAXBUF];
+    static double nndif[MAXP];
+
+    /* default 5:00 non-overlapping segments */
+    len = 300;
+
+    lflag = mflag = sflag = 0;
+    np = 1;
+    nndif[0] = NNDIF;
+
+    for (i=0; ++i < argc && *argv[i] == '-'; ) {
+        switch(argv[i][1]) {
+	    case 'l':  len = strtim(argv[++i]);
+	               break;
+	    case 'm':  mflag = 1;
+	               break;
+	    case 'p':  for (np=0; np<MAXP && i+1<argc && *argv[i+1]!='-'; np++)
+	                   nndif[np] = atof(argv[++i]);
+	               if (np==MAXP)
+			   while (i+1<argc && *argv[i+1] != '-')
+			       i++;
+	               break;
+            case 's':  sflag = 1;
+	               break;
+	    case 'L':  lflag = 1;
+	               break;
+            default :  usage(argv[0]);
+                       exit(1);
+        }
+    }
+
+    if (len < 60) {
+        fprintf(stderr, "%s : len must be 1:00 or greater\n", argv[0]);
+        exit(1);
+    }
+
+    for (j=0; j<np; j++) {
+        if (nndif[j] > 1)
+	    nndif[j] /= 1000;
+	if (mflag)
+	    nndif[j] *= 1000;
+    }
+
+    if (scanf("%lf %lf %c", &t, &rr, &ann) != 3) {
+        fprintf(stderr, "%s : improperly formatted data\n", argv[0]);
+        exit(2);
+    }
+
+    i = 0;
+    totrr = nrr = 1;
+    totnnn = totnn = nnn = 0;
+    rmssd = sum = sum2 = 0.0;
+
+    lastt = t;
+    lastrr = rr;
+    lastann[0] = 'X';
+    lastann[1] = ann;
+    end = t + len;
+
+    while (scanf("%lf %lf %c", &t, &rr, &ann) == 3) {
+
+	while (t > end+len) {
+	    i++;
+	    end += len;
+	}
+
+	if (t >= end) {
+	    if (nnn > 1) {
+		ratbuf[i] = (double)nnn/nrr;
+	        sdbuf[i] = sqrt((sdbuf[i] - avbuf[i]*avbuf[i]/nnn) / (nnn-1));
+	        avbuf[i] /= nnn;
+	    }
+	    i++;
+	    nnn = nrr = 0;
+	    end += len;
+        }
+
+        nrr++;
+	totrr++;
+
+        if (ann == 'N' && lastann[1] == 'N') {
+            nnn++;
+	    totnn++;
+	    avbuf[i] += rr;
+	    sum += rr;
+	    sdbuf[i] += rr*rr;
+	    sum2 += rr*rr;
+	    if (lastann[0] == 'N') {
+	      totnnn++;
+	        rmssd += (rr-lastrr)*(rr-lastrr);
+
+		/* modified 6-15-01 */
+		/* equal may not be equal do to round of error !! */
+	        /* if (ABS(rr-lastrr) > nndif) { */
+
+		for (j=0; j<np; j++) {
+		    if (ABS(rr-lastrr) - nndif[j] > ROFFERR) { 
+		        nnx[j]++;
+	            }
+		}
+	    }
+	}
+
+	lastrr = rr;
+	lastann[0] = lastann[1];
+	lastann[1] = ann;
+
+    }
+
+    if (nnn > 1) {
+	ratbuf[i] = (double)nnn/nrr;
+        sdbuf[i] = sqrt((sdbuf[i] - avbuf[i]*avbuf[i]/nnn) / (nnn-1));
+        avbuf[i] /= nnn;
+    }
+
+    nb = ++i;
+
+    if (lflag) {
+        /* print NN/RR AVNN SDNN SDANN SDNNIDX rMSSD pNNs on one line */
+        printf("%g", (double)totnn/totrr);
+        printf(" %g", sum/totnn);
+        printf(" %g", sqrt((sum2 - sum*sum/totnn) / (totnn-1)));
+    }
+    else {
+        printf("NN/RR = %g\n", (double)totnn/totrr);
+        printf("AVNN = %g\n", sum/totnn);
+        printf("SDNN = %g\n", sqrt((sum2 - sum*sum/totnn) / (totnn-1)));
+    }
+
+    if (!sflag) {
+        sum = sum2 = 0;
+        for (i=0, n=0; i<nb; i++) {
+            if (ratbuf[i] != 0) {
+                 n++;
+		 sum += avbuf[i];
+	         sum2 += avbuf[i]*avbuf[i];
+	    }
+        }
+        if (lflag) {
+            if (n > 1)
+                printf(" %g", sqrt((sum2 - sum*sum/n) / (n-1)));
+            else
+                printf(" -");
+        }
+        else {
+            if (n > 1)
+                printf("SDANN = %g\n", sqrt((sum2 - sum*sum/n) / (n-1)));
+            else
+                printf("SDANN = -\n");
+        }
+
+        sum = 0;
+        for (i=0, n=0; i<nb; i++) {
+            if (ratbuf[i] != 0) {
+	        n++;
+	        sum += sdbuf[i];
+	    }
+        }
+        if (lflag) {
+            if (n > 0)
+                printf(" %g", sum/n);
+            else
+                printf(" -");
+        }
+        else {
+            if (n > 0)
+                 printf("SDNNIDX = %g\n", sum/n);
+            else
+                printf("SDNNIDX = -\n");
+        }
+    }
+
+    if (lflag) {
+        printf(" %g", sqrt(rmssd/totnnn));
+        for (j=0; j<np; j++)
+            printf(" %g", mflag ? 100*(double)nnx[j]/totnnn
+                                : (double)nnx[j]/totnnn);
+        printf("\n");
+    }
+    else {
+        printf("rMSSD = %g\n", sqrt(rmssd/totnnn));
+        for (j=0; j<np; j++) {
+            printf("pNN%g = ", mflag ? nndif[j] : nndif[j]*1000);
+            printf("%g\n", mflag ? 100*(double)nnx[j]/totnnn
+		                 : (double)nnx[j]/totnnn);
+        }
+    }
+}
+
+
+/* convert string in [[HH:]MM:]SS format to seconds */
+
+strtim(buf)
+char *buf;
+{
+        int x, y, z;
+
+        switch (sscanf(buf, "%d:%d:%d", &x, &y, &z)) {
+                case 1: return (x);
+                case 2: return (60*x + y);
+                case 3: return (3600*x + 60*y + z);
+                default: return (-1);
+        }
+}
+
+
+usage(prog)
+char *prog;
+{
+        fprintf(stderr, "Usage: %s [options]\n", prog);
+        fprintf(stderr, "Reads stdin: time(sec), interval(sec), annotation\n");
+        fprintf(stderr, "and calculates nn interval stats\n");
+        fprintf(stderr, "Outputs nn/rr, avnn, sdnn, sdann,\n");
+        fprintf(stderr, "sdnnindx, rmssd and pnn on stderr\n");
+	fprintf(stderr, "  options :\n");
+        fprintf(stderr, "  [-l len] : window length (default 5:00)\n");
+        fprintf(stderr, "  [-m] : RR intervals in msec\n");
+        fprintf(stderr, "  [-p nndif ...] : nn diffenence for pnn\n");
+        fprintf(stderr, "  [-s] : short term stats\n");
+        fprintf(stderr, "         nn/rr, avnn, sdnn, rmssd and pnn on stderr\n");
+	fprintf(stderr, "  [-L] : print ratio avnn sdnn sdann sdnnindx");
+	fprintf(stderr, " rmssd pnns on one line\n");
+}