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

main(argc, argv)
int argc;
char *argv[];
{
    int btype, i, j, nbeats = 0, nsig, window;
    long **sum;
    double *dc;
    WFDB_Anninfo a;
    WFDB_Annotation annot;
    WFDB_Sample *v;
    WFDB_Siginfo *s;
    void *calloc();
    WFDB_Time stoptime=0;

    if (argc < 3) {
        fprintf(stderr,
                "usage: %s annotator record beat-type [from to]\n",
                argv[0]);
        exit(1);
    }

    //Get trial length from header file
    char *info;
    char *tlength,*tst0,*tnd;
    char *tag_st="<Trial Length (samples)>:";
    
     if (info = getinfo(argv[2]))
      do {
	tst0 = strstr(info,tag_st);
        tlength= strchr(tst0,':');
        tlength++;
	tnd = strchr(tlength,'<');
	*tnd='\0';
	window=atoi(tlength);
      } while (info = getinfo(NULL));

 
    a.name = argv[1]; a.stat = WFDB_READ;
    if ((nsig = isigopen(argv[2], NULL, 0)) < 1) exit(2);
    s = (WFDB_Siginfo *)malloc(nsig * sizeof(WFDB_Siginfo));
    v = (WFDB_Sample *)malloc(nsig * sizeof(WFDB_Sample));
    sum = (long **)malloc(nsig * sizeof(long *));
    dc = (double *)malloc(nsig * sizeof(double));
    if (s == NULL || v == NULL || dc == NULL || sum == NULL) {
	fprintf(stderr, "%s: insufficient memory\n", argv[0]);
	exit(2);
    }
    if (wfdbinit(argv[2], &a, 1, s, nsig) != nsig) exit(3);
    for (i = 0; i < nsig; i++)
        if ((sum[i]=(long *)calloc(window,sizeof(long))) == NULL) {
            fprintf(stderr, "%s: insufficient memory\n", argv[0]);
            exit(2);
	}

    if(argc > 3) btype = strann(argv[3]);
    if (argc > 4) iannsettime((WFDB_Time) atoi(argv[4]));
    if (argc > 5) stoptime=( (WFDB_Time) atoi(argv[5])); //stop time in samples
    
    while (getann(0, &annot) == 0 && annot.time < window)
        ;

   do {
        if (annot.anntyp != btype) continue;
        isigsettime(annot.time - 1);
        for (j = 0; j < window && getvec(v) > 0; j++)
	   sum[annot.chan][j] += v[annot.chan];
        nbeats++;
    } while (getann(0, &annot) == 0 &&
             (stoptime == 0L || annot.time < stoptime));
  
    if (nbeats < 1) {
        fprintf(stderr, "%s: no `%s' beats found\n",
                argv[0], annstr(btype));
        exit(4);
    }

    
    fprintf(stderr,"Average of %d `%s' beats. Trial length: %d\n", nbeats/nsig, annstr(btype),window);

    //Get DC component of the averages
    for (i=0;i<nsig;i++){
      for (j=0;j<window; j++)
	dc[i]=dc[i]+ sum[i][j];
      dc[i]=dc[i]/window;
    }

    for (j = 0; j < window; j++)
      for (i = 0; i < nsig; i++)
	printf("%g%c", (double)nsig*((double)sum[i][j]-dc[i])/nbeats,
                   (i == nsig-1) ? '\n' : '\t');
      
    exit(0);
}