/* file: edr.c		G. Moody	13 February 1986
			Last revised:	21 January 2004

-------------------------------------------------------------------------------
edr: Derive a respiration signal from an ECG
Copyright (C) 1986-2004 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, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.

You may contact the author by e-mail (george@mit.edu) 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 derives a sample of a respiratory signal for each QRS complex, by
measuring the mean electrical axis (in two-channel mode) or the projection of
that axis onto the lead axis (in single-channel mode).  See "Derivation of
respiratory signals from multi-lead ECGs", pp. 113-116, Computers in Cardiology
1985.  Input to this program includes an annotation file in which each beat to
be used for an EDR sample has been labelled, in addition to the signal and
header files.  The output is another annotation file, which is a copy of the
input annotation file except that the 'num' field of each beat annotation is
replaced by an EDR sample.

If the beat annotations are not located at the QRS peaks, it will be necessary
to set the window limits (the offsets relative to the annotations between
which the raw measurements for the EDR are taken), using the -d option.
By default, this program behaves as if the option "-d -0.04 0.04" is given
(in other words, measurements are taken over an 80 msec window beginning
40 msec -- .04 seconds -- before the annotation, and ending 40 msec after the
annotation).  If edr is supplied with annotations generated by sqrs, the
option "-d 0 0.08" is recommended.

Compile this program with the WFDB library (from http://www.physionet.org/)
and with the math library (-lm), e.g.:
	cc -o edr -O edr.c -lm -lwfdb
*/

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

#define VBL	2048	/* VBL must be a power of 2 */

#ifdef M_PI
#define PI M_PI	/* pi to machine precision */
#else
#define PI 3.141592653589793238462643383279502884197169399375105820974944
			/* quite a few more digits than we need :-) */
#endif

char *pname;
double x, y;
int blen = 0, nsig = 2, *sig;
WFDB_Sample dt1 = 0L, dt2;

int main(int argc, char **argv)
{
    char *record = NULL, *prog_name();
    int i, isiglist = 0, nasig = 0, s, vflag = 0, edr();
    WFDB_Sample dt1 = 0L, dt2, from = -1, to = -1;
    static WFDB_Annotation annot;
    static WFDB_Anninfo ai[2];
    static WFDB_Siginfo si[WFDB_MAXSIG];
    void getxy(), help();

    pname = prog_name(argv[0]);
    ai[0].stat = WFDB_READ;
    ai[1].name = "edr";
    ai[1].stat = WFDB_WRITE;
    for (i = 1; i < argc; i++) {
	if (*argv[i] == '-') switch (*(argv[i]+1)) {
	  case 'd':	/* window offsets from fiducial follow */
	    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':	/* start time follows */
	    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 'i':	/* input annotator follows */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: input annotator must follow -i\n",
			      pname);
		exit(1);
	    }
	    ai[0].name = argv[i];
	    break;
	  case 'o':	/* output annotator follows */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: output annotator must follow -o\n",
			      pname);
		exit(1);
	    }
	    ai[1].name = argv[i];
	    break;
	  case 'r':	/* record name follows */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: record name must follow -r\n",
			      pname);
		exit(1);
	    }
	    record = argv[i];
	    break;
	  case 's':	/* signal list follows */
	    isiglist = i+1; /* index of first argument containing a signal # */
	    while (i+1 < argc && *argv[i+1] != '-') {
		i++;
		nasig++;	/* number of elements in signal list */
	    }
	    if (nasig == 0) {
		(void)fprintf(stderr, "%s: signal list must follow -s\n",
			pname);
		exit(1);
	    }
	    break;
	  case 't':	/* end time follows */
	    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':	/* verbose mode */
	    vflag = 1;
	    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 || ai[0].name == NULL) {
	help();
	exit(1);
    }
    if ((nsig = isigopen(record, si, WFDB_MAXSIG)) <= 0) exit(2);
    if (nasig) {	/* analyze specified signals only */
	if ((sig = (int *)malloc((unsigned)nasig*sizeof(int))) == NULL) {
	    (void)fprintf(stderr, "%s: insufficient memory\n", pname);
	    exit(2);
	}
	for (i = 0; i < nasig; i++) {
	    if ((s = atoi(argv[isiglist+i])) < 0 || s >= nsig) {
		(void)fprintf(stderr, "%s: can't read signal %d\n", pname, s);
		exit(2);
	    }
	    sig[i] = s;
	}
	nsig = nasig;
    }
    else {	/* analyze all signals */
	if ((sig = (int *)malloc((unsigned)nsig*sizeof(int))) == NULL) {
	    (void)fprintf(stderr, "%s: insufficient memory\n", pname);
	    exit(2);
	}
	for (i = 0; i < nsig; i++)
	    sig[i] = i;
    }
    if (nsig > 2) {
	(void)fprintf(stderr,
		      "%s: warning: only signals %d and %d will be analyzed\n",
		      pname, sig[0], sig[1]);
	nsig = 2;
    }
    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");
    blen = strtim("1");	/* set half-length of baseline filter to 1 second */
    if (VBL < 2*blen+1) {
	fprintf(stderr, "%s: buffer length, %d, is too short\n", pname, VBL);
	exit(1);
    }
    if (annopen(record, ai, 2) < 0) exit(2);
    if (dt1 == 0L) dt1 = -(dt2 = strtim("0.04"));
    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) {
	    long temp = dt1;
	    
	    dt1 = dt2; dt2 = temp;
	}
    }
    if (dt1 == dt2) dt2++;
    if (iannsettime(from) < 0) exit(2);
    while (getann(0, &annot) == 0 && (to == 0L || annot.time < to)) {
	if (isqrs(annot.anntyp)) {
	    getxy(annot.time+dt1, annot.time+dt2);
	    annot.num = edr();
	    if (vflag)
		printf("%8ld %4d %g %g %s\n",
		       annot.time, annot.num, x, y, annstr(annot.anntyp));
	}
	if (putann(0, &annot) != 0) {
	    wfdbquit();
	    exit(3);
	}
    }
    wfdbquit();
    exit(0);
}

static long bbuf[2][VBL];

int baseline(int s, WFDB_Time t)
{
    int c, i;
    static WFDB_Time t0;

    if (s < 0 || s >= nsig) return (0);
    if (t0 == 0L)
	for (c = 0, bbuf[c][0] = sample(c, 0L) * (2*blen+1); c < nsig; c++)
	    for (i = 1; i < VBL; i++)
		bbuf[c][i] = bbuf[c][0];
    while (t0 < t)
	for (c = 0; c < nsig; c++)
	    bbuf[c][++t0 & (VBL-1)] += sample(c, t+blen) - sample(c, t-blen);
    return ((int)(bbuf[s][t & (VBL-1)]/blen));
}

void getxy(WFDB_Time t0, WFDB_Time t1)
{
    for (x = y = 0.0; t0 <= t1; t0++) {
	x += sample(0, t0) - baseline(0, t0);
	if (nsig > 1) y += sample(1, t0) - baseline(1, t0); 
    }
}

int edr()			/* calculate instantaneous EDR */
{
    double d, dn, r, theta;
    static int xc, yc, thc;
    static double xd, yd, td, xdmax, ydmax, tdmax, xm, ym, tm;

    if (x == 0 && y == 0) return (0);
    switch (nsig) {
      case 1:
	d = x - xm;
	if (xc < 500) dn = d/++xc;
	else if ((dn = d/xc) > xdmax) dn = xdmax;
	else if (dn < -xdmax) dn = -xdmax;
	xm += dn;
	xd += fabs(dn) - xd/xc;
	if (xd < 1.) xd = 1.;
	xdmax = 3.*xd/xc;
	r = d/xd;
	break;
      case 2:
	d = x - xm;
	if (xc < 500) dn = d/++xc;
	else if ((dn = d/xc) > xdmax) dn = xdmax;
	else if (dn < -xdmax) dn = -xdmax;
	xm += dn;
	xd += fabs(dn) - xd/xc;
	if (xd < 1.) xd = 1.;
	xdmax = 3.*xd/xc;
	d = y - ym;
	if (yc < 500) dn = d/++yc;
	else if ((dn = d/yc) > ydmax) dn = ydmax;
	else if (dn < -ydmax) dn = -ydmax;
	ym += dn;
	yd += fabs(dn) - yd/yc;
	if (yd < 1.) yd = 1.;
	ydmax = 3.*yd/yc;
	theta = atan2(x, y);
	d = theta - tm;
	if (d > PI) d -= 2.*PI;
	else if (d < -PI) d += 2.*PI;
	if (thc < 500) dn = d/++thc;
	else if ((dn = d/thc) > tdmax) dn = tdmax;
	else if (dn < -tdmax) dn = -tdmax;
	if ((tm += dn) > PI) tm -= 2.*PI;
	else if (tm < -PI) tm += 2.*PI;
	td += fabs(dn) - td/thc;
	if (td < 0.001) td = 0.001;
	tdmax = 3.*td/thc;
	r = d/td;
	break;
    }
    return ((int)(r*50.));
}

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 -i ANN [OPTIONS ...]\n",
 " where RECORD and ANN specify the inputs, and OPTIONS may include any of:",
 " -d DT1 DT2  set measurement window relative to QRS annotations",
 "              defaults: DT1 = 0.04 (seconds before annotation);",
 "              DT2 = 0.04 (seconds after annotation)",
 " -f TIME     begin at specified time",
 " -h          print this usage summary",
 " -o ANN      specify output annotator (default: 'edr')",
 " -s SIGNAL [SIGNAL ...]  analyze only the specified signal(s)",
 " -t TIME     stop at specified time",
 " -v          verbose mode: print individual measurements",
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]);
}