/*scanepl.c - frequency-domain EEG filter for vertical eye artefacts.
  Copyright (c) 1994 by Matthew Belmonte.

  The following files are included in this distribution:
	scanepl.c
	bigdata.c
	getopt.c
	scanepl.h
  If your system includes the getopt() system call, you do not need to
  compile and link the code in getopt.c.

  Licence is hereby granted to any and all individuals or organisations to
  use and to modify this software for the purpose of nonprofit research only.
  Use of this software in conjunction or association with any for-profit
  endeavour requires specific, written permission from Matthew Belmonte.

  Licence is hereby granted to copy and to distribute this software, provided
  that such copies and distributions include this copyright and licence
  notice in its entirety, and provided that the software is not altered,
  edited, amended, embellished, or otherwise modified in any way.
  Distribution of this software in other than its original form requires
  specific, written permission from Matthew Belmonte.  No fee may be charged
  for a distribution of this program or of any program based upon it, except
  for reasonable fees to cover the costs of storage media and handling.

  This software is provided 'as is', with no warranty, express or implied,
  including but not limited to warranties of merchantability or fitness for a
  particular purpose.  The user of this software assumes liability for any and
  all damages, whether direct or consequential, arising from its use.  The
  author of this software will not be liable for any such damages.

  This software makes use of the following copyrighted programs from
  _Numerical Recipes in C_, second edition, by Press, Teukolsky, Vetterling,
  and Flannery:
	realft.c
	four1.c
  These programs are not included in this distribution, and are available
  from Cambridge University Press.  The above licence for use, copying, and
  distribution DOES NOT apply to these programs.

  This software makes use of the following proprietary header file from the
  Scan386 EEG data processing system:
	sethead.h
  This file is not included in this distribution, and is available from
  NeuroScan, Inc.  The above licence for use, copying, and distribution
  DOES NOT apply to this file.

  This software makes use of the following header files from the 80x86
  implementation of the San Diego EPL system by John Hansen (San Diego
  Evoked Potential Laboratory, Neurosciences Department, University of
  California San Diego, La Jolla, CA 92093-0608, USA):
	header.h
	log.h
  These files are not included in this distribution.  The above licence for
  use, copying, and distribution DOES NOT apply to these files.  If you
  cannot obtain these files, then delete references to them and substitute
  the following dummy definitions:
#define RECSIZP 42
#define DGMAGIC 42
#define PAUZMRK 42
#define DELTMRK 42
struct header {
  short evtno, nchans, ctickt;
  char chndes[8*(CH_MAX+CH_EYE)], subdes[40], sbcdes[40], condes[1], expdes[21];
  };
*/

/*If you've found this software useful, please tell me so.  Hopefully you'll
  be able to contact me at one of the following addresses:
  mkb4@cornell.edu
  mbelmonte@ucsd.edu
  matthew@salk.edu
  Neurosciences Department, UCSD, La Jolla, CA 92093-0608, USA.
  Neuropsychology 5056, Children's Hospital, San Diego, CA 92123-4282, USA.
  Check http://salk.edu/~matthew/ for more information.
  Also, I would appreciate acknowledgement in, and a copy of, any publication
  based on data processed by this program or any derivative of it.*/

/*Define MSDOS if you're compiling for MS-DOS.  Define M_I86 if you're
  running on an Intel processor (or any processor with backwards
  byte-order).  (If you're using Microsoft QuickC, these preprocessor
  constants will be automatically defined for you.)  This code should port
  to any system but has been verified only on an 80486 MS-DOS system.
  Make sure your C compiler packs structures - the program includes a
  dynamic check for this.

  Input is in San Diego EPL (.RAW+.LOG), Scan386 continuous multiplexed
  (.CNT - but not SynAmps format). or Scan386 epoched (.EEG - not
  recommended) format.  Output in Scan386 epoched (.EEG) format.  Only the
  filename prefix is specified on the command line; given a prefix p, the
  program will search for p.RAW and p.LOG, then p.CNT, and finally p.EEG.
  The output file will be named p.EEG, except in the case where p.EEG is the
  name of the input file.  In this case the output file will be named
  OUTPUT.EEG.  To see a summary of options, execute the program with no
  parameters.  There are many parameters that can be set at compile-time -
  see scanepl.h for details.

  There is an option for Wiener filtering if you're worried about EEG
  contaminating the EOG - in practice I've not found this to be useful.
  Again, see scanepl.h for details.

  If you want to preserve your original data, make a copy of it before
  running this program.  The program will overwrite your original eye
  channels after computing the principal component vectors and receiving
  confirmation from the user.  The program will not overwrite anything
  unless you tell it explictly to do so.

  This program does not include windowing of the time-domain signal and thus
  is subject to spectral leakage during Fourier transformation.  If you want
  a window, it's not too difficult to put one in.  I have not bothered,
  because I plan to use time-domain corrections for my practical work.
  
  If your compiler can't handle the fact that the union inside the definition
  of TEEG in sethead.h doesn't have a field name, give it the name
  "teeg_union", and then in _this_ file define the preprocessor constant
  SETHEAD_HACK*/

/*Makefile for Microsoft QuickC:
OPTIONS=/c /G2 /FPi87 /Ox /W0 /AC /Zp

scanepl.exe: scanepl.obj bigdata.obj getopt.obj realft.obj four1.obj
	link /BATCH /STACK:0x2000 /EXEPACK /SEGMENTS:38 scanepl.obj bigdata.obj getopt.obj realft.obj four1.obj
scanepl.obj: scanepl.c scanepl.h sethead.h header.h
	qcl $(OPTIONS) scanepl.c
bigdata.obj: bigdata.c scanepl.h
	qcl $(OPTIONS) bigdata.c
getopt.obj: getopt.c
	qcl $(OPTIONS) getopt.c
realft.obj: realft.c
	qcl $(OPTIONS) realft.c
four1.obj: four1.c
	qcl $(OPTIONS) four1.c
*/

#include "scanepl.h"
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#ifndef MSDOS
# include <unistd.h>
# include <values.h>
# include <sys/param.h>
#else
# define MAXPATHLEN 1024
#endif
extern void realft(float[], unsigned long, int); /*_Numerical Recipes in C_*/
extern int getopt(int, char**, char*);
extern char *optarg;
extern int optind, opterr;

#define int short               /*int on PC is 2 bytes; on Sun it's 4*/
#include "header.h"             /*EPL header files*/
#include "log.h"
#undef int

#define N_ELECT 0
#define int long                /*NeuroScan header file denotes longs by "int"*/
#include "sethead.h"
#undef int
#undef N_ELECT
typedef struct {                /*this was left out of sethead.h*/
  char accept;  /*accept byte*/
  int ttype;    /*trial type*/
  int correct;  /*accuracy*/
  float rt;     /*reaction time*/
  int response; /*response type*/
  int reserved; /*not used*/
  } SWEEP_HEAD;
#ifdef SETHEAD_HACK
#define Offset teeg_union.Offset
/*jibes with a minor portability change to sethead.h*/
#endif

#ifdef M_I86
#  define msecs_to_ticks(t) (100L*(t)/tens_of_usecs_per_tick)
#  define ticks_to_msecs(t) ((((long)t)*tens_of_usecs_per_tick)/100.0)
   /*100*t takes more than one word*/
#else
#  define msecs_to_ticks(t) (100*(t)/tens_of_usecs_per_tick)
#  define ticks_to_msecs(t) ((t*tens_of_usecs_per_tick)/100.0)
#endif
typedef int Boolean;
#define true 1
#define false 0

#define RAW files[0]            /*EPL raw file*/
#define CNT files[0]            /*Scan386 continuous (raw data) file*/
#define LOG files[1]            /*event log*/
#define EEG files[2]            /*Scan386 .eeg file*/
#define NUM_FILES 3

/*Scan386 additions:*/
enum {EPL, NeuroScan, NeuroScanEpoched} input_mode;
TEEG_TYPE event_table_type;
long event_table_offset;
float reaction_time; /*set by response(), used by print_epoch()*/
long SynAmp_channel_offset,     /*block size in bytes*/
     samples_per_SynAmp_block;  /*number of clock ticks in each block*/
int points_per_waveform;        /*points per epoch in .EEG input file*/

FILE *(files[NUM_FILES]);
int
	tens_of_usecs_per_tick,
	num_samples,            /*number of values in an epoch (a power of 2)*/
	num_samples_in_output,  /*<= num_samples*/
	presample_length,       /*number of values in a presample*/
	saturation_window,      /*window length for amp. saturation test*/
	num_channels,           /*number of non-eye channels*/
	num_channels_used,      /*number of used non-eye channels*/
	num_events,             /*total number of events in the EPL log file*/
	num_cal_pulses,         /*number of calibration pulses*/
	num_eye_movements,      /*# of artefact epochs in Woestenburg filter*/
	num_quiet_epochs,       /*number of quiet epochs in Wiener filter*/
	num_rejects,            /*count of rejected epochs*/
	num_accepts,            /*count of accepted epochs*/
	event_number,           /*sequence number of current event*/
	event_code,             /*current event code*/
	event_condition;        /*condition in effect during current event*/
long
	event_time;             /*time of occurrence of current event*/
Boolean
	time_domain,            /*true for time output, false for frequency*/
	correct_artefacts,      /*true for correction, false for rejection*/
	ask_for_channels,       /*true iff no channels specified on cmd. line*/
	use_channel[CH_MAX];    /*tells which channels we're interested in*/
float
	rejection_threshold,
	calibration[CH_EYE+CH_MAX],
	avg_quiet_eog_power[(EP_LEN/EP_INC)/2-1],
	summed_Wiener_coefficients[(EP_LEN/EP_INC)/2-1],
	summed_squared_Wiener_coefficients[(EP_LEN/EP_INC)/2-1];
extern float blink_filter[CH_MAX][EP_LEN/EP_INC-2];

struct header hdr;
SETUP erp;
ELECTLOC channel_header[CH_EYE+CH_MAX];

/*the smallest power of 2 that's greater than or equal to the given N*/
int FFT_length(N)
int N;
  {
  int length;
  length = 1;
  while(length < N)
    length <<= 1;
  return(length);
  }

/*checks of constants defined at compile time*/
void check_static_definitions()
  {
  Boolean result;
  result = false;
  if((sizeof(SETUP)!=900)||(sizeof(ELECTLOC)!=75)||(sizeof(TEEG)!=9)||(sizeof(SWEEP_HEAD)!=13))
    {
    printf("Your compiler doesn't pack structures - if using Microsoft QuickC, give the /Zp option.\n");
    result = true;
    }
  if(input_mode != NeuroScanEpoched)
    {
    num_samples_in_output = EP_OUTLEN/EP_INC;
    num_samples = FFT_length(num_samples_in_output);
    if(num_samples > EP_LEN/EP_INC)
      {
      fprintf(stderr, "EP_LEN >= %d\n", num_samples*EP_INC);
      result = true;
      }
    }
  if(result)
    {
    fprintf(stderr, "Change the above static definition(s) and recompile.\n");
    exit(1);
    }
  }

/*close all files*/
void cleanup()
  {
  register int file;
  for(file = NUM_FILES-1; file >= 0; file--)
    fclose(files[file]);
  }

#ifdef M_I86
#  define swap_bytes(word) word
#  define reverse_long(double_word) double_word
#  define reverse_float(double_word) double_word
#  define reverse_double(quad_word) quad_word
#else

short swap_bytes(word)
short word;
  {
  return(((word >> 8) & 0xff)|(word << 8));
  }

void reverse_bytes(src, dest, len)
char *src, *dest;
int len;
  {
  register int i;
  for(i = 0; i != len; i++)
    dest[i] = src[len-1-i];
  }

long reverse_long(double_word)
long double_word;
  {
  long result;
  reverse_bytes((char *)&double_word, (char *)&result, sizeof(long));
  return(result);
  }

float reverse_float(double_word)
float double_word;
  {
  float result;
  reverse_bytes((char *)&double_word, (char *)&result, sizeof(long));
  return(result);
  }

double reverse_double(quad_word)
double double_word;
  {
  double result;
  reverse_bytes((char *)&quad_word, (char *)&result, sizeof(long));
  return(result);
  }
#endif

/*Massage data from NeuroScan format into a form that the rest of this program
  likes to deal with.  Set num_channels, tens_of_usecs_per_tick, calibration[],
  event_table_type, event_table_offset, SynAmp_channel_offset,
  samples_per_SynAmp_block, points_per_waveform.*/
void set_NeuroScan_mode(CNT_filename)
char *CNT_filename;
  {
  register int channel;
  TEEG event_table_header;
  if(fread(&erp, sizeof(erp), 1, CNT) != 1) /*experiment & subject labels*/
    exit(1);
  sprintf(hdr.subdes, "%s %c %d %cH", erp.patient, erp.sex, swap_bytes(erp.age), erp.hand);
  sprintf(hdr.sbcdes, "%s %s", erp.date, erp.time);
  hdr.condes[0] = '\0';
  strncpy(hdr.expdes, erp.label, 20);
  num_channels = swap_bytes(erp.nchannels)-CH_EYE;
  tens_of_usecs_per_tick = 100000/swap_bytes(erp.rate);
  SynAmp_channel_offset = reverse_long(erp.ChannelOffset);
  samples_per_SynAmp_block = SynAmp_channel_offset/((CH_EYE+num_channels)*sizeof(short));
  points_per_waveform = swap_bytes(erp.pnts);
  for(channel = 0; channel != CH_EYE+num_channels; channel++) /*chan. labels*/
    {
    if(fread(channel_header+channel, sizeof(*channel_header), 1, CNT) != 1)
      exit(1);
    calibration[channel] = reverse_float(channel_header[channel].sensitivity)*reverse_float(channel_header[channel].calib)/VOLTAGE_MULTIPLIER;
    strncpy(hdr.chndes+8*channel, channel_header[channel].lab, 8);
    }
  if(input_mode==NeuroScan)
    {
    LOG = fopen(CNT_filename,   /*tail of CNT file masquerades as LOG*/
#ifdef MSDOS
	"rb"
#else
	"r"
#endif
	);
    fseek(LOG, reverse_long(erp.EventTablePos), SEEK_SET);
    if(fread(&event_table_header, sizeof(event_table_header), 1, LOG) != 1)
      exit(1);
    event_table_type = event_table_header.Teeg;
    num_events = reverse_long(event_table_header.Size)/(event_table_type == TEEG_EVENT_TAB1 ? sizeof(EVENT1) : sizeof(EVENT2));
    fseek(LOG, reverse_long(event_table_header.Offset), SEEK_CUR);
    event_table_offset = ftell(LOG);
    }
  else
    {
    /*input_mode==NeuroScanEpoched*/
    num_events = swap_bytes(erp.compsweeps); /*in an epoched file each sweep is one event*/
    num_samples_in_output = (swap_bytes(erp.pnts)-msecs_to_ticks(-reverse_float(erp.xmin)*1000.0))/EP_INC;
    if(num_samples_in_output > EP_LEN/EP_INC)
      num_samples_in_output = EP_LEN/EP_INC;
    num_samples = FFT_length(num_samples_in_output);
    if(num_samples > EP_LEN/EP_INC)
      {
      fprintf(stderr, "Change EP_LEN to %d\n", num_samples*EP_INC);
      fclose(CNT);
      exit(1);
      }
    }
  }

/*Open inputs files and set num_events*/
void open_files(prefix)
char *prefix;
  {
#ifdef MSDOS
  static char extensions[][5] = {".raw", ".log", ".eeg"}, modes[][3] = {"rb", "rb", "wb"};
#else
  static char extensions[][10] = {".raw", ".log", ".eeg"}, modes[][2] = {"r", "r", "w"};
#endif
  register int file;
  char filename[MAXPATHLEN];
  struct stat stbuf;
  input_mode = EPL; /*until we find otherwise*/
  for(file = 0; file != NUM_FILES; file++)
    {
    if((files[file] = fopen(strcat(strncpy(filename, input_mode==NeuroScanEpoched? "output": prefix, MAXPATHLEN-9), extensions[file]), modes[file])) == NULL)
      {
      /*If we didn't find an EPL raw file, look for a NeuroScan CNT file:*/
      if((file == 0) && ((CNT = fopen(strcat(strncpy(filename, prefix, MAXPATHLEN-9), ".cnt"),
#ifdef MSDOS
	"rb"
#else
	"r"
#endif
	)) != NULL))
	{
	input_mode = NeuroScan;
	set_NeuroScan_mode(filename);
	file++;
	}
      else if((file == 0) && ((CNT = fopen(strcat(strncpy(filename, prefix, MAXPATHLEN-9), ".eeg"),
#ifdef MSDOS
	"rb"
#else
	"r"
#endif
	)) != NULL))
	{
	input_mode = NeuroScanEpoched;
	set_NeuroScan_mode("");
	file++;
	}
      else
	{
	perror(filename);
	while(--file >= 0)
	  fclose(files[file]);
	exit(1);
	}
      }
    }
  if(input_mode == EPL)
    {
    if(fstat(fileno(files[1]), &stbuf) == -1)
      {
      perror(prefix);
      cleanup();
      exit(1);
      }
    num_events = stbuf.st_size/LOGSIZE;
    }
  }

/*update header information, and then close files*/
void close_files()
  {
  register int file, channel;
  float epoch_length_in_seconds;
  rewind(EEG);
  epoch_length_in_seconds = ticks_to_msecs(EP_INC*num_samples_in_output)/1000.0;
  if(input_mode==EPL)
    {
    register int i;
    /*Scan386 setup fields aren't yet filled; get the subject's name from the
      EPL header and fill in the rest of the setup fields with default values:*/
    for(i = 0; i != sizeof(erp); i++)
      ((char *)&erp)[i] = 0;
    strcpy(erp.rev, "Version 3.0");
    hdr.subdes[19] = '\0';
    strcpy(erp.patient, hdr.subdes);
    erp.sex = 'U';
    erp.hand = 'U';
    }
  erp.type = 0; /*epoched EEG*/
  erp.mean_age = reverse_float(0.0);
  erp.stdev = reverse_float(0.0);
  erp.n = swap_bytes(num_accepts);
  erp.compfile[0] = '\0';
  erp.sortfile[0] = '\0';
  erp.NumEvents = reverse_long(0L);
  erp.compoper = 0;
  erp.avgmode = 0;
  erp.review = 1;
  erp.nsweeps = swap_bytes(num_accepts);
  erp.compsweeps = swap_bytes(num_accepts);
  erp.acceptcnt = swap_bytes(num_accepts);
  erp.rejectcnt = swap_bytes(0);
  erp.pnts = swap_bytes(num_samples_in_output); /*already subtracted 2 if frequency domain*/
  erp.nchannels = swap_bytes(time_domain? CH_EYE+num_channels_used: num_channels_used);
  erp.avgupdate = swap_bytes(1);
  erp.domain = 0;
  erp.variance = 0;
  erp.rate = swap_bytes((input_mode==EPL? (short)(100000L/tens_of_usecs_per_tick): swap_bytes(erp.rate))/EP_INC);
  erp.scale = reverse_double(1.0);
  erp.veogcorrect = 0;
  erp.heogcorrect = 0;
  erp.aux1correct = 0;
  erp.aux2correct = 0;
  erp.veogtrig = reverse_float(10.0);
  erp.heogtrig = reverse_float(10.0);
  erp.aux1trig = reverse_float(10.0);
  erp.aux2trig = reverse_float(10.0);
  erp.heogchnl = swap_bytes(0);
  erp.veogchnl = swap_bytes(0);
  erp.aux1chnl = swap_bytes(0);
  erp.aux2chnl = swap_bytes(0);
  erp.veogdir = 0;
  erp.heogdir = 0;
  erp.aux1dir = 0;
  erp.aux2dir = 0;
  erp.veog_n = swap_bytes(0);
  erp.heog_n = swap_bytes(0);
  erp.aux1_n = swap_bytes(0);
  erp.aux2_n = swap_bytes(0);
  erp.veogmaxcnt = swap_bytes(20);
  erp.heogmaxcnt = swap_bytes(20);
  erp.aux1maxcnt = swap_bytes(20);
  erp.aux2maxcnt = swap_bytes(20);
  erp.veogmethod = 0;
  erp.heogmethod = 0;
  erp.aux1method = 0;
  erp.aux2method = 0;
  erp.AmpSensitivity = reverse_float(0.0);
  erp.LowPass = 0;
  erp.HighPass = 0;
  erp.Notch = 0;
  erp.AutoClipAdd = 0;
  erp.baseline = 0;
  erp.offstart = reverse_float(0.0);
  erp.offstop = reverse_float(epoch_length_in_seconds);
  erp.reject = 0;
  erp.rejstart = reverse_float(0.0);
  erp.rejstop = reverse_float(epoch_length_in_seconds);
  erp.rejmin = reverse_float(-rejection_threshold/2.0);
  erp.rejmax = reverse_float(rejection_threshold/2.0);
  erp.trigtype = 2;
  erp.trigval = reverse_float(0.0);
  erp.trigchnl = 16;
  erp.trigmask = swap_bytes(0);
  erp.trigisi = reverse_float(1.0);
  erp.trigmin = reverse_float(0.0);
  erp.trigmax = reverse_float(5.0);
  erp.trigdir = 0;
  erp.Autoscale = 1;
  erp.n2 = swap_bytes(0);
  erp.dir = 1; /*positive up in plot*/
  erp.dispmin = reverse_float(-rejection_threshold/2.0);
  erp.dispmax = reverse_float(rejection_threshold/2.0);
  erp.xmin = reverse_float(0.0);
  erp.xmax = reverse_float(epoch_length_in_seconds);
  erp.AutoMin = reverse_float(0.0);
  erp.AutoMax = reverse_float(epoch_length_in_seconds);
  erp.zmin = reverse_float(0.0);
  erp.zmax = reverse_float(0.1);
  erp.lowcut = reverse_float(0.01);
  erp.highcut = reverse_float(100.0);
  erp.common = 0;
  erp.savemode = 0;
  erp.manmode = 0;
  strcpy(erp.ref, "REF");
  erp.Rectify = 0;
  erp.DisplayXmin = reverse_float(0.0);
  erp.DisplayXmax = reverse_float(epoch_length_in_seconds);
  erp.phase = 0;
  strcpy(erp.screen, "BLANK");
  erp.CalMode = swap_bytes(0);
  erp.CalMethod = swap_bytes(0);
  erp.CalUpdate = swap_bytes(5);
  erp.CalBaseline = swap_bytes(1);
  erp.CalSweeps = swap_bytes(20);
  erp.CalAttenuator = reverse_float(50000.0);
  erp.CalPulseVolt = reverse_float(1.0);
  erp.CalPulseStart = reverse_float(0.0);
  erp.CalPulseStop = reverse_float(1.0);
  erp.CalFreq = reverse_float(0.8333333);
  strcpy(erp.taskfile, "-------");
  strcpy(erp.seqfile, "--------");
  erp.SpectMethod = 0;
  erp.SpectScaling = 0;
  erp.SpectWindow = 0;
  erp.SpectWinLength = reverse_float(0.1);
  erp.SpectOrder = 0;
  erp.NotchFilter = 0;
  erp.FspStopMethod = swap_bytes(0);
  erp.FspStopMode = swap_bytes(0);
  erp.FspFValue = reverse_float(2.5);
  erp.FspPoint = swap_bytes(0);
  erp.FspBlockSize = swap_bytes(0);
  erp.FspP1 = swap_bytes(0);
  erp.FspP2 = swap_bytes(0);
  erp.FspAlpha = reverse_float(0.0);
  erp.FspNoise = reverse_float(0.0);
  erp.FspV1 = swap_bytes(0);
  erp.fratio = reverse_float(0.8338761925697327);
  erp.minor_rev = 9;
  erp.eegupdate = swap_bytes(1);
  erp.compressed = 0;
  erp.xscale = reverse_float(0.0);
  erp.yscale = reverse_float(0.0);
  erp.xsize = reverse_float(40.0);
  erp.ysize = reverse_float(20.0);
  erp.disptype = 0;
  erp.CommonChnl = 0;
  erp.Xtics = 0;
  erp.Xrange = 0;
  erp.Ytics = 0;
  erp.Yrange = 0;
  erp.XScaleValue = reverse_float(2500.0);
  erp.XScaleInterval = reverse_float(1280.0);
  erp.YScaleValue = reverse_float(100.0);
  erp.YScaleInterval = reverse_float(50.0);
  erp.ScaleToolX1 = reverse_float(474.2577514648437);
  erp.ScaleToolY1 = reverse_float(41.55833435058594);
  erp.ScaleToolX2 = reverse_float(536.7577514648437);
  erp.ScaleToolY2 = reverse_float(21.55833435058594);
  erp.port = swap_bytes(0x2CB);
  erp.NumSamples = swap_bytes(0);
  erp.FilterFlag = 0;
  erp.LowCutoff = reverse_float(5.12);
  erp.LowPoles = swap_bytes(2);
  erp.HighCutoff = reverse_float(64.0);
  erp.HighPoles = swap_bytes(2);
  erp.FilterType = 3; /*low-pass*/
  erp.FilterDomain = 1; /*time domain*/
  erp.SnrFlag = 0;
  erp.CoherenceFlag = 0;
  erp.ContinousType = 1;
  erp.EventTablePos = reverse_long(0L); /*no event table in epoched file*/
  erp.ContinousSeconds = reverse_float(0.0);
  erp.ChannelOffset = swap_bytes(0); /*used only w/ SynAmps*/
  erp.AutoCorrectFlag = 0;
  erp.DCThreshold = 0;
  fwrite(&erp, sizeof(erp), 1, EEG);
  for(channel = time_domain? 0: CH_EYE; channel != num_channels+CH_EYE; channel++)
    if((channel < CH_EYE) || use_channel[channel-CH_EYE])
      {
      if(input_mode == EPL)
	{
	strncpy(channel_header[channel].lab, hdr.chndes+8*(channel), 8);
	channel_header[channel].lab[8] = '\0';
	channel_header[channel].reference = -1;
	channel_header[channel].skip = 0;
	channel_header[channel].reject = 0;
	channel_header[channel].display = 1;
	channel_header[channel].bad = 0;
	channel_header[channel].avg_reference = 0;
	channel_header[channel].ClipAdd = 0;
	channel_header[channel].x_coord = reverse_float(9.0+(channel % 4)*128.0);
	channel_header[channel].y_coord = reverse_float(35.0+(channel / 4)*42.0);
	channel_header[channel].veog_wt = reverse_float(0.0);
	channel_header[channel].veog_std = reverse_float(0.0);
	channel_header[channel].snr = reverse_float(0.0);
	channel_header[channel].heog_wt = reverse_float(0.0);
	channel_header[channel].heog_std = reverse_float(0.0);
	channel_header[channel].Filtered = 0;
	channel_header[channel].Fsp = 0;
	channel_header[channel].aux1_wt = reverse_float(0.0);
	channel_header[channel].aux1_std = reverse_float(0.0);
	channel_header[channel].Gain = 0;
	channel_header[channel].HiPass = 0;
	channel_header[channel].LoPass = 0;
	channel_header[channel].Page = 0;
	channel_header[channel].Size = 0;
	channel_header[channel].Impedance = 0;
	channel_header[channel].PhysicalChnl = channel;
	channel_header[channel].Rectify = 0;
	}
      channel_header[channel].n = swap_bytes(1);
      channel_header[channel].baseline = swap_bytes(0); /*output is baseline-corrected*/
      channel_header[channel].sensitivity = reverse_float(1.0);
      channel_header[channel].calib = reverse_float(time_domain? 1.0: 1.0/AMPLITUDE_SCALE);
      fwrite(channel_header+channel, sizeof(*channel_header), 1, EEG);
      }
  cleanup();
  }

/*generate the .lis file*/
void generate_report(prefix)
char *prefix;
  {
  register int channel, sample;
  FILE *lis;
  char lis_name[MAXPATHLEN];
  if((lis = fopen(strcat(strncpy(lis_name, prefix, MAXPATHLEN-5), ".lis"), "w")) == NULL)
    perror(lis_name);
  else
    {
    fprintf(lis, "%.40s\n%.40s\n%.40s\n%.40s\n", hdr.subdes, hdr.sbcdes, hdr.condes, hdr.expdes);
    if(input_mode == EPL)
      {
      fprintf(lis, "%d calibration pulses.\n", num_cal_pulses);
      if(num_cal_pulses > CA_PULSES)
	fprintf(lis, "(Used only %d of them - recompile w/ larger CA_PULSES.)\n", CA_PULSES);
      }
    if(correct_artefacts)
      fprintf(lis, "%d blink epochs in Woestenburg regression computation;\n%d quiet epochs averaged into Wiener filter.\n", num_eye_movements, num_quiet_epochs);
    fprintf(lis, "%d epochs, %d accepted, %d (%d%%) rejected at threshold of %fuV.\n", num_rejects+num_accepts, num_accepts, num_rejects, (int)(100L*num_rejects/(num_rejects+num_accepts)), rejection_threshold);
    fprintf(lis, "%d Hertz ", 100000/tens_of_usecs_per_tick);
    fprintf(lis, "(%.2f ms/sample) ", tens_of_usecs_per_tick/100.0);
    fprintf(lis, "on %d channels:\n", CH_EYE+num_channels);
    for(channel = 0; channel != CH_EYE; channel++)
      fprintf(lis, " %2d %-8.8s %3.0f%%\n", channel+1, hdr.chndes+8*channel, (input_mode == EPL ? 100.0/CA_SCALE : 20.0)/calibration[channel]);
    for(channel = 0; channel != num_channels; channel++)
      fprintf(lis, "%c%2d %-8.8s %3.0f%%\n", use_channel[channel]? '+' : ' ', channel+CH_EYE+1, hdr.chndes+8*(channel+CH_EYE), (input_mode == EPL ? 100.0/CA_SCALE : 20.0)/calibration[channel+CH_EYE]);
    if(correct_artefacts)
      for(channel = 0; channel != num_channels; channel++)
	if(use_channel[channel])
	  {
	  fprintf(lis, "CHANNEL %d (%.8s) FILTER VALUES:\n", channel+1+CH_EYE, hdr.chndes+8*(CH_EYE+channel));
	  for(sample = 0; sample != num_samples/2-1; sample++)
	    fprintf(lis, "%2.3f Hz : %f%+fi * %f SD=%f\n", 1e+5*(sample+1)/(
#ifdef M_I86
(long)
#endif
		num_samples*tens_of_usecs_per_tick*EP_INC), blink_filter[channel][2*sample], blink_filter[channel][2*sample+1], summed_Wiener_coefficients[sample]/num_accepts, sqrt((summed_squared_Wiener_coefficients[sample]-summed_Wiener_coefficients[sample]*summed_Wiener_coefficients[sample]/num_accepts)/num_accepts));
	  putc('\n', lis);
	  }
    fclose(lis);
    }
  }

void fatal_error(message)
char *message;
  {
  if(message == (char *)0)
    perror((char *)0);
  else
    fprintf(stderr, "%s\n", message);
  cleanup();
  exit(1);
  }

/*read the header from the EPL file into hdr, set num_channels and
  tens_of_usecs_per_tick, and zero the Scan386 headers*/
void read_header()
  {
  register int i, channel;
  if(fread((char *)&hdr, sizeof(hdr), 1, RAW) != 1)
    fatal_error((char *)0);
  if(swap_bytes(hdr.evtno) != DGMAGIC)
    fatal_error("bad magic number");
  num_channels = swap_bytes(hdr.nchans)-CH_EYE;
  tens_of_usecs_per_tick = swap_bytes(hdr.ctickt);
  for(i = 0; i != sizeof(erp); i++)
    ((char *)&erp)[i] = 0;
  for(channel = 0; channel != CH_EYE+num_channels; channel++)
    for(i = 0; i != sizeof(*channel_header); i++)
      ((char *)(channel_header+channel))[i] = 0;
  }

/*Checks that depend upon information from the EPL raw file header.  Set
  presample_length and saturation_window.*/
void check_dynamic_definitions()
  {
  if(num_channels <= CH_EYE)
    fatal_error("no head channels");
  if(num_channels > CH_EYE+CH_MAX)
    fatal_error("too many channels");
  presample_length = msecs_to_ticks(RE_PRESAMPLE);
  if(input_mode == NeuroScanEpoched)
    {
    int max_presample_length;
    max_presample_length = msecs_to_ticks((int)(-reverse_float(erp.xmin)*1000.0));
    if(presample_length > max_presample_length)
      presample_length = max_presample_length;
    }
/*don't divide by EP_INC above, since it's pointless to decimate the presample*/
  if(presample_length > RE_PSLEN)
    fatal_error("presample too long");
  saturation_window = msecs_to_ticks(RE_WINSIZE);
  }

/*Set time_domain, correct_artefacts, use_channel[], and num_channels_used.*/
void process_options(argc, argv)
int argc;
char *argv[];
  {
  register int arg, channel;
  time_domain = true;
  correct_artefacts = false;
  rejection_threshold = RE_THRESHOLD;
  while((arg = getopt(argc, argv, "tfcr:")) != -1)
    switch(arg)
      {
      case 't': time_domain = true; break;
      case 'f': time_domain = false; break;
      case 'c': correct_artefacts = true; break;
      case 'r': correct_artefacts = false;
		sscanf(optarg, "%f", &rejection_threshold);
      }
  if(optind == argc)
    {
    fprintf(stderr, "usage: %s [-tfrc] [channel numbers] file_prefix\n", argv[0]);
    fprintf(stderr, "options:\n");
    fprintf(stderr, "\t-t: time domain [default]\n");
    fprintf(stderr, "\t-f: frequency domain\n");
    fprintf(stderr, "\t-r <threshold>: reject artefacts [default at %fuV]\n", RE_THRESHOLD);
    fprintf(stderr, "\t-c: correct artefacts\n");
    exit(2);
    }
  if(!time_domain)
    {
    if(num_samples != num_samples_in_output)
      {
      fprintf(stderr, "for frequency-domain output, EP_LEN/EP_INC must be a power of 2\n");
      exit(1);
      }
    num_samples_in_output = num_samples-2; /*omit DC & Nyquist components*/
    }
  num_channels_used = 0;
  ask_for_channels = optind == argc-1;
  if(!ask_for_channels)
    {
    for(channel = 0; channel != CH_MAX; channel++)
      use_channel[channel] = false;
    for(arg = optind; arg != argc-1; arg++)
      {
      channel = atoi(argv[arg])-CH_EYE-1;
      if((channel < 0) || (channel >= CH_MAX))
	fatal_error("bogus channel");
      if(!use_channel[channel])
	{
	use_channel[channel] = true;
	num_channels_used++;
	}
      }
    }
  }

void ask_channels()
  {
  register int channel, arg;
  if(ask_for_channels)
    for(channel = 0; channel != num_channels; channel++)
      {
      printf("Use channel %d (%-8.8s)? ", channel+CH_EYE+1, hdr.chndes+8*(channel+CH_EYE));
      do arg=getchar(); while((arg != 'y') && (arg != 'Y') && (arg != 'n') && (arg != 'N'));
      while(getchar() != '\n')
	;
      use_channel[channel] = (arg == 'y') || (arg == 'Y');
      if(use_channel[channel]) num_channels_used++;
      }
  else
    for(channel = num_channels; channel != CH_MAX; channel++)
      if(use_channel[channel])
	fatal_error("bogus channel");
  }

short read_swapped_word(stream)
FILE *stream;
  {
  short result;
  result = getc(stream);
  result |= getc(stream) << 8;
  return(result);
  }

/*the EEG exactly as it appears in the EPL or Scan386 raw data file*/
int raw_eeg(time, channel)
long time;
int channel;
  {
  fseek(RAW,
    input_mode==EPL?
	(long)(sizeof(struct header)
	+ (time - time % RECSIZP) * (CH_EYE+1+num_channels) * sizeof(short)
	+ RECSIZP * sizeof(short)
	+ (time % RECSIZP) * (CH_EYE+num_channels) * sizeof(short)
	+ (CH_EYE+channel) * sizeof(short)):
      input_mode==NeuroScan? 
	(long)(sizeof(SETUP)
	+ (CH_EYE+num_channels)*sizeof(ELECTLOC)
	+ (samples_per_SynAmp_block==0L?
	/*Scan386 without SynAmps*/
	  (time*(CH_EYE+num_channels)+CH_EYE+channel)*sizeof(short):
	/*Scan386 with SynAmps*/
	  time/samples_per_SynAmp_block*SynAmp_channel_offset
	  + (CH_EYE+channel)*sizeof(short)*samples_per_SynAmp_block
	  + time%samples_per_SynAmp_block*sizeof(short))):
      /*input_mode==NeuroScanEpoched*/
	(long)(sizeof(SETUP)
	+ (CH_EYE+num_channels)*sizeof(ELECTLOC)
	+ (time/points_per_waveform)*(sizeof(SWEEP_HEAD)+(CH_EYE+num_channels)*sizeof(short)*points_per_waveform)
	+ sizeof(SWEEP_HEAD) + ((time%points_per_waveform)*(CH_EYE+num_channels)+CH_EYE+channel)*sizeof(short)),
    SEEK_SET);
  return((int)(read_swapped_word(RAW)));
  }

/*Retrieve the next event, code, time, and condition from the EPL log file, and
  increment event_number.*/
void next_event()
  {
  if(input_mode == EPL)
    {
    event_code = (int)(read_swapped_word(LOG));
    event_time = ((long)(read_swapped_word(LOG))) << 16;
    event_time |= ((long)(read_swapped_word(LOG))) & 0xffff;
    event_condition = getc(LOG);
    getc(LOG);
    }
  else if(input_mode == NeuroScan)
    {
    event_code = (int)(read_swapped_word(LOG));
    getc(LOG); getc(LOG);
    event_time = ((long)(read_swapped_word(LOG))) & 0xffff;
    event_time |= ((long)(read_swapped_word(LOG))) << 16;
    event_time = (event_time - sizeof(SETUP) - (CH_EYE+num_channels)*sizeof(ELECTLOC))/((CH_EYE+num_channels)*sizeof(short));
    /*calculation of event_time will be different for SynAmps, but I'm not sure how*/
    event_condition = 0; /*a dummy value; don't know what to put here*/
    }
  else /*input_mode == NeuroScanEpoched*/
    {
    SWEEP_HEAD sweep_header;
    fseek(CNT, sizeof(SETUP)
	+ (CH_EYE+num_channels)*sizeof(ELECTLOC)
	+ (long)event_number*(sizeof(SWEEP_HEAD)+(CH_EYE+num_channels)*sizeof(short)*points_per_waveform),
      SEEK_SET);
    fread(&sweep_header, sizeof(sweep_header), 1, CNT);
    event_code = swap_bytes(sweep_header.ttype);
    event_time = event_number*(long)(swap_bytes(erp.pnts)) + msecs_to_ticks((int)(-reverse_float(erp.xmin)*1000.0));
    event_condition = 0; /*a dummy value; don't know what to put here*/
    reaction_time = reverse_float(sweep_header.rt);
    }
  event_number++;
  }

/*Seek to the specified event, and set all variables accordingly.*/
void goto_event(evt)
int evt;
  {
  if(evt > 0)
    {
    event_number = evt-1;
    if((input_mode==EPL)||(input_mode==NeuroScan)) fseek(LOG,
      input_mode == EPL ?
	(long)(event_number*LOGSIZE) :
	event_table_offset + event_number*(event_table_type == TEEG_EVENT_TAB1 ? sizeof(EVENT1) : sizeof(EVENT2)),
      SEEK_SET);
    next_event();
    }
  else
    {
    event_number = 0;
    if((input_mode==EPL)||(input_mode==NeuroScan))
      fseek(LOG, input_mode == EPL? 0L: event_table_offset, SEEK_SET);
    }
  }

void swap(b, i, j)
int *b;
int i, j;
  {
  register int temp;
  temp = b[i];
  b[i] = b[j];
  b[j] = temp;
  }

/*Find the median of an array b of integers.  This is accomplished using a
  partial quicksort; thus the contents of b end up as a permutation of the
  original sequence.  A median of the calibration values is used instead of
  their mean, because it is insensitive to amplifier saturation that might
  shrink some of the calibration pulses.
  This routine isn't an important part of the program, but it is pretty cool;
  it runs in linear expected time (cN + cN/2 +cN/4 + cN/8 + ... = 2cN).  The key
  to this is deciding what parts of the quicksort really need to be done (see
  the invariant below).
*/
int median(b, num_pulses)
int *b, num_pulses;
  {
  register int h, k, m, n;
  m = 0;
  n = num_pulses-1;
  /*inv: m <= num_pulses/2 < n and b[m..n] is a permutation of the subsequence
    b[m..n] of sorted sequence b[0..num_pulses-1]*/
  do
    {
    h = m+1;
    k = n;
    /*inv: b[m+1..h-1] <= b[m] and b[k+1..n] > b[m]*/
    while(h <= k)
      {
      while((h <= k) && (b[h] <= b[m]))
	h++;
      while((h <= k) && (b[k] > b[m]))
	k--;
      if(h < k)
	swap(b, h++, k--);
      } /*elihw*/
    /*k=h-1*/
    if(k==m)
      m++;
    else if(h==n+1)
      swap(b, m, n--);
    else
      {
      swap(b, m, k);
      if(k <= num_pulses/2)
	m=k;
      else
	n=k;
      } /*fi*/
    /*b[k] is now in sorted position*/
    } while(k != num_pulses/2);
  return(b[k]);
  }

/*Set calibration[] s.t. the numbers produced by eeg() will be in microvolts,
  and set num_cal_pulses.*/
void calibrate()
  {
  register int channel;
  extern int pulse_values[CH_EYE+CH_MAX][CA_PULSES];
  event_number = 0;
/*no events have been read from the log file*/
  next_event();
  while((event_number != 1+CA_PULSES) && (event_code == EV_CALIBRATE))
    {
    for(channel = 0; channel != CH_EYE+num_channels; channel++)
      pulse_values[channel][event_number-1] = -raw_eeg(event_time+msecs_to_ticks(CA_LATENCY), channel-CH_EYE);
    for(channel = 0; channel != CH_EYE+num_channels; channel++)
      pulse_values[channel][event_number-1] += raw_eeg(event_time+msecs_to_ticks(CA_LATENCY+CA_LENGTH/2), channel-CH_EYE);
    next_event();
    }
/*the current event is the event following the last used calibration pulse*/
  if(event_number == 1)
  /*in the absence of calibration pulses, assume that each channel is nominal*/
    for(channel = 0; channel != CH_EYE+num_channels; channel++)
      calibration[channel] = 1.0/CA_SCALE;
  else
    {
    for(channel = 0; channel != CH_EYE+num_channels; channel++)
      calibration[channel] = CA_AMPLITUDE/median(pulse_values[channel], event_number-1);
    if(event_code == EV_CALIBRATE)
    /*there were more pulses than could be handled, so discard the extras*/
      {
      do
	next_event();
      while(event_code == EV_CALIBRATE);
      }
    }
/*the current event is the event following the last calibration pulse*/
  num_cal_pulses = event_number-1;
  }

/*the EEG scaled by the appropriate calibration value*/
float eeg(time, channel)
long time;
int channel;
  {
  return(raw_eeg(time, channel)*calibration[CH_EYE+channel]);
  }

/*the difference between the maximum and minimum peak amplitudes in the given
  interval*/
float ppa(signal, length)
float *signal;
int length;
  {
  register int t;
  float min, max;
  min = signal[0];
  max = signal[0];
  for(t = 1; t != length; t++)
    {
    if(signal[t] < min) min = signal[t];
    if(signal[t] > max) max = signal[t];
    }
  return(max-min);
  }

/*the average of the signal over the given interval*/
float avg(signal, length)
float *signal;
int length;
  {
  register int t;
  float sum;
  sum = 0.0;
  for(t = 0; t != length; t++)
    sum += signal[t];
  return(sum/length);
  }

/*true iff the signal is constant for at least saturation_window ticks
  (RE_WINSIZE milliseconds)*/
Boolean saturated(signal, length)
float *signal;
int length;
  {
  int t, flatline_length;
  if(length < saturation_window)
    return(false);
  else
    {
    flatline_length = 1;
    for(t = 1; (flatline_length != saturation_window) && (t != length); t++)
      flatline_length = signal[t]==signal[t-1]? 1+flatline_length : 1;
    return(flatline_length == saturation_window);
    }
  }

void ComplexMultiply(m0, m1, prod)
float *m0, *m1, *prod;
  {
  *(prod++) = m0[0]*m1[0] - m0[1]*m1[1];
  *prod = m0[1]*m1[0] + m0[0]*m1[1];
  }

void load_eeg(data, time)
float data[][EP_LEN/EP_INC];
long time;
  {
  register int channel, sample;
  for(sample = 0; sample != num_samples; sample++)
    for(channel = -CH_EYE; channel != num_channels; channel++)
      data[channel][sample] = eeg(time+EP_INC*sample, channel);
  }

/*Determine Woestenburg filter coefficients for each EEG channel and each
  frequency by computing complex correlation coefficients over all epochs in
  which the EOG amplitude exceeds the rejection threshold.  Also compute the
  average amplitude spectrum of non-blink EOG, for use in Wiener filtering.*/
void find_regression_coefficients()
  {
  register int sample, channel;
  int epoch, num_epochs; /*used only for printing progress report*/
  Boolean not_saturated;
  long time;
  float term[2], dynamic_range;
  float veog_data[EP_LEN/EP_INC],
    summed_eog_power[(EP_LEN/EP_INC)/2-1];
  extern float data[CH_EYE+CH_MAX][EP_LEN/EP_INC];
  num_eye_movements = 0;
  num_quiet_epochs = 0;
  for(sample = 0; sample != num_samples/2-1; sample++)
    {
    avg_quiet_eog_power[sample] = 0.0;
    summed_eog_power[sample] = 0.0;
    summed_Wiener_coefficients[sample] = 0.0;
    summed_squared_Wiener_coefficients[sample] = 0.0;
    for(channel = 0; channel != num_channels; channel++)
      {
      blink_filter[channel][2*sample] = 0.0;
      blink_filter[channel][2*sample+1] = 0.0;
      }
    }
  do
    {
    time = event_time;
    do
      {
      next_event();
      printf("\rSEARCHING EVENT %d OF %d", event_number, num_events);
      }
    while((event_code != PAUZMRK) && (event_code != DELTMRK) && (event_number != num_events));
  /*EEG has been recorded continuously from time to event_time*/
    putchar('\n');
    num_epochs = (int)((event_time-time)/(EP_INC*num_samples));
    epoch = 0;
    while(time <= event_time - EP_INC*num_samples)
      {
      printf("\rCOMPUTING EPOCH %d OF %d", ++epoch, num_epochs);
      load_eeg(data+CH_EYE, time);
      if(!saturated(data[CH_EYE+CH_LOWER_EYE], num_samples) && !saturated(data[CH_EYE+CH_UPPER_EYE], num_samples))
      /*neither of the VEOG channels is saturated*/
	{
	for(sample = 0; sample != num_samples; sample++)
	  veog_data[sample] = data[CH_EYE+CH_UPPER_EYE][sample] - data[CH_EYE+CH_LOWER_EYE][sample];
	dynamic_range = ppa(veog_data, num_samples);
	if(dynamic_range < RE_BRAIN_THRESHOLD)
	/*this is not an eye movement*/
	  {
	  realft(veog_data-1, num_samples, 1);
	  for(sample = 0; sample != num_samples/2-1; sample++)
	    avg_quiet_eog_power[sample] += veog_data[2*sample+2]*veog_data[2*sample+2]+veog_data[2*sample+3]*veog_data[2*sample+3];
	  num_quiet_epochs++;
	  }
	else if(dynamic_range > RE_EYE_THRESHOLD)
	/*this is an eye movement*/
	  {
	  not_saturated = true;
	  for(channel = 0; channel != num_channels; channel++)
	    not_saturated = not_saturated && !(use_channel[channel] && saturated(data[CH_EYE+channel], num_samples));
	  if(not_saturated)
	    {
	    realft(veog_data-1, num_samples, 1);
	    for(sample = 1; sample != num_samples/2; sample++)
	      {
	      veog_data[2*sample+1] *= -1.0;
	      summed_eog_power[sample-1] += veog_data[2*sample]*veog_data[2*sample]+veog_data[2*sample+1]*veog_data[2*sample+1];
	      }
	    for(channel = 0; channel != num_channels; channel++)
	      if(use_channel[channel])
		{
		realft(data[CH_EYE+channel]-1, num_samples, 1);
		for(sample = 1; sample != num_samples/2; sample++)
		  {
		  ComplexMultiply(data[CH_EYE+channel]+2*sample, veog_data+2*sample, term);
		  blink_filter[channel][2*sample-2] += term[0];
		  blink_filter[channel][2*sample-1] += term[1];
		  }
		}
	    num_eye_movements++;
	    }
	  }
	}
      time += EP_INC*num_samples;
      }
    }
  while(event_number != num_events);
/*normalise:*/
  for(sample = 0; sample != num_samples/2-1; sample++)
    {
    for(channel = 0; channel != num_channels; channel++)
      if(use_channel[channel])
	{
	blink_filter[channel][2*sample] /= summed_eog_power[sample];
	blink_filter[channel][2*sample+1] /= summed_eog_power[sample];
	}
    if(num_quiet_epochs != 0)
      avg_quiet_eog_power[sample] /= num_quiet_epochs;
    }
  }

/*Search the log file for a response that occurs between RS_BEGIN ms
  RS_END ms after presentation of the target.*/
Boolean response()
  {
  static int cached_event = 0;
  static Boolean cached_result;
  long target_time;
  if((input_mode==EPL)||(input_mode==NeuroScan))
    {
    if(cached_event != event_number)
      {
      cached_event = event_number;
      target_time = event_time;
      while((event_number != num_events) && (event_time < target_time+msecs_to_ticks(RS_BEGIN)))
	next_event();
      while((event_number != num_events) && (event_time <= target_time+msecs_to_ticks(RS_END)) && (event_code != EV_RESPONSE))
	next_event();
      cached_result = (event_time <= target_time+msecs_to_ticks(RS_END)) && (event_time >= target_time+msecs_to_ticks(RS_BEGIN));
      reaction_time = cached_result? ticks_to_msecs(event_time - target_time)/1000.0: 0.0; /*used by print_epoch to generate NeuroScan output*/
      goto_event(cached_event); /*cover our tracks*/
      }
    }
  else
    /*input_mode==NeuroScanEpoched*/
    {
    cached_result = (reaction_time <= RS_END/1000.0) && (reaction_time >= RS_BEGIN/1000.0);
    /*reaction_time has already been set by next_event()*/
    }
  return(cached_result);
  }

void load_presample(data)
float data[][RE_PSLEN];
  {
  register int channel, t;
  for(t = 0; t != presample_length; t++)
    for(channel = -CH_EYE; channel != num_channels; channel++)
      data[channel][t] = eeg(event_time-presample_length+t, channel);
  }

/*Check for saturation of eye channels during the presample.  If artefact
  correction is off, always check for eye artefacts, too.*/
Boolean presample_clear(data)
float data[][RE_PSLEN];
  {
  Boolean result;
  result = (
	(correct_artefacts || (ppa(data[CH_LOWER_EYE], presample_length) < rejection_threshold))
     && !saturated(data[CH_LOWER_EYE], presample_length)
     && (correct_artefacts || (ppa(data[CH_UPPER_EYE], presample_length) < rejection_threshold))
     && !saturated(data[CH_UPPER_EYE], presample_length));
  if(!result) num_rejects++;
  return(result);
  }

/*Apply a Woestenburg-Wiener filter to the EEG.  Use the average quiet EOG
  power spectrum to estimate the portion of the EOG at each frequency that is
  due to eye activity, as distinct from the portion due to contamination by EEG
  (Wiener filtering).  Then apply the complex regression coefficient to this EOG
  fraction and subtract the product from the EEG (Woestenburg filtering).
  Woestenburg's F-test is not included in this implementation.
  Any negative Wiener coefficients are zeroed.*/
void filter(data)
float data[][EP_LEN/EP_INC];
  {
  register int sample, channel;
  float artefact[2],
	Wiener_coefficients[(EP_LEN/EP_INC)/2-1],
	veog_data[EP_LEN/EP_INC];
  for(sample = 0; sample != num_samples; sample++)
/*compute the VEOG*/
    veog_data[sample] = data[CH_UPPER_EYE][sample] - data[CH_LOWER_EYE][sample];
/*transform the VEOG to the frequency domain*/
  realft(veog_data-1, num_samples, 1);
/*transform & filter each used channel of EEG*/
  for(sample = 1; sample != num_samples/2; sample++)
    {
    Wiener_coefficients[sample-1] = num_quiet_epochs==0? 1.0: /*skip Wiener stage if no quiet epochs to base it on*/
      1.0 - avg_quiet_eog_power[sample-1]/(veog_data[2*sample]*veog_data[2*sample]+veog_data[2*sample+1]*veog_data[2*sample+1]);
    if(Wiener_coefficients[sample-1] > 0.0)
      {
      summed_Wiener_coefficients[sample-1] += Wiener_coefficients[sample-1];
      summed_squared_Wiener_coefficients[sample-1] += Wiener_coefficients[sample-1]*Wiener_coefficients[sample-1];
      }
    }
  for(channel = 0; channel != num_channels; channel++)
    if(use_channel[channel])
      {
      realft(data[channel]-1, num_samples, 1);
      for(sample = 1; sample != num_samples/2; sample++)
	if(Wiener_coefficients[sample-1] > 0.0)
	  {
	  ComplexMultiply(blink_filter[channel]+2*(sample-1), veog_data+2*sample, artefact);
	  data[channel][2*sample] -= artefact[0]*Wiener_coefficients[sample-1];
	  data[channel][2*sample+1] -= artefact[1]*Wiener_coefficients[sample-1];
	  }
      }
  }

/*Check for saturation of eye channels and of all used EEG channels.  If
  artefact correction is off, also check for blinks on the eye channels.  If
  artefact correction is on, filter out blinks.  If the time_domain flag is set,
  leave the EEG data in the time domain; otherwise leave it in the frequency
  domain.  Correct baseline to the average level of the presample.*/
Boolean sample_clear(data, presample)
float data[][EP_LEN/EP_INC], presample[][RE_PSLEN];
  {
  register int channel, sample;
  float baseline_correction;
  Boolean result;
  result = true;
  for(channel = 0; channel != CH_EYE; channel++)
    result = result && !saturated(data[channel], num_samples) && (correct_artefacts || (ppa(data[channel], num_samples) < rejection_threshold));
  for(channel = 0; channel != num_channels; channel++)
    result = result && !(use_channel[channel] && saturated(data[CH_EYE+channel], num_samples));
  if(result)
    {
    if(correct_artefacts)
      {
      filter(data+CH_EYE);
      if(time_domain)
	{
	for(channel = 0; channel != num_channels; channel++)
	  if(use_channel[channel])
	    {
	    data[CH_EYE+channel][0] = avg(presample[CH_EYE+channel], presample_length);
	    realft(data[CH_EYE+channel]-1, num_samples, -1);
	    for(sample = 0; sample != num_samples; sample++)
	      data[CH_EYE+channel][sample] /= num_samples/2;
	    }
	}
      }
    else
      {
      if(time_domain)
	{
	for(channel = CH_EYE; channel != CH_EYE+num_channels; channel++)
	  if(use_channel[channel-CH_EYE])
	    {
	    baseline_correction = avg(presample[channel], presample_length);
	    for(sample = 0; sample != num_samples; sample++)
	      data[channel][sample] -= baseline_correction;
	    }
	}
      else
	for(channel = 0; channel != num_channels; channel++)
	  if(use_channel[channel])
	    realft(data[CH_EYE+channel]-1, num_samples, 1);
      }
    num_accepts++;
    }
  else
    num_rejects++;
  return(result);
  }

/*If the time_domain flag is set, print the epoch as num_samples_in_output samples.
  Otherwise print it as num_samples/2-1 amplitudes followed by num_samples/2-1
  phase angles (the DC offset and the Nyquist frequency don't really matter).*/
void print_epoch(data)
float data[][EP_LEN/EP_INC];
  {
  register int sample, channel;
  SWEEP_HEAD sweep_header;
  short buf;
  sweep_header.accept = 1; /*rejected epochs have already been thrown out*/
  sweep_header.ttype = swap_bytes(event_code);
  sweep_header.correct = swap_bytes(reaction_time != 0.0);
  sweep_header.rt = reverse_float(reaction_time);
  sweep_header.response = swap_bytes(EV_RESPONSE);
  fwrite(&sweep_header, sizeof(sweep_header), 1, EEG);
  if(time_domain)
    {
    for(sample = 0; sample != num_samples_in_output; sample++)
      for(channel = -CH_EYE; channel != num_channels; channel++)
      /*in the time domain, include eye channels in the output*/
	if((channel < 0) || use_channel[channel])
	  {
	  float temp;
	  temp = VOLTAGE_MULTIPLIER*data[channel][sample];
	  buf = swap_bytes(temp>32767.0? 32767: temp<-32768.0? -32768: (short)temp);
	  fwrite(&buf, sizeof(short), 1, EEG);
	  }
    }
  else
    {
    /*frequency domain*/
    for(sample = 1; sample != num_samples/2; sample++)
      for(channel = 0; channel != num_channels; channel++)
	if(use_channel[channel])
	  {
	  buf = swap_bytes((short)(VOLTAGE_MULTIPLIER*AMPLITUDE_SCALE*hypot(data[channel][2*sample], data[channel][2*sample+1])));
	  fwrite(&buf, sizeof(short), 1, EEG);
	  }
    for(sample = 1; sample != num_samples/2; sample++)
      for(channel = 0; channel != num_channels; channel++)
	if(use_channel[channel])
	  {
	  buf = swap_bytes((short)(VOLTAGE_MULTIPLIER*AMPLITUDE_SCALE*PHASE_SCALE*atan2(data[channel][2*sample+1], data[channel][2*sample])));
	  fwrite(&buf, sizeof(short), 1, EEG);
	  }
    }
  }

/*Split the continuous record into epochs locked to stimulus presentation.*/
void epoch()
  {
  extern float data[CH_EYE+CH_MAX][EP_LEN/EP_INC];
  extern float presample[CH_EYE+CH_MAX][RE_PSLEN];
  num_rejects = 0;
  num_accepts = 0;
  while(event_number != num_events)
    {
    next_event();
    printf("\rEVENT %d OF %d", event_number, num_events);
    if((event_code <= EV_STIMULUS) && (event_code >= EV_LOW_STIMULUS))
      /*this is a stimulus, and not a response or pause or other event*/
      {
      load_presample(presample+CH_EYE);
	if(presample_clear(presample+CH_EYE))
	  {
	  load_eeg(data+CH_EYE, event_time);
	  if(sample_clear(data, presample))
	    {
	    response();
	    print_epoch(data+CH_EYE);
	    }
	  }
      }
    }
  }

void main(argc, argv)
int argc;
char *argv[];
  {
  check_static_definitions();
  process_options(argc, argv);
  open_files(argv[argc-1]);
  if(input_mode == EPL) read_header();
  check_dynamic_definitions();
  ask_channels();
  if(fseek(EEG, (long)(sizeof(erp)+(time_domain? CH_EYE+num_channels_used: num_channels_used)*sizeof(ELECTLOC)), SEEK_SET) != 0)
    fatal_error(".EEG seek failed");
  printf("calibrating..."); fflush(stdout);
  if(input_mode == EPL) calibrate(); else num_cal_pulses = 0;
  putchar('\n');
  if(correct_artefacts)
    {
    printf("computing filters...\n"); fflush(stdout);
    find_regression_coefficients();
    putchar('\n');
    }
  printf("epoching...\n"); fflush(stdout);
  goto_event(num_cal_pulses);
  epoch();
  putchar('\n');
  close_files();
  generate_report(argv[argc-1]);
  exit(0);
  }