/*extreme.c - compute the range, in a moving window, of an averaged GnuroScan
  EEG data file.
  Copyright (c) 1996 Matthew Belmonte

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

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

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

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include "intel.h"
#include "fmedian.h"
#include "gnuro.h"
#ifdef NOFLOATS
#  define fabsf(x) ((float)(fabs((float)x)))
#endif

#define MAX_WINDOW_LENGTH 256

int read_and_copy_swapped_word(in, out)
FILE *in, *out;
  {
  register int word;
  int c;
  word = getc(in);
  putc(word, out);
  c = getc(in);
  putc(c, out);
  word |= c<<8;
  return(word);
  }

float extract_extreme(t)
btree *t;
  {
  register float min, max;
  min = extract_minimum(t);
  max = extract_maximum(t);
  return((fabsf(min) > fabsf(max))? min: max);
  }

void main(argc, argv)
int argc;
char **argv;
  {
  register int i, sample, channel, window_length, num_samples, num_channels;
  int srate, window_length_ms;
  FILE *in, *out;
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc != 4)
    {
    fprintf(stderr, "usage: %s <infile.avg> <outfile.avg> <window length (ms)>\n", *argv);
    exit(1);
    }
  if((in = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
  if((out = fopen(argv[2],
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    )) == NULL)
    {
    perror(argv[2]);
    exit(errno);
    }
  for(i = 0; i != S_pnts_offset; i++)
    putc(getc(in), out);
  num_samples = read_and_copy_swapped_word(in, out);
  for(i = S_pnts_offset+2; i != S_nchans_offset; i++)
    putc(getc(in), out);
  num_channels = read_and_copy_swapped_word(in, out);
  for(i = S_nchans_offset+2; i != S_variance_offset; i++)
    putc(getc(in), out); 
  getc(in); putc(0, out); /*variance data not included*/
  for(i = S_variance_offset+1; i != S_rate_offset; i++)
    putc(getc(in), out);
  srate = read_and_copy_swapped_word(in, out);
  for(i = S_rate_offset+2; i != packed_sizeof_SETUP+num_channels*packed_sizeof_ELECTLOC; i++)
    putc(getc(in), out);
  window_length_ms = atoi(argv[3]);
  window_length = (window_length_ms*srate+500)/1000;
  if(window_length > MAX_WINDOW_LENGTH)
    {
    fclose(out);
    unlink(argv[2]);
    fclose(in);
    fprintf(stderr, "Recompile %s with MAX_WINDOW_LENGTH=%d\n", *argv, window_length);
    exit(1);
    }
  printf("%d points at %dHz on %d channels - using a %d-point window\n", num_samples, srate, num_channels, window_length);
  if(window_length > num_samples)
    window_length = num_samples;
  for(channel = 0; channel != num_channels; channel++)
    {
    btree *t;
  /*5-byte sweep header*/
    for(i = 0; i != 5; i++)
      putc(getc(in), out);
    t = create_btree(window_length);
  /*accumulate the initial contents of the window*/
    for(sample = 0; sample != (window_length+1)/2; sample++)
      insert_newest(read_Intel_float(in), t);
    for(sample = (window_length+1)/2; sample != window_length; sample++)
      {
      insert_newest(read_Intel_float(in), t);
      write_Intel_float(extract_extreme(t), out);
      }
  /*slide the range window across the average time series*/
    for(sample = window_length; sample != num_samples; sample++)
      {
      delete_oldest(t);
      insert_newest(read_Intel_float(in), t);
      write_Intel_float(extract_extreme(t), out);
      }
    for(sample = 0; sample != (window_length+1)/2; sample++)
      {
      delete_oldest(t);
      write_Intel_float(extract_extreme(t), out);
      }
    destroy_btree(t);
    }
  fclose(out);
  fclose(in);
  exit(0);
  }
