/*addavg.c - add GnuroScan averaged EEG data files.
  Copyright (c) 1996 Matthew Belmonte

  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.

  If you find this program useful, please send mail to Matthew Belmonte.
  <mkb4@Cornell.edu>.  If you base a publication on data processed by this
  program, please notify Matthew Belmonte and include the following citation
  in your publication:

	Matthew Belmonte, `A Software System for Analysis of
	Steady-State Evoked Potentials', Association for Computing
	Machinery SIGBIO Newsletter 17:1:9-14 (April 1997).
*/

#define MAX_NUM_INPUTS 22
#define MAX_NUM_CHANNELS 16
/*Might want to add support for truncation to a common epoch length and
  resampling to a common sampling rate.*/

#include <stdio.h>
#include <unistd.h>
#include <errno.h>
#include "intel.h"
#include "gnuro.h"

/*Open the GnuroScan .avg file 'filename'.  Store the resulting file pointer in
  the location pointed to by 'fpp'.  Extract the number of channels, the number
  of samples, and the variance-data-included flag from the file header and store
  them in the locations pointed to by 'num_channels', 'num_samples', and
  'variance_flag', respectively.  For each channel, extract the baseline and add
  it into the vector pointed to by 'baseline_sum' and extract the calibration
  coefficient and store it in the vector pointed to by 'calibration'.  Return -1
  for failure to open the file, 0 for success.*/
int open_input(filename, fpp, num_channels, num_samples, baseline_sum, calibration, variance_flag)
char *filename;
FILE **fpp;
int *num_channels, *num_samples;
short *baseline_sum;
float *calibration;
int *variance_flag;
  {
  register int channel;
  if((*fpp = fopen(filename,
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    fprintf(stderr, "couldn't open %s\n", filename);
    return(-1);
    }
  fseek(*fpp, (long)S_pnts_offset, SEEK_SET);
  *num_samples = read_Intel_short(*fpp);
  fseek(*fpp, (long)S_nchans_offset, SEEK_SET);
  *num_channels = read_Intel_short(*fpp);
  fseek(*fpp, (long)S_variance_offset, SEEK_SET);
  *variance_flag = getc(*fpp);
  fseek(*fpp, (long)packed_sizeof_SETUP+EL_baseline_offset, SEEK_SET);
  if(*num_channels > MAX_NUM_CHANNELS)
    /*This error condition is detected and reported in the calling routine.*/
    *num_channels = MAX_NUM_CHANNELS;
  for(channel = 0; channel != *num_channels; channel++)
    {
    baseline_sum[channel] += read_Intel_short(*fpp);
    fseek(*fpp, EL_calib_offset-EL_baseline_offset-2L, SEEK_CUR);
    calibration[channel] = read_Intel_float(*fpp);
    fseek(*fpp, packed_sizeof_ELECTLOC+EL_baseline_offset-EL_calib_offset-sizeof(float), SEEK_CUR);
    }
  return(0);
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int sample, channel, file;
  int num_samples, num_channels, variance_flag;
  short baseline[MAX_NUM_CHANNELS];
  float calibration_ratio[MAX_NUM_INPUTS][MAX_NUM_CHANNELS];
  FILE *out, *(in[MAX_NUM_INPUTS]);
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc < 4)
    {
    fprintf(stderr, "usage: %s <infile.avg> ... <infile.avg> <outfile.avg>\n", *argv);
    exit(1);
    }
  if(argc > 2+MAX_NUM_INPUTS)
    {
    fprintf(stderr, "Recompile %s with MAX_NUM_INPUTS=%d\n", *argv, argc-2);
    exit(1);
    }
  for(channel = 0; channel != MAX_NUM_CHANNELS; channel++)
    baseline[channel] = 0;
  if(open_input(argv[1], in, &num_channels, &num_samples, baseline, calibration_ratio[0], &variance_flag))
    exit(errno);
  if(num_channels > MAX_NUM_CHANNELS)
    {
    fclose(*in);
    fprintf(stderr, "Recompile %s with MAX_NUM_CHANNELS=%d\n", *argv, num_channels);
    exit(1);
    }
/*open input files and check consistency of #channels and #samples*/
  for(file = 1; file != argc-2; file++)
    {
    int chns, pnts, vflag;
    if(open_input(argv[1+file], in+file, &chns, &pnts, baseline, calibration_ratio[file], &vflag))
      {
      while(file != -1)
	fclose(in[file--]);
      exit(errno);
      }
    if(chns != num_channels)
      {
      while(file != -1)
	fclose(in[file--]);
      fprintf(stderr, "incompatible number of channels %d in average %s\n", chns, argv[1+file]);
      exit(1);
      }
    if(pnts != num_samples)
      {
      while(file != -1)
	fclose(in[file--]);
      fprintf(stderr, "incompatible number of points %d in average %s\n", pnts, argv[1+file]);
      exit(1);
      }
    variance_flag &= vflag;
    }
/*compute ratios of calibration coefficients to those of the zeroeth file*/
  for(channel = 0; channel != num_channels; channel++)
    {
    for(file = 1; file != argc-2; file++)
      calibration_ratio[file][channel] /= calibration_ratio[0][channel];
    calibration_ratio[0][channel] = 1.0;
    }
/*open the output*/
  if((out = fopen(argv[argc-1],
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    )) == NULL)
    {
    fprintf(stderr, "couldn't open output %s\n", argv[argc-1]);
    exit(errno);
    }
/*copy the header from the zeroeth input file into the output file*/
  rewind(*in);
  for(sample = 0; sample != packed_sizeof_SETUP; sample++)
    putc(getc(*in), out);
  for(channel = 0; channel != num_channels; channel++)
    {
    for(sample = 0; sample != EL_baseline_offset; sample++)
      putc(getc(*in), out);
    putc(baseline[channel]&0xff, out);
    putc((baseline[channel]>>8)&0xff, out);
    fseek(*in, 2L, SEEK_CUR);
    for(sample = EL_baseline_offset+2; sample != packed_sizeof_ELECTLOC; sample++)
      putc(getc(*in), out);
    }
/*use the calibration ratios to add the averages*/
  for(file = 0; file != argc-2; file++)
    fseek(in[file], (long)(packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC), SEEK_SET);
  for(channel = 0; channel != num_channels; channel++)
    {
    int i;
  /*copy the sweep header*/
    for(i = 0; i != 5; i++)
      putc(getc(*in), out);
    for(file = 1; file != argc-2; file++)
      fseek(in[file], 5L, SEEK_CUR);
    for(sample = 0; sample != num_samples; sample++)
      {
      float sum;
      sum = 0.0;
      for(file = 0; file != argc-2; file++)
	sum += read_Intel_float(in[file])*calibration_ratio[file][channel];
      write_Intel_float(sum, out);
      }
    }
  if(variance_flag)
  /*If the standard deviations were included at the end of the file, do the same
    as above, except no header.*/
    for(channel = 0; channel != num_channels; channel++)
      for(sample = 0; sample != num_samples; sample++)
	{
	float sum;
	sum = 0.0;
	for(file = 0; file != argc-2; file++)
	  sum += read_Intel_float(in[file])*calibration_ratio[file][channel];
	write_Intel_float(sum, out);
	}
  fclose(out);
  for(file = argc-3; file != -1; file--)
    fclose(in[file]);
  exit(0);
  }
