/* file: sigavg.c		G. Moody		25 November 2002

-------------------------------------------------------------------------------
sigavg: Calculate averages of annotated waveforms
Copyright (C) 2002 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/).
_______________________________________________________________________________

This program is based on Example 9 in the WFDB Programmer's Guide.
*/

#include <stdio.h>
#include <wfdb/wfdb.h>
#include <wfdb/ecgmap.h>

char *irec;	/* input record name */
char *pname;	/* name by which this program was invoked */
double sps;	/* sampling frequency (Hz) */
double **sum;
int dt1 = 0, dt2 = 0, flag[ACMAX+1], Hflag, nsig, vflag = 0, zflag = 0;
WFDB_Anninfo ai;
WFDB_Sample *v, *vb;
WFDB_Siginfo *s;
WFDB_Time from = 0L, to = 0L;

char *prog_name();
void help(), init(), memerr();

main(argc, argv)
int argc;
char *argv[];
{
    int i, j, nbeats = 0;
    WFDB_Annotation annot;

    init(argc, argv);	/* read and interpret command line */

    while (getann(0, &annot) == 0 && annot.time < from)
        ;
    do {
        if (flag[annot.anntyp] == 0) continue;
        isigsettime(annot.time + dt1 - 1);
        getvec(vb);
        for (j = dt1; j <= dt2  && getvec(v) > 0; j++) {
	    if (zflag)
		for (i = 0; i < nsig; i++)
		    sum[i][j-dt1] += v[i] - vb[i];
	    else
		for (i = 0; i < nsig; i++)
		    sum[i][j-dt1] += v[i];
	}
        nbeats++;
    } while (getann(0, &annot) == 0 &&
             (to == 0L || annot.time < to));
    if (nbeats < 1) {
        fprintf(stderr, "%s: no beats found\n", pname);
        exit(4);
    }
    if (vflag) {
	printf("# Average of %d beats:\n", nbeats);
	printf("#     Time\t");
	for (i = 0; i < nsig; i++)
	    printf("%10s%c", s[i].desc, (i == nsig-1) ? '\n' : '\t');
	printf("#      sec\t");
	for (i = 0; i < nsig; i++) {
	    if (s[i].units == NULL) s[i].units = "mV";
	    printf("%10s%c", s[i].units, (i == nsig-1) ? '\n' : '\t');
	}
    }
    for (j = dt1; j < dt2; j++) {
	printf("%10.5lf\t", (double)j/sps);
        for (i = 0; i < nsig; i++) {
	    double m;

	    m = sum[i][j-dt1]/nbeats;
	    if (!zflag) m -= s[i].baseline;
	    m /= s[i].gain;
            printf("%10.5lf%c", m, (i == nsig-1) ? '\n' : '\t');
	}
    }
    exit(0);
}

void init(int argc, char **argv)
{
    int i, j;

    pname = prog_name(argv[0]);
    for (i = 1; i < ACMAX; i++)
	flag[i] = isqrs(i);
    for (i = 1; i < argc; i++) {
	if (*argv[i] == '-') switch (*(argv[i]+1)) {
	  case 'a':	/* annotator follows */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: annotator must follow -a\n",
			      pname);
		exit(1);
	    }
	    ai.name = argv[i];
	    ai.stat = WFDB_READ;
	    break;
	  case 'd':
	    if (++i >= argc-1) {
		(void)fprintf(stderr, "%s: window offsets must follow -w\n",
			      pname);
		exit(1);
	    }
	    /* save arg list indices, convert argument to samples later */
	    dt1 = i;
	    dt2 = ++i;
	    break;
	  case 'f':	/* starting time */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: time must follow -f\n", pname);
		exit(1);
	    }
	    from = i;
	    break;
	  case 'h':	/* help requested */
	    help();
	    exit(0);
	    break;
	  case 'H':   /* read multifrequency records in high resolution mode */
	    Hflag = 1;
	    break;
	  case 'p':	/* annotation mnemonic(s) follow */
	    for (j = 1; j < ACMAX; j++)
		flag[j] = 0;
	    if (++i >= argc || !isann(j = strann(argv[i]))) {
		(void)fprintf(stderr,
			      "%s: annotation mnemonic(s) must follow -p\n",
			      pname);
		exit(1);
	    }
	    flag[j] = 1;
	    /* The code above not only checks that there is a mnemonic where
	       there should be one, but also allows for the possibility that
	       there might be a (user-defined) mnemonic beginning with `-'.
	       The following lines pick up any other mnemonics, but assume
	       that arguments beginning with `-' are options, not mnemonics. */
	    while (++i < argc && argv[i][0] != '-')
		if (isann(j = strann(argv[i]))) flag[j] = 1;
	    if (i == argc || argv[i][0] == '-') i--;
	    flag[0] = 0;
	    break;
	  case 'r':	/* record name */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: record name must follow -r\n",
			      pname);
		exit(1);
	    }
	    irec = argv[i];
	    break;
	  case 't':	/* end time */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: time must follow -t\n",pname);
		exit(1);
	    }
	    to = i;
	    break;
	  case 'v':	/* verbose mode */
	    vflag = 1;
	    break;
	  case 'z':	/* set the baseline to zero */
	    zflag = 1;
	    break;
	  default:
	    (void)fprintf(stderr, "%s: unrecognized option %s\n", pname,
			  argv[i]);
	    exit(1);
	}
	else {
	    (void)fprintf(stderr, "%s: unrecognized argument %s\n", pname,
			  argv[i]);
	    exit(1);
	}
    }
    if (irec == NULL || ai.name == NULL) {
	help();
	exit(1);
    }
    if (Hflag)
	setgvmode(WFDB_HIGHRES);
    if ((nsig = isigopen(irec, NULL, 0)) <= 0)
	exit(2);
    s = malloc((size_t)nsig * sizeof(WFDB_Siginfo));
    v = (WFDB_Sample *)malloc((size_t)nsig * sizeof(WFDB_Sample));
    vb = (WFDB_Sample *)malloc((size_t)nsig * sizeof(WFDB_Sample));
    sum = (double **)malloc(nsig * sizeof(double *));
    if (s == NULL || v == NULL || vb == NULL || sum == NULL)
	memerr();
    if (wfdbinit(irec, &ai, 1, s, nsig) != nsig)
	exit(2);
    sps = sampfreq(NULL);
    dt1 = strtim(dt1 ? argv[dt1] : "-.05");
    dt2 = strtim(dt2 ? argv[dt2] : ".05");
    if (dt1 > dt2) {
	int tmp;
	tmp = dt1; dt1 = dt2; dt2 = tmp;
    }
    else if (dt1 == dt2)
	dt2 += strtim(".1");
    for (i = 0; i < nsig; i++)
        if ((sum[i]=(double *)calloc(dt2-dt1+1,sizeof(double))) == NULL) {
	    wfdbquit();
	    memerr();
        }
    if (from > 0L) {
	from = strtim(argv[from]);
	if (from < (WFDB_Time)0) from = -from;
	(void)isigsettime(from);
    }
    if (to > 0L) {
	to = strtim(argv[to]);
	if (to < (WFDB_Time)0) to = -to;
    }
}

void memerr()
{
    (void)fprintf(stderr, "%s: insufficient memory\n", pname);
    exit(2);
}

char *prog_name(s)
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 -a ANNOTATOR [OPTIONS ...]\n",
 "where RECORD and ANNOTATOR specify the input, and OPTIONS may include:",
 " -d DT1 DT2          calculate the average over a window defined by offsets",
 "                      DT1 and DT2 from the input annotations (defaults:",
 "                      DT1 = -0.05, DT2 = 0.05 (seconds from annotations)",
 " -f TIME             begin at specified time (default: beginning of RECORD)",
 " -h                  print this usage summary",
 " -H                  read multifrequency records in high resolution mode",
 " -p TYPE [TYPE ...]  include annotations of specified TYPEs only (default:",
 "                      include all QRS annotations",
 " -t TIME             stop at specified time (default: end of RECORD)",
 " -v                  verbose mode",
 " -z                  set the baseline to zero before averaging",
NULL
};

void help()
{
    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]);
}