/*average.c - average a scored GnuroScan continuous file for the fastshft.c
  experiment.
  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 <sys/types.h>
#include <sys/stat.h>
/*#include <unistd.h> - def. of 'select' conflicts w/ nr.h*/ extern long fork(); extern int execl();
#include <string.h>
#ifdef MSDOS
#  define strcasecmp(s, t) strcmpi(s, t)
#endif
#include <math.h>
#include <setjmp.h>
#include "median.h"

extern void mrqmin(float x[], float y[], float sig[], int ndata, float a[],
	const int ia[], int ma, float**, float**, float*,
	void (*funcs)(double, float[], float*, float[], int), float*);
extern void drealft(double[], unsigned long, int);

#include <errno.h>
#include "fastshft.h"
#include "intel.h"
#include "gnuro.h"
#define PI 3.141592653589793

#define ST_INFINITY (1+ST_LONGEST_ISI/ST_MSECS_PER_TICK)
/*greater than the greatest possible number of ticks*/

#define LEFT 0
#define RIGHT 1

#define HIT 0
#define MISS 1
#define IGNORE 2

#define MAX_NUM_OBSERVATIONS 1023 /*max. # of observations at a single point*/
#define MAX_NUM_SAMPLES (1+2*(int)(((long)(ST_LONGEST_ISI))*AV_MAX_SRATE/1000L))

#ifndef HISTO_ONLY
#  define NUM_COMPONENTS 2	/* < AV_MAX_FFT_LENGTH/2 */
#else
#  define NUM_COMPONENTS 1
#endif

int behav_iti_bin_latency(behav_iti_bin)
int behav_iti_bin;
  {
  switch(behav_iti_bin)
    {
    case 0: return(ST_MSECS_PER_TICK/2);
    case 1: return(4*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    case 2: return(7*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    case 3: return(11*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    case 4: return(14*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    case 5: return(18*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    case 6: return(21*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    default: return(28*ST_MSECS_PER_TICK+ST_MSECS_PER_TICK/2);
    }
  }

jmp_buf env; int in_mrqmin=0;
FILE *evt, *cnt;
short num_channels;
int	num_samples,		/*number of samples in the epoch (always odd)*/
	stimulus_latency,	/*# of samples between event code & stimulus*/
	bytes_per_sample,	/*#of bytes in file for each sample*/
	window_length,		/*# of samples in a single stimulus cycle*/
	fft_length,		/*a power of two close to 'window_length'*/
	num_components,		/*# of freq. components saved (<=fft_length/2)*/
	heog_channel;		/*channel # of HEOG channel*/
double heog_scale;		/*conversion to uV for HEOG channel*/
long	prev_shift;		/*event table offset of last shift*/
char evtrec[packed_sizeof_EVENT2];
int is_eye_channel[CH_MAX];
double cosine_bell[AV_MAX_FFT_LENGTH];
long tsum[2][3][AV_NUM_BINS][CH_MAX][MAX_NUM_SAMPLES];
double t_ss[2][3][AV_NUM_BINS][CH_MAX][MAX_NUM_SAMPLES]; /*too big for a long*/
double fsum[2][3][AV_NUM_BINS][NUM_COMPONENTS][CH_MAX][MAX_NUM_SAMPLES], f_ss[2][3][AV_NUM_BINS][NUM_COMPONENTS][CH_MAX][MAX_NUM_SAMPLES];
short divisor[2][3][AV_NUM_BINS][MAX_NUM_SAMPLES];
float noise_correction_param[CH_MAX];
long onset_tsum[2][CH_MAX][MAX_NUM_SAMPLES];
double onset_t_ss[2][CH_MAX][MAX_NUM_SAMPLES];
double onset_fsum[2][NUM_COMPONENTS][CH_MAX][MAX_NUM_SAMPLES], onset_f_ss[2][NUM_COMPONENTS][CH_MAX][MAX_NUM_SAMPLES];
short onset_divisor[2][MAX_NUM_SAMPLES];
unsigned int histogram[2][2][(1+2*ST_LONGEST_ISI*AV_MAX_SRATE/1000)/2][1+2*AV_HISTO_MAXV];

/*Numerical Recipes error handler*/
void nrerror(mesg)
char *mesg;
  {
  fflush(stdout);
  fprintf(stderr, "%s\n", mesg);
  if(in_mrqmin)
    {
    fflush(stderr);
    longjmp(env, 1);
    }
  else
    exit(1);
  }

/*Initialise HEOG histograms.*/
void init_histograms(ch_hdr)
char (*ch_hdr)[packed_sizeof_ELECTLOC];
  {
  register int v, t, accuracy, laterality;
  for(heog_channel = 0; (heog_channel != num_channels) && (strcasecmp("heog", ch_hdr[heog_channel]) != 0); heog_channel++)
    ;
  if(heog_channel != num_channels)
    {
    heog_scale = EL_sensitivity(ch_hdr[heog_channel])*EL_calib(ch_hdr[heog_channel])/204.8;
    for(laterality = LEFT; laterality <= RIGHT; laterality++)
      for(accuracy = HIT; accuracy <= MISS; accuracy++)
	for(t = 0; t != num_samples/2; t++)
	  for(v = -AV_HISTO_MAXV; v <= AV_HISTO_MAXV; v++)
	    histogram[laterality][accuracy][t][AV_HISTO_MAXV+v] = 0;
    }
  }

/*Normalise HEOG histograms for easier viewing.*/
void normalise_histograms()
  {
  register unsigned int sum, elt;
  register int v, t, accuracy, laterality;
  for(laterality = LEFT; laterality <= RIGHT; laterality++)
    for(accuracy = HIT; accuracy <= MISS; accuracy++)
      for(t = 0; t != num_samples/2; t++)
	{
	sum = 0;
	for(v = -AV_HISTO_MAXV; v <= AV_HISTO_MAXV; v++)
	  sum += histogram[laterality][accuracy][t][AV_HISTO_MAXV+v];
	if(sum != 0)
	  for(v = -AV_HISTO_MAXV; v <= AV_HISTO_MAXV; v++)
	    {
	    elt = histogram[laterality][accuracy][t][AV_HISTO_MAXV+v]*10*255/sum;
	    histogram[laterality][accuracy][t][AV_HISTO_MAXV+v] = (elt>255)?255:elt;
	    }
	}
  }

/*Output histograms in simple binary row-major format, & in RLE if available.*/
void write_histograms()
  {
  static char frombin[] = "/usr/NeWS/bin/frombin", laterality_code[] = {'l', 'r'}, accuracy_code[] = {'h', 'm'};
  struct stat stbuf;
  register int v, t, accuracy, laterality;
  int sgi_flag;
  FILE *fp;
  char num_samples_ascii[6], num_vsteps_ascii[6], bin_img_file[13], rle_img_file[12];
  sgi_flag = (stat(frombin, &stbuf) == 0);
  strcpy(bin_img_file+2, "_histo.bin");
  for(laterality = LEFT; laterality <= RIGHT; laterality++)
    {
    bin_img_file[0] = laterality_code[laterality];
    for(accuracy = HIT; accuracy <= MISS; accuracy++)
      {
      bin_img_file[1] = accuracy_code[accuracy];
      if((fp = fopen(bin_img_file,
#ifdef MSDOS
	"wb"
#else
	"w"
#endif
	)) != NULL)
	for(v = -AV_HISTO_MAXV; v <= AV_HISTO_MAXV; v++)
	  for(t = 0; t != num_samples/2; t++)
	    putc(histogram[laterality][accuracy][t][AV_HISTO_MAXV+v], fp);
      fclose(fp);
      }
    }
  if(sgi_flag)
  /*translate binary image files to RLE files*/
    {
    fflush(stdout);
    fflush(stderr);
    sprintf(num_samples_ascii, "%d", num_samples/2);
    sprintf(num_vsteps_ascii, "%d", 1+2*AV_HISTO_MAXV);
    strcpy(rle_img_file, "rm_histo.bw");
    if(fork() == 0)
      execl(frombin, "frombin", bin_img_file, rle_img_file, num_samples_ascii, num_vsteps_ascii, (char *)0);
    bin_img_file[1] = rle_img_file[1] = 'h';
    if(fork() == 0)
      execl(frombin, "frombin", bin_img_file, rle_img_file, num_samples_ascii, num_vsteps_ascii, (char *)0);
    bin_img_file[0] = rle_img_file[0] = 'l';
    if(fork() == 0)
      execl(frombin, "frombin", bin_img_file, rle_img_file, num_samples_ascii, num_vsteps_ascii, (char *)0);
    bin_img_file[1] = rle_img_file[1] = 'm';
    execl(frombin, "frombin", bin_img_file, rle_img_file, num_samples_ascii, num_vsteps_ascii, (char *)0);
    exit(errno);
    }
  }

int nearest_power_of_two(x)
int x;
  {
  register int result;
  for(result = 1; result < x; result <<= 1)
    ;
  return(result-x<=x-result/2? result: result/2);
  }

/*Stretch or squeeze the circular buffer s[offset%ns..(offset+ns-1)%ns] onto
  d[0..nd-1].*/
void resample(s, ns, d, nd, offset)
short *s;
int ns;
double *d;
int nd, offset;
  {
  register int i, j, acc;
  if(ns == nd)
    for(i = 0; i != nd; i++)
      d[i] = (double)(s[(i+offset)%nd]);
  else
    for(i = 0; i != nd; i++)
      {
      acc = 0;
      for(j = i*ns; j != (i+1)*ns; j++)
	acc += s[(j/nd+offset)%ns];
      d[i] = ((double)acc)/ns;
      }
  }

/*Subtract any DC offset and linear trend (based upon the difference between
  the endpoints y[0] and y[n-1]) from the series y[0..n-1].*/
void detrend(y, n)
double *y;
int n;
  {
  register int i;
  register float decr, dy;
  dy = (y[n-1]-y[0])/(n-1);
  decr = y[0];
  for(i = 0; i != n; i++)
    {
    y[i] -= decr;
    decr += dy;
    }
  }


/*Convolve y[0..n-1] with x[0..n-1], and store the result in x[0..n-1].*/
void convolve(x, y, n)
double *x, *y;
int n;
  {
  register int i;
  for(i = 0; i != n; i++)
    x[i] *= y[i];
  }

/*Slave function to Numerical Recipes' Levenberg-Marquardt curve-fitting
  routine.  Evaluates the predicted local range

          1+tanh(n-p[1])      p[2]        1-tanh(n-p[1])    p[2]
  R(n) =  -------------- * ----------  +  -------------- * -------
                2          sqrt(p[1])           2          sqrt(n)

  and the partial derivatives wrt the two parameters,

 dR      p[2]    1+tanh(n-p[1])
-----(n)=---- * (-------------- + (sech(n-p[1]))^2 * (1/sqrt(n) - 1/sqrt(p[1])))
dp[1]     2        p[1]^(3/2)

  and

   dR         1      1+tanh(n-p[1])   1-tanh(n-p[1])
  -----(n) = --- * ( -------------- + -------------- )
  dp[2]       2        sqrt(p[1])        sqrt(n)

  In the definition of R(n) above, the left-hand factor of each term determines
  the relative contribution of the corresponding right-hand factor.  The two
  left-hand factors are complementary, having unit sum.  They are highly
  nonlinear functions that implement a rather abrupt, yet continuous (and so
  differentiable) switch between dominance of the left-hand term at high n and
  dominance of the right-hand term at low n.  When the value of n, the number
  of observations in the average, is high, it has little effect on the
  amplitude of the relatively low-frequency signal being measured.  Conversely,
  when n is low, noise biases the amplitude upward with an inverse dependency
  on the square root of n.*/
void funcs(n, p, R, partial, two)
double n;
/*must be double instead of float to prevent stack allignment problems due to
  default argument promotion - corresponding changes must be made to the
  Numerical Recipes files mrqmin.c and mrqcof.c.*/
float *p, *R, *partial;
int two;
  {
  register double squash, squashd, sqrt_p1, sqrt_n;
  squash = tanh((double)(n-p[1]));
  squashd = 1.0/cosh((double)(n-p[1]));
  sqrt_p1 = sqrt((double)(p[1]));
  sqrt_n = sqrt((double)n);
  *R = (float)(((1.0+squash)*p[2]/sqrt_p1 + (1.0-squash)*p[2]/sqrt_n)/2.0);
  partial[1] = (float)((p[2]/2.0)*((1.0+squash)/(p[1]*sqrt_p1) + squashd*squashd*(1.0/sqrt_n-1.0/sqrt_p1)));
  partial[2] = (float)(((1.0+squash)/sqrt_p1 + (1.0-squash)/sqrt_n)/2.0);
  }

/*Since the number of points in the average drops towards the ends of the epoch,
  the signal is noisier at the ends of the epoch than in the middle.  This
  must be compensated for, or else the moving-window amplitude calculated from
  the Fourier transform will be artificially bowed upwards at the ends of the
  epoch.  This routine models the local range of the measured signal as a
  nonlinear combination of a constant component due to the true signal and a
  noise component that varies inversely with the square root of the number of
  points in the average.  The range is calculated for each point in a moving
  window whose length is one complete stimulus period.  Levenberg-Marquardt
  error minimisation is then used to zero in on optimal values for the model's
  two parameters, described above.  This will allow the computation of a
  correction factor to eliminate the bias due to noise.*/
void squelch_noise()
  {
  register int num_obs, sample, num_points, num_iterations, channel, iti_bin, accuracy, laterality;
  register short local_range;
  register float prev_lambda;
  register btree *signal;
  float lambda, chi_squared,
	params[2],
	num_observations[MAX_NUM_OBSERVATIONS], /*abscissa*/
	range[MAX_NUM_OBSERVATIONS], /*ordinate*/
	sd_of_range[MAX_NUM_OBSERVATIONS],
	*(covar[2]), *(alpha[2]), storage[8]; /*work space*/
  long summed_range[1+MAX_NUM_OBSERVATIONS], summed_squared_range[1+MAX_NUM_OBSERVATIONS];
  short n[1+MAX_NUM_OBSERVATIONS];
  static const int var_flags[] = {1, 1};
/*initialise storage*/
  for(num_points = 0; num_points != 2; num_points++)
    {
    covar[num_points] = storage+num_points*2-1;
    alpha[num_points] = storage+num_points*2+3;
    }
/*form a noise correction for each non-eye channel separately*/
  for(channel = 0; channel != num_channels; channel++)
    if(!(is_eye_channel[channel]))
      {
    /*initialise statistics*/
      for(num_obs = 1; num_obs <= MAX_NUM_OBSERVATIONS; num_obs++)
	{
	summed_range[num_obs] = summed_squared_range[num_obs] = 0L;
	n[num_obs] = 0;
	}
    /*compute average local range as a function of number of observations*/
      for(laterality = LEFT; laterality <= RIGHT; laterality++)
	for(accuracy = HIT; accuracy <= IGNORE; accuracy++)
	  for(iti_bin = 0; iti_bin != AV_NUM_BINS; iti_bin++)
	    {
	  /*initialise a moving window with samples of the averaged signal*/
	    signal = create_btree(window_length);
	    for(sample = num_samples/2+window_length/2; (sample != num_samples/2-(window_length+1)/2+1) && (divisor[laterality][accuracy][iti_bin][sample] != 0); sample--)
	      insert_newest((short)(tsum[laterality][accuracy][iti_bin][channel][sample]/divisor[laterality][accuracy][iti_bin][sample]), signal);
	  /*'signal' contains the rightmost window_length-1 points in the
	    prestimulus time-domain average of the current channel*/
	   while((sample >= 0) && (divisor[laterality][accuracy][iti_bin][sample] != 0))
	    /*inv: 'signal' contains the averaged signal from points
	      sample+1 to sample+window_length-1, and sample_num, range, and
	      sd_of_range have been set for each point considered.*/
	      {
	      insert_newest((short)(tsum[laterality][accuracy][iti_bin][channel][sample]/divisor[laterality][accuracy][iti_bin][sample]), signal);
	      local_range = extract_maximum(signal)-extract_minimum(signal);
	      delete_oldest(signal);
	      if(divisor[laterality][accuracy][iti_bin][sample] <= MAX_NUM_OBSERVATIONS)
		{
		summed_range[divisor[laterality][accuracy][iti_bin][sample]] += local_range;
		summed_squared_range[divisor[laterality][accuracy][iti_bin][sample]] += local_range*(long)local_range;
		n[divisor[laterality][accuracy][iti_bin][sample--]]++;
		}
	      }
	    destroy_btree(signal);
	    }
    /*'summed_range', 'summed_squared_range', and 'n' have been accumulated from
      each sample across all ITI bins, accuracies, and lateralities.  Now
      extract average values and standard deviations for the local range at each
      sample.*/
      num_points = 0;
    /*extract average local ranges for each number of observations*/
      for(num_obs = MAX_NUM_OBSERVATIONS; num_obs != 0; num_obs--)
	if(n[num_obs] > 1)
	  {
	  num_observations[num_points] = (float)num_obs;
	  range[num_points] = ((float)(summed_range[num_obs]))/n[num_obs];
	  sd_of_range[num_points] = (summed_squared_range[num_obs]-((float)(summed_range[num_obs]*summed_range[num_obs]))/n[num_obs])/(n[num_obs]-1);
	/*force SD estimate to be something reasonable*/
	  if(sd_of_range[num_points] < range[num_points])
	    sd_of_range[num_points] = range[num_points];
	  num_points++;
	  }
/*DEBUG*/ for(sample = 0; sample != num_points; sample++) printf("%.0f %f %f\n",num_observations[sample],range[sample],sd_of_range[sample]);
    /*initialise Levenberg-Marquardt parameters*/
      params[0] = 100.0;	   /*switch at ~100 observations*/
      params[1] = range[0]*10.0; /*~minimum range times sqrt(params[0])*/
      lambda = -1.0;
    /*do at least 10 iterations*/
      num_iterations = 10;
    /*flag the error handler that we're inside mrqmin*/
      in_mrqmin = 1;
    /*if a singular matrix arises, use a default cutoff of zero*/
      if(setjmp(env))
	noise_correction_param[channel] = 0.0;
      else
	{
	do
	  {
	  prev_lambda = lambda;
	  mrqmin(num_observations-1, range-1, sd_of_range-1, num_points, params-1, var_flags-1, 2, covar-1, alpha-1, &chi_squared, funcs, &lambda);
	  }
	while((--num_iterations > 0) || (lambda < prev_lambda));
      /*clean up*/
	lambda = 0.0;
	mrqmin(num_observations-1, range-1, sd_of_range-1, num_points, params-1, var_flags-1, 2, covar-1, alpha-1, &chi_squared, funcs, &lambda);
	noise_correction_param[channel] = params[0];
	}
    /*clear the mrqmin flag*/
      in_mrqmin = 0;
      printf("channel %d squelch cutoff at %d observations\n", channel, (int)(0.5+noise_correction_param[channel]));
      }
  }

/*Compute a noise correction factor for a given number of observations.
  'ncrit' here is params[0] from squelch_noise().*/
double find_noise_correction(ncrit, n)
double ncrit, n;
  {
  register double squash;
  squash = tanh(n-ncrit);
  return(2.0/(1.0+squash+(1.0-squash)*sqrt(ncrit/n)));
  }

/*Smooth the given vector of noise correction factors.*/
void smooth_noise_corrections(noise_correction)
double *noise_correction;
  {
  register double summed_noise_corrections;
  register int sample;
  summed_noise_corrections = 0.0;
/*initialise*/
  for(sample = 0; sample != window_length/2; sample++)
    summed_noise_corrections += noise_correction[sample];
/*smooth left boundary*/
  while(sample != window_length-1)
    {
    summed_noise_corrections += noise_correction[sample];
    noise_correction[sample-window_length/2] = summed_noise_corrections/(1+sample);
    sample++;
    }
/*smooth central interval*/
  while(sample != num_samples)
  /*inv: noise_correction[0 .. sample-window_length/2-1] has been smoothed*/
    {
    summed_noise_corrections += noise_correction[sample];
    noise_correction[sample-window_length/2] = summed_noise_corrections/window_length;
    summed_noise_corrections -= noise_correction[++sample-window_length];
    }
/*smooth right boundary*/
  while(sample != num_samples+window_length/2)
  /*inv: noise_correction[0 .. sample-window_length/2-1] has been smoothed*/
    {
    noise_correction[sample-window_length/2] = summed_noise_corrections/(num_samples+window_length-sample);
    summed_noise_corrections -= noise_correction[++sample-window_length];
    }
  }

/*Write the given average to a GnuroScan file whose name is as given and whose
  extension is '.avg'.  Use the given header information.  If 'lsum' is
  non-null, then use its entries to generate the average and the entries of
  'l_ss' to generate the variance.  Otherwise, use the entries of 'dsum' and
  'd_ss'.*/
void write_avg(name, erp_hdr, ch_hdr, lsum, l_ss, dsum, d_ss, divisor, num_samples, num_channels)
char *name, *erp_hdr, (*ch_hdr)[packed_sizeof_ELECTLOC];
long (*lsum)[MAX_NUM_SAMPLES];
double (*l_ss)[MAX_NUM_SAMPLES], (*dsum)[MAX_NUM_SAMPLES], (*d_ss)[MAX_NUM_SAMPLES];
short *divisor;
int num_samples, num_channels;
  {
  register int channel, sample;
  float eeg, min_eeg, max_eeg, epoch[MAX_NUM_SAMPLES];
  FILE *avg;
  char filename[13];
  double noise_correction[MAX_NUM_SAMPLES];
  sprintf(filename, "%s.avg", name);
  avg = fopen(filename,
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    );
  fwrite(erp_hdr, packed_sizeof_SETUP, 1, avg);
  for(channel = 0; channel != num_channels; channel++)
    fwrite(ch_hdr[channel], packed_sizeof_ELECTLOC, 1, avg);
  for(channel = 0; channel != num_channels; channel++)
    {
  /*compute a vector of noise correction factors to be convolved with the
    time-domain averages of all non-eye channels*/
    if(!(is_eye_channel[channel]))
      {
      noise_correction[num_samples/2] = find_noise_correction((double)(noise_correction_param[channel]), (double)(divisor[num_samples/2]));
      for(sample = 1+num_samples/2; sample < num_samples; sample++)
	noise_correction[sample] = (divisor[sample]? find_noise_correction((double)(noise_correction_param[channel]), (double)(divisor[sample])): noise_correction[sample-1]);
      for(sample = num_samples/2-1; sample >= 0; sample--)
	noise_correction[sample] = (divisor[sample]? find_noise_correction((double)(noise_correction_param[channel]), (double)(divisor[sample])): noise_correction[sample+1]);
      smooth_noise_corrections(noise_correction);
      }
  /*5-byte dummy header - no longer used*/
    for(sample = 0; sample != 5; sample++)
      putc(0, avg);
  /*vectors of floats, in channel-major order*/
    min_eeg = 1e+10;
    max_eeg = -1e+10;
    for(sample = 0; sample != num_samples; sample++)
      {
      epoch[sample] = EL_sensitivity(ch_hdr[channel])/204.8*
	(divisor[sample]?
	  (is_eye_channel[channel]?
	    (float)(((double)(lsum[channel][sample]))/divisor[sample]-EL_baseline(ch_hdr[channel])):
	    (dsum == (double (*)[MAX_NUM_SAMPLES])0)?
	      (float)(noise_correction[sample]*((double)(lsum[channel][sample]))/divisor[sample]-EL_baseline(ch_hdr[channel])):
	    (float)(dsum[channel][sample]/divisor[sample])):
	  0.0);
      if(noise_correction[sample] > 0.5)
	{
	if(epoch[sample] > max_eeg)
	  max_eeg = epoch[sample];
	if(epoch[sample] < min_eeg)
	  min_eeg = epoch[sample];
	}
      }
    if(min_eeg == 1e+10)
      min_eeg = -1e+10;
    if(max_eeg == -1e+10)
      max_eeg = 1e+10;
  /*clip high amplitudes in intervals that are noisy (as indicated by a low
    squelch factor), and write the resulting vector to disk*/
    for(sample = 0; sample != num_samples; sample++)
      {
      if(noise_correction[sample] < 0.5)
	{
	if(epoch[sample] > max_eeg)
	  epoch[sample] = max_eeg;
	else if(epoch[sample] < min_eeg)
	  epoch[sample] = min_eeg;
	}
      store_Intel_float((char *)&eeg, 0, epoch[sample]);
      fwrite(&eeg, sizeof(eeg), 1, avg);
      }
    }
/*same as above but generate standard deviation instead of average*/
  for(channel = 0; channel != num_channels; channel++)
    for(sample = 0; sample != num_samples; sample++)
      {
      store_Intel_float((char *)&eeg,
	0,
	(divisor[sample]>1)?
	  (is_eye_channel[channel]?
	    (float)(sqrt((l_ss[channel][sample] - (double)(lsum[channel][sample])*(double)(lsum[channel][sample])/divisor[sample])/(divisor[sample]-1))):
	    (dsum == (double (*)[MAX_NUM_SAMPLES])0)?
	      (float)(noise_correction[sample]*(sqrt((l_ss[channel][sample] - (double)(lsum[channel][sample])*(double)(lsum[channel][sample])/divisor[sample])/(divisor[sample]-1)))):
	    (float)(sqrt((d_ss[channel][sample] - dsum[channel][sample]*dsum[channel][sample]/divisor[sample])/(divisor[sample]-1)))):
	  0.0);
      fwrite(&eeg, sizeof(eeg), 1, avg);
      }
  fclose(avg);
  }

/*Write the given vector of sample sizes to a text file whose name is as given
  and whose extension is '.txt'.*/
void write_sample_sizes(name, divisor, num_samples)
char *name;
short *divisor;
int num_samples;
  {
  register int sample;
  FILE *vec;
  char filename[13];
  sprintf(filename, "%s.txt", name);
  vec = fopen(filename, "w");
  for(sample = 0; sample != num_samples; sample++)
    fprintf(vec, "%d\n", divisor[sample]);
  fclose(vec);
  }

/*Add into the running sums 'tsum' (time-domain) and 'fsum' (frequency-domain)
  and the running sums of squares 't_ss' and 'f_ss' the epoch surrounding the
  event whose record is in 'evtrec', and update the given HEOG distribution
  histogram.*/
void average(tsum, t_ss, fsum, f_ss, divisor, heog_histogram)
long tsum[CH_MAX][MAX_NUM_SAMPLES];
double t_ss[CH_MAX][MAX_NUM_SAMPLES], fsum[AV_MAX_FFT_LENGTH/2][CH_MAX][MAX_NUM_SAMPLES], f_ss[AV_MAX_FFT_LENGTH/2][CH_MAX][MAX_NUM_SAMPLES];
short *divisor;
unsigned int (*heog_histogram)[1+2*AV_HISTO_MAXV];
  {
  register int sample, channel, accept, prev_accept, window_start, i;
  register long offset, prev_offset, offset_to_stimulus;
  register double f_elt;
  register long t_elt;
  register short heog_baseline;
  long savedpos;
  short window[CH_MAX][ST_MSECS_PER_TICK*AV_MAX_SRATE/1000];
  double transform[AV_MAX_FFT_LENGTH];
  if(E_Accept(evtrec))
    {
    offset_to_stimulus = E_Offset(evtrec);
  /*remember the current position in the event table*/
    savedpos = ftell(evt);
    do
      {
      fseek(evt, -2L*packed_sizeof_EVENT2, SEEK_CUR);
      fread(evtrec, packed_sizeof_EVENT2, 1, evt);
      } while((E_Accept(evtrec) || (E_Code(evtrec) != CC_SACCADE))
	 && (ftell(evt) != prev_shift)
	 && ((offset_to_stimulus - E_Offset(evtrec))/bytes_per_sample <= num_samples/2-stimulus_latency));
    if(!E_Accept(evtrec) && (E_Code(evtrec) == CC_SACCADE))
      fread(evtrec, packed_sizeof_EVENT2, 1, evt);
  /*'evt' contains the event that bounds the end of this epoch (i.e., the
    saccade preceding the current stimulus, the stimulus that evoked the most
    recent shift, or the event just beyond the end of the averaging window,
    whichever is latest).  The event following this is about to be read.*/
    offset = E_Offset(evtrec);
    accept = E_Accept(evtrec);
  /*make sure the initial offset lies within the averaging window*/
    if((offset_to_stimulus - offset)/bytes_per_sample > num_samples/2-stimulus_latency)
      offset = offset_to_stimulus - (num_samples/2-stimulus_latency)*bytes_per_sample;
    fseek(cnt, offset_to_stimulus+stimulus_latency*bytes_per_sample+heog_channel*sizeof(short), SEEK_SET);
    heog_baseline = read_Intel_short(cnt);
    do
    /*inv: 'accept' and 'offset' are set for the most recent event processed*/
      {
      int items_read;
      prev_offset = offset;
      prev_accept = accept;
      items_read = fread(evtrec, packed_sizeof_EVENT2, 1, evt);
      offset = E_Offset(evtrec);
      accept = (items_read==1) && E_Accept(evtrec);
      if((offset - offset_to_stimulus)/bytes_per_sample > 1+num_samples/2+stimulus_latency)
	offset = offset_to_stimulus + (1+num_samples/2+stimulus_latency)*bytes_per_sample;
    /*if interval between events hasn't been rejected, include it in the sum*/
      if(prev_accept && accept)
	{
	fseek(cnt, prev_offset-(window_length/2)*bytes_per_sample, SEEK_SET);
	window_start = 0;
	for(sample = 0; sample != window_length-1; sample++)
	  for(channel = 0; channel != num_channels; channel++)
	    window[channel][sample] = read_Intel_short(cnt);
	for(sample = num_samples/2 - stimulus_latency - (offset_to_stimulus-prev_offset)/bytes_per_sample; sample != num_samples/2 - stimulus_latency - (offset_to_stimulus-offset)/bytes_per_sample; sample++)
	  {
	/*update HEOG histogram*/
	  if((heog_histogram != (unsigned int (*)[1+2*AV_HISTO_MAXV])0) && (sample >= num_samples/2) && (sample < num_samples))
	    {
	    i = (int)rint((window[heog_channel][(window_start+window_length/2)%window_length]-heog_baseline)*heog_scale);
	    if(i > AV_HISTO_MAXV)
	      i = AV_HISTO_MAXV;
	    else if(i < -AV_HISTO_MAXV)
	      i = -AV_HISTO_MAXV;
	    heog_histogram[sample-num_samples/2][AV_HISTO_MAXV+i]++;
	    }
	  for(channel = 0; channel != num_channels; channel++)
	    {
	  /*update the window*/
	    window[channel][(window_start+window_length-1)%window_length] = read_Intel_short(cnt);
#ifndef HISTO_ONLY
	  /*include the window's new centre-point in the time-domain average*/
	    t_elt = window[channel][(window_start+window_length/2)%window_length];
	    tsum[channel][sample] += t_elt;
	    t_ss[channel][sample] += t_elt*t_elt;
	  /*resample onto a window whose length is a power of 2*/
	    resample(window[channel], window_length, transform, fft_length, window_start);
	  /*detrend below base frequency, to prevent spectral leakage*/
	    detrend(transform, fft_length);
	  /*convolve a cosine-bell window, to prevent spectral leakage*/
	    convolve(transform, cosine_bell, fft_length);
	  /*FFT*/
	    drealft(transform-1, fft_length, 1);
	  /*extract components of interest*/
	    for(i = 0; i != (num_components==fft_length/2? num_components-1: num_components); i++)
	      {
	      f_elt = transform[2*i+2]*transform[2*i+2]+transform[2*i+3]*transform[2*i+3];
	    /*compute RMS voltage*/
	      fsum[i][channel][sample] += sqrt(2.0*f_elt)/fft_length;
	      f_ss[i][channel][sample] += 2.0*f_elt/(fft_length*fft_length);
	      }
	    if(num_components == fft_length/2)
	    /*Nyquist term has no imaginary part - so if num_components is set
	      to the maximum possible number of components, we have to do the
	      Nyquist term as a special case.*/
	      {
	      f_elt = transform[1]*transform[1];
	      fsum[num_components-1][channel][sample] += sqrt(2.0*f_elt)/fft_length;
	      f_ss[num_components-1][channel][sample] += 2.0*f_elt/(fft_length*fft_length);
	      }
#endif
	    }
	  divisor[sample]++;
	  window_start++;
	  }
        }
      }
    while(
      !((offset > offset_to_stimulus) && ((E_StimType(evtrec) == 2) || (E_StimType(evtrec) == 4)) && (E_EpochEvent(evtrec) == 1) && E_Accuracy(evtrec)) /*next attended target has not yet occurred*/
       && (E_Accept(evtrec) || (E_Code(evtrec) != CC_SACCADE)) /*saccade has not occurred*/
       && (offset < offset_to_stimulus + (1+num_samples/2+stimulus_latency)*bytes_per_sample) /*maximum epoch length has not been reached*/
       && (prev_offset != offset) /*EOF has not been reached in event table*/ );
    fseek(evt, savedpos, SEEK_SET);
    }
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int channel, sample, iti, iti_bin, behav_iti_bin, frequency_bin, accuracy, laterality, count;
  register long event;
  short srate;
  long num_events;
  char filename[12];
  float half_epoch;	/*one half of epoch length, in seconds*/
/*hit_count[laterality][iti_bin] <= event_count[laterality][HIT][iti_bin],
  because some hits have no response latency information recorded and therefore
  are omitted from the calculation of average response latency which is based
  on the value of hit_count[][].*/
  int event_count[2][3][2*AV_NUM_BINS], hit_count[2][2*AV_NUM_BINS];
  double summed_response_latencies[2][2*AV_NUM_BINS], summed_squared_response_latencies[2][2*AV_NUM_BINS];
/*must define these as char arrays since IRIX cc doesn't pack structures*/
  char erp[packed_sizeof_SETUP];
  char channel_header[CH_MAX][packed_sizeof_ELECTLOC];
  char evt_tbl_hdr[packed_sizeof_TEEG];
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc != 2)
    {
    fprintf(stderr, "Usage: %s <file.cnt>\n", *argv);
    exit(1);
    }
  if((evt = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
  fread(&erp, sizeof(erp), 1, evt);
  num_channels = S_nchannels(erp);
  srate = S_rate(erp);
  num_samples = 1+2*(int)(((long)(ST_LONGEST_ISI))*srate/1000L);
  if((num_channels > CH_MAX) || (srate > AV_MAX_SRATE))
    {
    fclose(evt);
    fprintf(stderr, "Recompile %s with", *argv);
    if(num_channels > CH_MAX)
      fprintf(stderr, " CH_MAX=%d", num_channels);
    if(srate > AV_MAX_SRATE)
      fprintf(stderr, " AV_MAX_SRATE=%d", srate);
    putc('\n', stderr);
    exit(1);
    }
  stimulus_latency = (int)(ST_LATENCY/1000.0*srate);
  bytes_per_sample = num_channels*sizeof(short);
  window_length = (ST_MSECS_PER_TICK*srate+500)/1000;
  fft_length = nearest_power_of_two(window_length);
  num_components = fft_length < NUM_COMPONENTS? fft_length: NUM_COMPONENTS;
  for(sample = 1; sample != fft_length-1; sample++)
    cosine_bell[sample] = (cos(2.0*PI/(fft_length-1))-cos(2.0*PI*sample/(fft_length-1)))/(2.0*cos(2.0*PI/(fft_length-1)));
  cosine_bell[0] = cosine_bell[fft_length-1] = 0.0;
  for(channel = 0; channel != num_channels; channel++)
    {
    fread(channel_header+channel, sizeof(*channel_header), 1, evt);
    is_eye_channel[channel] = ((strcasecmp("heog", EL_lab(channel_header[channel])) == 0) || (strcasecmp("blink", EL_lab(channel_header[channel])) == 0) || (strcasecmp("upe", EL_lab(channel_header[channel])) == 0) || (strcasecmp("veog", EL_lab(channel_header[channel])) == 0));
    }
  init_histograms(channel_header);
  fseek(evt, S_EventTablePos(erp), SEEK_SET);
  fread(evt_tbl_hdr, packed_sizeof_TEEG, 1, evt);
  if(T_Teeg(evt_tbl_hdr) != 2)
    {
    fclose(evt);
    fprintf(stderr, "You must first score %s\n", argv[1]);
    exit(1);
    }
  fseek(evt, T_Offset(evt_tbl_hdr), SEEK_CUR);
#ifndef HISTO_ONLY
/*initialise arrays*/
  for(laterality = LEFT; laterality <= RIGHT; laterality++)
    for(iti_bin = 0; iti_bin != AV_NUM_BINS; iti_bin++)
      {
      summed_response_latencies[laterality][2*iti_bin] = summed_squared_response_latencies[laterality][2*iti_bin] = 0.0;
      summed_response_latencies[laterality][2*iti_bin+1] = summed_squared_response_latencies[laterality][2*iti_bin+1] = 0.0;
      hit_count[laterality][2*iti_bin] = 0;
      hit_count[laterality][2*iti_bin+1] = 0;
      for(accuracy = HIT; accuracy <= IGNORE; accuracy++)
	{
	event_count[laterality][accuracy][2*iti_bin] = 0;
	event_count[laterality][accuracy][2*iti_bin+1] = 0;
	for(sample = 0; sample != num_samples; sample++)
	  {
	  for(channel = 0; channel != num_channels; channel++)
	    {
	    tsum[laterality][accuracy][iti_bin][channel][sample] = 0L;
	    t_ss[laterality][accuracy][iti_bin][channel][sample] = 0.0;
	    for(frequency_bin = 0; frequency_bin != num_components; frequency_bin++)
	      {
	      fsum[laterality][accuracy][iti_bin][frequency_bin][channel][sample] = 0.0;
	      f_ss[laterality][accuracy][iti_bin][frequency_bin][channel][sample] = 0.0;
	      }
	    }
	  divisor[laterality][accuracy][iti_bin][sample] = 0;
	  }
	}
      }
  for(laterality = LEFT; laterality <= RIGHT; laterality++)
    for(sample = 0; sample != num_samples; sample++)
      {
      for(channel = 0; channel != num_channels; channel++)
	{
        onset_tsum[laterality][channel][sample] = 0L;
        onset_t_ss[laterality][channel][sample] = 0.0;
	for(frequency_bin = 0; frequency_bin != num_components; frequency_bin++)
	  {
          onset_fsum[laterality][frequency_bin][channel][sample] = 0.0;
          onset_f_ss[laterality][frequency_bin][channel][sample] = 0.0;
	  }
	}
      onset_divisor[laterality][sample] = 0;
      }
#endif

/*open 2nd pointer into file - eases reading EEG & event data simultaneously*/
  cnt = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    );
/*process events*/
  num_events = T_Size(evt_tbl_hdr) / packed_sizeof_EVENT2;
  event = 0L;
  do
    {
    fread(evtrec, packed_sizeof_EVENT2, 1, evt);
    event++;
    }
  while((E_StimType(evtrec) != EV_START) && (E_StimType(evtrec) != EV_START+1));
  laterality = E_StimType(evtrec) - EV_START;
  prev_shift = ftell(evt);
  iti = -1;
  fread(evtrec, packed_sizeof_EVENT2, 1, evt);
  event++;
  average(onset_tsum[laterality], onset_t_ss[laterality], onset_fsum[laterality], onset_f_ss[laterality], onset_divisor[laterality], (unsigned int (*)[1+2*AV_HISTO_MAXV])0);
  while(event < num_events)
  /*inv: 'event' events have been processed, and 'evtrec' contains the record
    for the next event to be processed.*/
    {
    fread(evtrec, packed_sizeof_EVENT2, 1, evt);
    event++;
    count = E_StimType(evtrec);
    iti += (((iti+1)/2!=ST_INFINITY) && (count<4));
    if((count == 2) || (count == 4))
      {
      laterality = (E_StimType(evtrec)==2? LEFT: RIGHT);
      accuracy = (E_EpochEvent(evtrec)==1? E_Accuracy(evtrec)? HIT: MISS: IGNORE);
      iti_bin = ((iti+1)/2)*(AV_NUM_BINS-1)/ST_INFINITY;
    /*behav_iti_bin is a finer-grained division than iti_bin --
      the former is used for behavioural measures; the latter is used for
      ERP averages, which are quite sensitive to noise.
      If you want to make the bins for the behavioural responses the same as
      the bins for the ERPs, then replace the following statement with
      'behav_iti_bin = iti_bin' and update the readout routine accordingly.*/
      behav_iti_bin = ((iti+1)/2==ST_INFINITY)? 2*AV_NUM_BINS-2: (iti-iti/(2*AV_NUM_BINS-2))/(2*AV_NUM_BINS-2);
    /*skip processing of unattended targets that have been flagged as possible
      sources of false-alarm responses*/
      if(E_EpochEvent(evtrec) != 2)
	{
        average(tsum[laterality][accuracy][iti_bin], t_ss[laterality][accuracy][iti_bin], fsum[laterality][accuracy][iti_bin], f_ss[laterality][accuracy][iti_bin], divisor[laterality][accuracy][iti_bin], (((iti_bin==0)||(accuracy==IGNORE))? ((unsigned int (*)[1+2*AV_HISTO_MAXV])0): histogram[laterality][accuracy]));
        event_count[laterality][accuracy][behav_iti_bin]++;
	if(accuracy == HIT)
	  {
	  float latency;
	  latency = E_Latency(evtrec);
	  if(latency != 0.0)
	    {
	    summed_response_latencies[laterality][behav_iti_bin] += latency;
	    summed_squared_response_latencies[laterality][behav_iti_bin] += latency*latency;
	    hit_count[laterality][behav_iti_bin]++;
	    }
	  prev_shift = ftell(evt);
	  iti = -1;
	  }
	}
      }
    }
/*close both pointers into the input file*/
  fclose(cnt);
  fclose(evt);
/*generate output*/
#ifndef HISTO_ONLY
  for(behav_iti_bin = 0; behav_iti_bin != 2*AV_NUM_BINS-1; behav_iti_bin++)
    for(laterality = LEFT; laterality <= RIGHT; laterality++)
      {
      for(accuracy = HIT; accuracy <= IGNORE; accuracy++)
	printf("%d %s %s in ITI bin %d (%dms .. %dms)\n", event_count[laterality][accuracy][behav_iti_bin], (laterality==LEFT? "left": "right"), (accuracy==HIT? "hits": accuracy==MISS? "misses": "ignores"), behav_iti_bin, behav_iti_bin_latency(behav_iti_bin), behav_iti_bin_latency(1+behav_iti_bin)-ST_MSECS_PER_TICK);
      printf("\thit ratio %f\n\tresponse latency average=%fms SD=%fms\n", event_count[laterality][HIT][behav_iti_bin]/(float)(event_count[laterality][HIT][behav_iti_bin]+event_count[laterality][MISS][behav_iti_bin]), 1000.0*(summed_response_latencies[laterality][behav_iti_bin]/hit_count[laterality][behav_iti_bin]), 1000.0*sqrt((summed_squared_response_latencies[laterality][behav_iti_bin] - summed_response_latencies[laterality][behav_iti_bin]*summed_response_latencies[laterality][behav_iti_bin]/hit_count[laterality][behav_iti_bin])/hit_count[laterality][behav_iti_bin]));
      }
  squelch_noise();
  half_epoch = num_samples/(float)srate/2.0;
  S_type(erp) = 1;		/*averaged*/
  S_variance(erp) = 1;		/*standard deviation data included*/
  put_S_pnts(erp, num_samples);
  put_S_rejstart(erp, -half_epoch);
  put_S_rejstop(erp, half_epoch);
  put_S_xmin(erp, -half_epoch);
  put_S_xmax(erp, half_epoch);
  put_S_AutoMin(erp, -half_epoch);
  put_S_AutoMax(erp, half_epoch);
  put_S_DisplayXmin(erp, -half_epoch);
  put_S_DisplayXmax(erp, half_epoch);
  for(channel = 0; channel != num_channels; channel++)
    put_EL_n(channel_header[channel], 1);
  for(laterality = LEFT; laterality <= RIGHT; laterality++)
    {
    *filename = (laterality==LEFT? 'l': 'r');
    for(accuracy = HIT; accuracy <= IGNORE; accuracy++)
      {
      filename[1] = (accuracy==HIT? 'h': accuracy==MISS? 'm': 'i');
      for(iti_bin = 0; iti_bin != AV_NUM_BINS; iti_bin++)
	{
	sprintf(filename+2, "%dn", iti_bin);
	write_sample_sizes(filename, divisor[laterality][accuracy][iti_bin], num_samples);
	filename[3] = 't';
	write_avg(filename, erp, channel_header, tsum[laterality][accuracy][iti_bin], t_ss[laterality][accuracy][iti_bin], (double (*)[MAX_NUM_SAMPLES])0, (double (*)[MAX_NUM_SAMPLES])0, divisor[laterality][accuracy][iti_bin], num_samples, num_channels);
	filename[3] = 'f';
	for(frequency_bin = 0; frequency_bin != num_components; frequency_bin++)
	  {
	  sprintf(filename+4, "%03d", (1000*(1+frequency_bin)+(ST_MSECS_PER_TICK/2))/ST_MSECS_PER_TICK);
	  write_avg(filename, erp, channel_header, tsum[laterality][accuracy][iti_bin], t_ss[laterality][accuracy][iti_bin], fsum[laterality][accuracy][iti_bin][frequency_bin], f_ss[laterality][accuracy][iti_bin][frequency_bin], divisor[laterality][accuracy][iti_bin], num_samples, num_channels);
	  }
	}
    /*now do the same as above for averages across all ITI bins*/
      filename[3] = '\0';
      for(iti_bin = 1; iti_bin < AV_NUM_BINS; iti_bin++)
	for(sample = 0; sample != num_samples; sample++)
	  {
	  for(channel = 0; channel != num_channels; channel++)
	    {
	    tsum[laterality][accuracy][0][channel][sample] += tsum[laterality][accuracy][iti_bin][channel][sample];
	    t_ss[laterality][accuracy][0][channel][sample] += t_ss[laterality][accuracy][iti_bin][channel][sample];
	    for(frequency_bin = 0; frequency_bin != num_components; frequency_bin++)
	      {
	      fsum[laterality][accuracy][0][frequency_bin][channel][sample] += fsum[laterality][accuracy][iti_bin][frequency_bin][channel][sample];
	      f_ss[laterality][accuracy][0][frequency_bin][channel][sample] += f_ss[laterality][accuracy][iti_bin][frequency_bin][channel][sample];
	      }
	    }
	  divisor[laterality][accuracy][0][sample] += divisor[laterality][accuracy][iti_bin][sample];
	  }
      filename[2] = '_';
      filename[3] = 'n';
      filename[4] = '\0';
      write_sample_sizes(filename, divisor[laterality][accuracy][0], num_samples);
      filename[3] = 't';
      write_avg(filename, erp, channel_header, tsum[laterality][accuracy][0], t_ss[laterality][accuracy][0], (double (*)[MAX_NUM_SAMPLES])0, (double (*)[MAX_NUM_SAMPLES])0, divisor[laterality][accuracy][0], num_samples, num_channels);
      filename[3] = 'f';
      for(frequency_bin = 0; frequency_bin != num_components; frequency_bin++)
	{
	sprintf(filename+4, "%03d", (1000*(1+frequency_bin)+(ST_MSECS_PER_TICK/2))/ST_MSECS_PER_TICK);
	write_avg(filename, erp, channel_header, tsum[laterality][accuracy][0], t_ss[laterality][accuracy][0], fsum[laterality][accuracy][0][frequency_bin], f_ss[laterality][accuracy][0][frequency_bin], divisor[laterality][accuracy][0], num_samples, num_channels);
	}
      }
    }
/*finally, store the averages of the trial onset periods, in order to get a
  picture of the buildup of the resonance*/
  put_S_rejstart(erp, 0.0);
  put_S_rejstop(erp, 2.0*half_epoch);
  put_S_xmin(erp, 0.0);
  put_S_xmax(erp, 2.0*half_epoch);
  put_S_AutoMin(erp, 0.0);
  put_S_AutoMax(erp, 2.0*half_epoch);
  put_S_DisplayXmin(erp, 0.0);
  put_S_DisplayXmax(erp, 2.0*half_epoch);
  filename[1] = 'o';
  filename[2] = '_';
  for(laterality = LEFT; laterality <= RIGHT; laterality++)
    {
    *filename = (laterality==LEFT? 'l': 'r');
    filename[3] = 't';
    write_avg(filename, erp, channel_header, onset_tsum[laterality], onset_t_ss[laterality], (double (*)[MAX_NUM_SAMPLES])0, (double (*)[MAX_NUM_SAMPLES])0, onset_divisor[laterality], num_samples, num_channels);
    filename[3] = 'f';
    for(frequency_bin = 0; frequency_bin != num_components; frequency_bin++)
      {
      sprintf(filename+4, "%03d", (1000*(1+frequency_bin)+(ST_MSECS_PER_TICK/2))/ST_MSECS_PER_TICK);
      write_avg(filename, erp, channel_header, onset_tsum[laterality], onset_t_ss[laterality], onset_fsum[laterality][frequency_bin], onset_f_ss[laterality][frequency_bin], onset_divisor[laterality], num_samples, num_channels);
      }
    }
#endif
  normalise_histograms();
  write_histograms();
  exit(0);
  }
