/* file: sigamp.c	G. Moody	30 November 1991
			Last revised:	27 April 2020

-------------------------------------------------------------------------------
sigamp: Measure signal amplitudes
Copyright (C) 1991-2009 George B. Moody

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, see <http://www.gnu.org/licenses/>.

You may contact the author by e-mail (wfdb@physionet.org) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA).  For updates to this software,
please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________

*/

#include <stdio.h>
#include <math.h>
#include <wfdb/wfdb.h>
#define isqrs
#define map2
#define ammap
#define mamap
#define annpos
#include <wfdb/ecgmap.h>

#define NMAX 300
#define _STRING(X) #X
#define STRING(X) _STRING(X)

/* values for timeunits */
#define SECONDS     1
#define MINUTES     2
#define HOURS       3
#define TIMSTR      4
#define MSTIMSTR    5
#define SHORTTIMSTR 6
#define HHMMSS	    7
#define SAMPLES     8

char *pname;
int namp;
int nsig;
int pflag;		/* if non-zero, print physical units */
int qflag;		/* if non-zero, suppress output of trimmed means */
int timeunits = SECONDS;
int vflag;		/* if non-zero, print individual measurements */
int *v0, *vmax, *vmin, *vv;
double **amp, *vmean, *vsum;
WFDB_Time dt1, dt2, dtw;
WFDB_Anninfo ai;
WFDB_Frequency sfreq;
WFDB_Siginfo *si;
WFDB_Time from = 0L, to = 0L, t;

main(int argc, char **argv)
{
    char *p, *record = NULL, *prog_name(char *s);
    int gvmode = 0, i, j, jlow, jhigh, nmax = NMAX,
	ampcmp(), getptp(WFDB_Time t), getrms(WFDB_Time t);
    void help(void), printamp(WFDB_Time t);

    /* Interpret command-line arguments. */
    pname = prog_name(argv[0]);
    for (i = 1; i < argc; i++) {
	if (*argv[i] == '-') switch(*(argv[i]+1)) {
	  case 'a':
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: annotator name must follow -a\n",
			      pname);
		exit(1);
	    }
	    ai.name = argv[i];
	    break;
	  case 'd':
	    if (++i >= argc-1) {
		(void)fprintf(stderr, "%s: DT1 and DT2 must follow -d\n",
			      pname);
		exit(1);
	    }
	    /* save arg list index, convert arguments to samples later (after
	       record has been opened and sampling frequency is known) */
	    dt1 = i++;
	    break;
	  case 'f':
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: start time must follow -f\n",
			      pname);
		exit(1);
	    }
	    /* save arg list index, convert argument to samples later */
	    from = i;
	    break;
	  case 'h':	/* help requested */
	    help();
	    exit(0);
	    break;
	  case 'H':	/* operate in WFDB_HIGHRES mode */
	    gvmode = WFDB_HIGHRES;
	    break;
	  case 'n':
	    if (++i >= argc) {
		(void)fprintf(stderr,
			      "%s: number of measurements must follow -n\n",
			      pname);
		exit(1);
	    }
	    if ((nmax = atoi(argv[i])) < 1) nmax = NMAX;
	    break;
	  case 'p':
	    pflag = 1;
	    if (*(argv[i]+2) == 'd') timeunits = TIMSTR;
	    else if (*(argv[i]+2) == 'e') timeunits = HHMMSS;
	    else if (*(argv[i]+2) == 'h') timeunits = HOURS;
	    else if (*(argv[i]+2) == 'm') timeunits = MINUTES;
	    else if (*(argv[i]+2) == 'S') timeunits = SAMPLES;
	    else timeunits = SECONDS;
	    break;
          case 'q':
	    qflag = 1;
	    break;
	  case 'r':
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: record name must follow -r\n",
			      pname);
		exit(1);
	    }
	    record = argv[i];
	    break;
	  case 't':
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: stop time must follow -t\n",
			      pname);
		exit(1);
	    }
	    /* save arg list index, convert argument to samples later */
	    to = i;
	    break;
	  case 'v':
	    vflag = 1;
	    break;
	  case 'w':
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: window size must follow -w\n",
			      pname);
		exit(1);
	    }
	    /* save arg list index, convert argument to samples later */
	    dtw = i;
	    break;
	  default:
	    (void)fprintf(stderr, "%s: unrecognized option %s\n",
			  pname, argv[i]);
	    exit(1);
	    break;
	}
	else {
	    (void)fprintf(stderr, "%s: unrecognized argument %s\n",
			  pname, argv[i]);
	    exit(1);
	    break;
	}
    }
    if (record == NULL) {
	help();
	exit(1);
    }

    /* Finish initialization. */
    if (gvmode == 0 && (p = getenv("WFDBGVMODE")))
	gvmode = atoi(p);
    setgvmode(gvmode|WFDB_GVPAD);

    if ((nsig = isigopen(record, NULL, 0)) <= 0) exit(2);
    if ((si = malloc(nsig * sizeof(WFDB_Siginfo))) == NULL ||
	(v0 = malloc(nsig * sizeof(int))) == NULL ||
	(vmax = malloc(nsig * sizeof(int))) == NULL ||
	(vmin = malloc(nsig * sizeof(int))) == NULL ||
	(vv = malloc(nsig * sizeof(int))) == NULL ||
	(amp = malloc(nsig * sizeof(double *))) == NULL ||
	(vmean = malloc(nsig * sizeof(double))) == NULL ||
	(vsum = malloc(nsig * sizeof(double))) == NULL) {
	(void)fprintf(stderr, "%s: insufficient memory\n", pname);
	exit(3);
    }
    if (isigopen(record, si, nsig) != nsig) exit(2);
    for (i = 0; i < nsig; i++) {
	if (si[i].gain == 0.0) si[i].gain = WFDB_DEFGAIN;
	if ((amp[i] = malloc(nmax * sizeof(double))) == NULL) {
	    (void)fprintf(stderr, "%s: insufficient memory\n", pname);
	    exit(3);
	}
    }
    if ((sfreq = sampfreq(NULL)) <= 0.0)
	sfreq = WFDB_DEFFREQ;

    /* Adjust timeunits if starting time is undefined. */
    if (timeunits == TIMSTR || timeunits == HHMMSS) {
	char *p = timstr(0L);
	if (*p != '[') timeunits = HHMMSS;
    }
    if (timeunits == HOURS) sfreq *= 3600.;
    else if (timeunits == MINUTES) sfreq *= 60.;

    if (from > 0L) from = strtim(argv[(int)from]);
    if (to > 0L) to = strtim(argv[(int)to]);
    if (from < 0L) from = -from;
    if (to < 0L) to = -to;
    if (from >= to) to = strtim("e");

    /* Collect amplitudes, print them if qflag or vflag set. */
    if (ai.name) {	/* windows are relative to annotations */
	WFDB_Annotation annot;

	if (dtw > 0L) {
	    (void)fprintf(stderr,
			  "%s: -a and -w options cannot be used together;\n",
			  pname);
	    (void)fprintf(stderr, "  -w option is being ignored\n");
	}
	if (annopen(record, &ai, 1) < 0) exit(2);
	if (dt1 == 0L) dt1 = -(dt2 = strtim("0.05"));
	else {
	    char *p = argv[(int)dt1+1];

	    if (*p == '-') dt2 = -strtim(p+1);
	    else dt2 = strtim(p);
	    p = argv[(int)dt1];
	    if (*p == '-') dt1 = strtim(p+1);
	    else dt1 = -strtim(p);
	    if (dt1 > dt2) {
		WFDB_Time temp = dt1;
		
		dt1 = dt2; dt2 = temp;
	    }
	}
	if (dt1 == dt2) dt2++;
	if (iannsettime(from) < 0) exit(2);
	while (getann(0, &annot) == 0 && namp < nmax &&
	       (to == 0L || annot.time < to)) {
	    if (map1(annot.anntyp) == NORMAL) {
		if (getptp(annot.time) < 0) break;
		if (vflag || qflag) printamp(annot.time);
		namp++;
	    }
	}
    }
    else {		/* windows are consecutive nonoverlapping segments */
	if (dt1 > 0L) {
	    (void)fprintf(stderr,
			  "%s: -d option must be used with -a;\n", pname);
	    (void)fprintf(stderr, "  -d option is being ignored\n");
	}
	if (dtw == 0L || (dtw = strtim(argv[(int)dtw])) <= 0L)
	    dtw = strtim("1");
	if (from > 0L && isigsettime(from) < 0L) exit(2);
	for (t = from; namp < nmax && (to == 0L || t < to); namp++, t += dtw) {
	    if (getrms(t) < 0) break;
	    if (vflag || qflag) printamp(t);
	}	
    }

    /* Calculate trimmed means of collected amplitudes unless -q option set. */
    if (qflag == 0) {
	jlow = namp/20;
	jhigh = namp - jlow;
	if (vflag)
	    (void)printf("Trimmed mean");
	for (i = 0; i < nsig; i++) {
	    double a;

	    qsort((char*)amp[i], namp, sizeof(double), ampcmp);
	    for (a = 0.0, j = jlow; j < jhigh; j++)
		a += amp[i][j];
	    a /= jhigh - jlow;
	    (void)printf("\t%g", pflag ? a/si[i].gain : a);
	}
	(void)printf("\n");
    }
    exit(0);
}

int getrms(WFDB_Time t)
{
    int i, v;
    WFDB_Time tt;

    if (getvec(vv) < nsig) return (-1);

    for (i = 0; i < nsig; i++) {
	vmean[i] = v0[i] = vv[i];
	vsum[i] = 0.0;
    }

    for (tt = 1; tt < dtw; tt++) {
	if (getvec(vv) < nsig) return (-1);
	for (i = 0; i < nsig; i++)
	    vmean[i] += vv[i];
    }

    if (isigsettime(t) < 0 || getvec(vv) < nsig) return (-1);

    /* The quantity vv[i] - v0[i] is normally zero, but may be non-zero if
       the signals are stored in difference format (since isigsettime may
       introduce a DC error in this case).  In the calculation of vmean[i]
       below, this quantity is used to correct any such error. */
    for (i = 0; i < nsig; i++) {
	vmean[i] = vmean[i]/dtw + vv[i] - v0[i];
	v = vv[i] - vmean[i];
	vsum[i] += (double)v*v;
    }

    for (tt = 1; tt < dtw; tt++) {
	if (getvec(vv) < nsig) return (-1);
	for (i = 0; i < nsig; i++) {
	    v = vv[i] - vmean[i];
	    vsum[i] += (double)v*v;
	}
    }

    for (i = 0; i < nsig; i++)
	amp[i][namp] = sqrt(vsum[i]/dtw);

    return (0);
}

int getptp(WFDB_Time t)
{
    int i;
    WFDB_Time tt;

    if (isigsettime(t + dt1) < 0) return (-1);

    if (getvec(vmin) < nsig) return (-1);

    for (i = 0; i < nsig; i++)
	vmax[i] = vmin[i];

    for (tt = dt1; tt < dt2; tt++) {
	if (getvec(vv) < nsig) return (-1);
	for (i = 0; i < nsig; i++) {
	    if (vv[i] < vmin[i]) vmin[i] = vv[i];
	    else if (vv[i] > vmax[i]) vmax[i] = vv[i];
	}
    }

    for (i = 0; i < nsig; i++)
	amp[i][namp] = vmax[i] - vmin[i];

    return (0);
}

void printamp(WFDB_Time t)
{
    int i;

    switch (timeunits) {
      case TIMSTR:  printf("%s",mstimstr(-t)); break;
      case HHMMSS:    printf("%s", mstimstr(t)); break;
      case SAMPLES:   printf("%"WFDB_Pd_TIME, t); break;
      default:	    printf("%lf", t/sfreq); break;
    }
    for (i = 0; i < nsig; i++)
	(void)printf("\t%g", pflag ? amp[i][namp]/si[i].gain :
		     amp[i][namp]);
    (void)printf("\n");
}

int ampcmp(double *p1, double *p2)
{
    if (*p1 > *p2) return (1);
    else if (*p1 == *p2) return (0);
    else return (-1);
}

char *prog_name(char *s)
{
    char *p = s + strlen(s);

#ifdef MSDOS
    while (p >= s && *p != '\\' && *p != ':') {
	if (*p == '.')
	    *p = '\0';		/* strip off extension */
	if ('A' <= *p && *p <= 'Z')
	    *p += 'a' - 'A';	/* convert to lower case */
	p--;
    }
#else
    while (p >= s && *p != '/')
	p--;
#endif
    return (p+1);
}

static char *help_strings[] = {
 "usage: %s -r RECORD [OPTIONS ...]\n",
 "where OPTIONS may include any of:",
 " -a ANN      specify annotator, measure amplitudes near QRS annotations",
 " -d DT1 DT2  set measurement window relative to QRS annotations",
 "              defaults: DT1 = 0.05 (seconds before annotation);",
 "              DT2 = 0.05 (seconds after annotation)",
 " -f TIME     begin at specified time",
 " -h          print this usage summary",
 " -H          read multifrequency signals in high resolution mode",
 " -n NMAX     make up to NMAX measurements per signal (default: NMAX = "
     STRING(NMAX) ")",
 " -p          print results in physical units (default: ADC units)",
 "              -p may be followed by a character to choose a time format:",
 "  -pd         print time of day and date if known",
 "  -pe         print elapsed time as <hours>:<minutes>:<seconds>",
 "  -ph         print elapsed time in hours",
 "  -pm         print elapsed time in minutes",
 "  -ps         print elapsed time in seconds (default)",
 "  -pS         print elapsed time in sample intervals",
 " -q          quick mode: print individual measurements only",
 " -t TIME     stop at specified time",
 " -v          verbose mode: print individual measurements and trimmed means",
 " -w DTW      set RMS amplitude measurement window",
 "              default: DTW = 1 (second)",
NULL
};

void help(void)
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++)
	(void)fprintf(stderr, "%s\n", help_strings[i]);
}