/* 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");

}