/*detrend.c - detrend a GnuroScan averaged 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 <stdlib.h>
#include <string.h>
#ifdef MSDOS
#  define strcasecmp strcmpi
#endif
#include <errno.h>
#include "intel.h"
#include "gnuro.h"

#define MAX_WINDOW_LENGTH 256
#define CH_MAX 16

void main(argc, argv)
int argc;
char **argv;
  {
  register float sum;
  register int dx, i, channel, window_length, num_samples, num_channels;
  int srate, variance_flag;
  FILE *in, *out;
  int detrend_flag[CH_MAX];
  float buf[MAX_WINDOW_LENGTH];
  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]);
    fclose(in);
    exit(errno);
    }
  for(i = 0; i != S_pnts_offset; i++)
    putc(getc(in), out);
  num_samples = getc(in);
  putc(num_samples, out);
  i = getc(in);
  putc(i, out);
  num_samples |= i<<8;
  for(i = 0; i != S_nchans_offset-(S_pnts_offset+2); i++)
    putc(getc(in), out);
  num_channels = getc(in);
  putc(num_channels, out);
  i = getc(in);
  putc(i, out);
  num_channels |= i<<8;
  for(i = 0; i != S_variance_offset-(S_nchans_offset+2); i++)
    putc(getc(in), out);
  variance_flag = getc(in);
  putc(variance_flag, out);
  for(i = 0; i != S_rate_offset-(S_variance_offset+1); i++)
    putc(getc(in), out);
  srate = getc(in);
  putc(srate, out);
  i = getc(in);
  putc(i, out);
  srate |= i<<8;
  for(i = 0; i != packed_sizeof_SETUP-(S_rate_offset+2); i++)
    putc(getc(in), out);
  window_length = (atoi(argv[3])*srate+500)/1000;
  if(window_length > MAX_WINDOW_LENGTH)
    {
    fclose(out);
    fclose(in);
    fprintf(stderr, "Recompile %s with MAX_WINDOW_LENGTH=%d\n", *argv, window_length);
    exit(1);
    }
  if(window_length&1)
  /*make sure that window_length is even, for symmetry*/
    window_length--;
  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++)
    {
    char channel_label[EL_lab_length];
    for(i = 0; i != EL_lab_offset; i++)
      putc(getc(in), out);
    fread(channel_label, 1, EL_lab_length, in);
    detrend_flag[channel] = ((strcasecmp("heog", channel_label) != 0)
      &&(strcasecmp("veog", channel_label) != 0)
      &&(strcasecmp("blink", channel_label) != 0)
      &&(strcasecmp("upe", channel_label) != 0)
      &&(strcasecmp("loe", channel_label) != 0)
      &&(strcasecmp("lc", channel_label) != 0)
      &&(strcasecmp("rc", channel_label) != 0));
    fwrite(channel_label, 1, EL_lab_length, out);
    for(i = EL_lab_offset+EL_lab_length; i != packed_sizeof_ELECTLOC; i++)
      putc(getc(in), out);
    }
  for(channel = 0; channel != num_channels; channel++)
    {
  /*copy the sweep header*/
    for(i = 0; i != 5; i++)
      putc(getc(in), out);
  /*if this is an eye channel, do not detrend it; rather, just copy it as is*/
    if(detrend_flag[channel])
      {
    /*copy the left-hand half-window of boundary points*/
      for(i = 0; i != window_length/2-1; i++)
	{
	buf[i] = read_Intel_float(in);
	write_Intel_float(buf[i], out);
	}
    /*accumulate the rest of the initial contents of the window*/
      for(i = window_length/2-1; i != window_length-1; i++)
	buf[i] = read_Intel_float(in);
  /*
  An average of the effects of all the linear trends in windows overlapping the
  sample Y(i) is added to Y(i):

                 w/2-1
                  ___
                  \             dx
                  /  (Y(i-dx) + --- * (Y(i-dx+w/2)-Y(i-dx)))
                  ---           w/2
                 dx=1
  Y'(i) = Y(i) - -------------------------------------------
                                   w/2-1
  */
      for(i = window_length-1; i != num_samples; i++)
	{
	sum = 0.0;
	for(dx = 1; dx != window_length/2; dx++)
	  {
	  sum += buf[(i-window_length/2-dx)%window_length] + (dx/(window_length/2))*(buf[(i-dx)%window_length]-buf[(i-window_length/2-dx)%window_length]);
	  }
	write_Intel_float(buf[(i-window_length/2)%window_length]-sum/(window_length/2-1), out);
	buf[i%window_length] = read_Intel_float(in);
	}
    /*copy the right-hand half-window of boundary points*/
      for(i = num_samples-window_length/2; i != num_samples; i++)
	write_Intel_float(buf[i%window_length], out);
      }
    else
      for(i = 0; i != num_samples*sizeof(float); i++)
	putc(getc(in), out);
    }
  if(variance_flag)
    for(channel = 0; channel != num_channels; channel++)
      for(i = 0; i != num_samples*sizeof(float); i++)
	putc(getc(in), out);
  fclose(out);
  fclose(in);
  exit(0);
  }
