/*ztrans.c - z-transform a GnuroScan averaged EEG file, with variance and mean
  calculated from another (possibly identical) GnuroScan averaged EEG file.
  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).
*/

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

void main(argc, argv)
int argc;
char **argv;
  {
  register int sample, channel;
  register double d, s, ss;
  register float calibration;
  int num_samples, num_channels;
  double *mean, *stdev;
  FILE *avg, *zfile;
  char **ch_hdr, erp[packed_sizeof_SETUP];
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc != 4)
    {
    fprintf(stderr, "usage: %s <baseline.avg> <voltages.avg> <zscores.avg>\n", *argv);
    exit(1);
    }
/*open the baseline file*/
  if((avg = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[1]);
    exit(1);
    }
/*read the header from the baseline file*/
  if(fread(erp, 1, packed_sizeof_SETUP, avg) != packed_sizeof_SETUP)
    {
    perror(argv[1]);
    fclose(avg);
    exit(errno);
    }
  num_channels = S_nchannels(erp);
  num_samples = S_pnts(erp);
/*read the baseline channel headers*/
  ch_hdr = (char **)calloc(num_channels, sizeof(char *));
  for(channel = 0; channel != num_channels; channel++)
    {
    ch_hdr[channel] = (char *)malloc(packed_sizeof_ELECTLOC);
    if(fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg) != packed_sizeof_ELECTLOC)
      {
      perror(argv[1]);
      fclose(avg);
      exit(errno);
      }
    }
/*scan baseline data and compute mean and standard deviation for each channel*/
  mean = (double *)calloc(num_channels, sizeof(double));
  stdev = (double *)calloc(num_channels, sizeof(double));
  for(channel = 0; channel != num_channels; channel++)
    {
    calibration = EL_calib(ch_hdr[channel]);
    s = 0.0;
    ss = 0.0;
    fseek(avg, 5L, SEEK_CUR);
    for(sample = 0; sample != num_samples; sample++)
      {
      d = read_Intel_float(avg);
      s += d;
      ss += d*d;
      }
    mean[channel] = calibration*s/num_samples;
    stdev[channel] = calibration*sqrt((ss-s*s/num_samples)/(num_samples-1));
    }
  fclose(avg);
/*open the file whose data are to be z-transformed*/
  if((avg = fopen(argv[2],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[2]);
    exit(1);
    }
/*read the file header*/
  if(fread(erp, 1, packed_sizeof_SETUP, avg) != packed_sizeof_SETUP)
    {
    perror(argv[2]);
    fclose(avg);
    exit(errno);
    }
/*check that the number of channels is the same as in the baseline file*/
  if(S_nchannels(erp) != num_channels)
    {
    fprintf(stderr, "%d channels in %s incompatible with %d channels in %s\n", num_channels, argv[1], S_nchannels(erp), argv[2]);
    fclose(avg);
    exit(1);
    }
  num_samples = S_pnts(erp);
/*read the channel headers*/
  for(channel = 0; channel != num_channels; channel++)
    if(fread(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, avg) != packed_sizeof_ELECTLOC)
      {
      perror(argv[2]);
      fclose(avg);
      exit(errno);
      }
/*open the output file*/
  if((zfile = fopen(argv[3],
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    )) == NULL)
    {
    perror(argv[3]);
    fclose(avg);
    exit(errno);
    }
/*variance data will not be saved, so clear the variance flag*/
  S_variance(erp) = 0;
/*write the file header*/
  fwrite(erp, 1, packed_sizeof_SETUP, zfile);
/*write the channel headers*/
  for(channel = 0; channel != num_channels; channel++)
    {
  /*calibration has been multiplied in, so set the calibration factor to 1*/
    calibration = EL_calib(ch_hdr[channel]);
    put_EL_calib(ch_hdr[channel], 1.0);
    fwrite(ch_hdr[channel], 1, packed_sizeof_ELECTLOC, zfile);
  /*restore original calibration factor, for use later*/
    put_EL_calib(ch_hdr[channel], calibration);
    }
  for(channel = 0; channel != num_channels; channel++)
    {
  /*retrieve calibration factor*/
    calibration = EL_calib(ch_hdr[channel]);
  /*copy 5-byte trace header*/
    for(sample = 0; sample != 5; sample++)
      putc(getc(avg), zfile);
  /*z-transform each sample*/
    for(sample = 0; sample != num_samples; sample++)
      write_Intel_float((read_Intel_float(avg)*calibration-mean[channel])/stdev[channel], zfile);
    }
/*close files and release storage*/
  fclose(zfile);
  fclose(avg);
  free(stdev);
  free(mean);
  for(channel = num_channels-1; channel >= 0; channel--)
    free(ch_hdr[channel]);
  free(ch_hdr);
  }
