/*lderiv.c - construct a linear derivation of a GnuroScan continuous 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 <string.h>
#include <math.h>
#include <errno.h>
#include "intel.h"
#include "gnuro.h"

#define CH_MAX 16

static int channel_number(name, channel_hdr, num_channels)
char *name, channel_hdr[][packed_sizeof_ELECTLOC];
int num_channels;
  {
  register int len, chn;
  chn = 0;
  do
    {
    for(len = 0; (len != 10) && (channel_hdr[chn][len]==name[len]) && (name[len] != '\0'); len++)
      ;
    chn++;
    }
  while((chn != num_channels) && (len != 10) && !((channel_hdr[chn-1][len]==name[len]) && (name[len]=='\0')));
  return(((len == 10) || ((channel_hdr[chn-1][len]==name[len]) && name[len]=='\0'))? chn-1: -1);
  }

void main(argc, argv)
int argc;
char *argv[];
  {
  register int src_channel, obj_channel;
  register long sample;
  int num_src_channels, num_obj_channels;
  long num_samples, alt_num_samples, num_events;
  FILE *in, *out, *ldr;
  int channel[CH_MAX];
  short eeg[CH_MAX];
  int src_baseline[CH_MAX];
  double max_abs_eeg[CH_MAX], obj_baseline[CH_MAX];
  double trans_matrix[CH_MAX][CH_MAX];
  char erp[packed_sizeof_SETUP], channel_hdr[CH_MAX][packed_sizeof_ELECTLOC];
/*re-used storage:*/
#define acc src_channel
#define teeg_type obj_channel
#define event sample
#define evtpos num_events
  printf("Copyright (c) 1996 Matthew Belmonte <mkb4@Cornell.edu>.  Please cite.\n");
  if(argc != 4)
    {
    fprintf(stderr, "usage: %s <spec.ldr> <source.cnt> <derived.cnt>\n", *argv);
    exit(1);
    }
/*open the .ldr file*/
  if((ldr = fopen(argv[1],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[1]);
    exit(errno);
    }
/*get the column and row dimensions of the transformation matrix*/
  if(fscanf(ldr, "%d %d", &num_obj_channels, &num_src_channels) != 2)
    {
    fclose(ldr);
    fprintf(stderr, "Couldn't read dimensions from LDR file %s\n", argv[1]);
    exit(1);
    }
  if((num_src_channels <= 0) || (num_obj_channels <= 0))
    {
    fclose(ldr);
    fprintf(stderr, "bad number of channels in parameter file %s\n", argv[1]);
    exit(1);
    }
  if((num_src_channels > CH_MAX) || (num_obj_channels > CH_MAX))
    {
    fclose(ldr);
    fprintf(stderr, "Recompile %s with CH_MAX=%d\n", *argv, (num_src_channels > num_obj_channels)? num_src_channels: num_obj_channels);
    }
/*open the input .cnt file*/
  if((in = fopen(argv[2],
#ifdef MSDOS
    "rb"
#else
    "r"
#endif
    )) == NULL)
    {
    perror(argv[2]);
    exit(errno);
    }
/*read the file header*/
  fread(erp, 1, packed_sizeof_SETUP, in);
/*verify the number of channels in the input*/
  src_channel = S_nchannels(erp);
  if(src_channel != num_src_channels)
    {
    fclose(in);
    fprintf(stderr, "%d channels specified in %s, %d channels found in %s\n", num_src_channels, argv[1], src_channel, argv[2]);
    exit(1);
    }
/*read the channel headers*/
  for(src_channel = 0; src_channel != num_src_channels; src_channel++)
    fread(channel_hdr[src_channel], 1, packed_sizeof_ELECTLOC, in);
/*Make channel[] a mapping from matrix column indices to .cnt record numbers.
  Also, record the baseline value of each .cnt channel.*/
  for(src_channel = 0; src_channel != num_src_channels; src_channel++)
    {
    char ch_name[11];
    if(fscanf(ldr, "%s", ch_name) != 1)
      {
      fprintf(stderr, "bad format in LDR file %s\n", argv[1]);
      fclose(in);
      fclose(ldr);
      exit(1);
      }
    channel[src_channel] = channel_number(ch_name, channel_hdr, num_src_channels);
    if(channel[src_channel] == -1)
      {
      fclose(in);
      fclose(ldr);
      fprintf(stderr, "Couldn't find channel %s in file %s\n", ch_name, argv[2]);
      exit(1);
      }
    src_baseline[src_channel] = EL_baseline(channel_hdr[src_channel]);
    }
/*Using the channel[] mapping, set up a transformation matrix.  Also, record the
  new channel names.*/
  for(obj_channel = 0; obj_channel != num_obj_channels; obj_channel++)
    {
    char ch_name[11];
    if(fscanf(ldr, "%s", ch_name) != 1)
      {
      fprintf(stderr, "bad format in LDR file %s\n", argv[1]);
      fclose(in);
      fclose(ldr);
      exit(1);
      }
    strncpy(channel_hdr[obj_channel], ch_name, 9);
    channel_hdr[obj_channel][9] = '\0';
    for(src_channel = 0; src_channel != num_src_channels; src_channel++)
      if(fscanf(ldr, "%lf", &(trans_matrix[channel[src_channel]][obj_channel])) != 1)
	{
	fclose(in);
	fclose(ldr);
	fprintf(stderr, "couldn't read transformation coefficient [%d, %d] from file %s\n", obj_channel, src_channel, argv[1]);
	exit(1);
	}
    }
/*Done with LDR file.*/
  fclose(ldr);
/*Pull out # of samples and event table offset from file header.*/
  num_samples = S_NumSamples(erp);
  evtpos = S_EventTablePos(erp);
/*check that the number of samples makes sense---if not, attempt to correct it*/
  alt_num_samples = (evtpos-(long)(num_src_channels*packed_sizeof_ELECTLOC)-(long)packed_sizeof_SETUP)/(long)(num_src_channels*sizeof(short));
  if(alt_num_samples != num_samples)
    {
    fprintf(stderr, "Warning: %ld samples specified in header doesn't match %ld samples calculated from file size - ", num_samples, alt_num_samples);
    if((num_samples < S_rate(erp) && num_samples < alt_num_samples))
      {
      num_samples = alt_num_samples;
      put_S_NumSamples(erp, alt_num_samples);
      }
    fprintf(stderr, "using %ld\n", num_samples);
    }
/*For each source channel, scale the corresponding row of the transformation
  matrix by its calibration values.*/
  for(obj_channel = 0; obj_channel != num_obj_channels; obj_channel++)
    {
    for(src_channel = 0; src_channel != num_src_channels; src_channel++)
      trans_matrix[src_channel][obj_channel] *= EL_sensitivity(channel_hdr[src_channel])*EL_calib(channel_hdr[src_channel]);
    }
/*Apply the transformation matrix, but don't save the results.  Rather, keep
  only the maximum absolute value of EEG for each derived channel.  This will
  be used in the computation of virtual baselines and calibration factors that
  optimise the resolution of the digital record without exceeding the range of
  the short integers used to compute and to store the record.*/
  for(obj_channel = 0; obj_channel != num_obj_channels; obj_channel++)
    {
    max_abs_eeg[obj_channel] = -32768.0;
    obj_baseline[obj_channel] = 0;
    for(src_channel = 0; src_channel != num_src_channels; src_channel++)
      obj_baseline[obj_channel] += trans_matrix[src_channel][obj_channel]*src_baseline[src_channel];
    }
  for(sample = 0L; sample != num_samples; sample++)
    {
    for(src_channel = 0; src_channel != num_src_channels; src_channel++)
      eeg[src_channel] = read_Intel_short(in);
    for(obj_channel = 0; obj_channel != num_obj_channels; obj_channel++)
      {
      double new_eeg;
      new_eeg = 0.0;
      for(src_channel = 0; src_channel != num_src_channels; src_channel++)
	new_eeg += trans_matrix[src_channel][obj_channel]*eeg[src_channel];
      if(new_eeg > max_abs_eeg[obj_channel])
	max_abs_eeg[obj_channel] = new_eeg;
      else if(-new_eeg > max_abs_eeg[obj_channel])
	max_abs_eeg[obj_channel] = -new_eeg;
      new_eeg += obj_baseline[obj_channel];
      if(new_eeg > max_abs_eeg[obj_channel])
	max_abs_eeg[obj_channel] = new_eeg;
      else if(-new_eeg > max_abs_eeg[obj_channel])
	max_abs_eeg[obj_channel] = -new_eeg;
      }
    }
/*Copy new # of channels and new event table offset into file header.*/
  put_S_nchannels(erp, num_obj_channels);
  put_S_EventTablePos(erp, packed_sizeof_SETUP+num_obj_channels*(packed_sizeof_ELECTLOC+num_samples*sizeof(short)));
/*Open output file.*/
  if((out = fopen(argv[3],
#ifdef MSDOS
    "wb"
#else
    "w"
#endif
    )) == NULL)
    {
    perror(argv[3]);
    fclose(in);
    exit(errno);
    }
/*Write the file header.*/
  fwrite(erp, 1, packed_sizeof_SETUP, out);
/*Compute virtual baselines and calibration factors, store them in the channel
  headers, and output them.  Multiply the calibration factors into the
  transformation matrix.*/
  for(obj_channel = 0; obj_channel != num_obj_channels; obj_channel++)
    {
    double scale;
    scale = max_abs_eeg[obj_channel]/32767.0;
    for(src_channel = 0; src_channel != num_src_channels; src_channel++)
      trans_matrix[src_channel][obj_channel] /= scale;
    put_EL_baseline(channel_hdr[obj_channel], (short)(obj_baseline[obj_channel]/scale));
    put_EL_sensitivity(channel_hdr[obj_channel], 1.0);
    put_EL_calib(channel_hdr[obj_channel], scale);
    fwrite(channel_hdr[obj_channel], 1, packed_sizeof_ELECTLOC, out);
    }
/*Apply the transformation matrix.*/
  fseek(in, (long)(packed_sizeof_SETUP+num_src_channels*packed_sizeof_ELECTLOC), SEEK_SET);
  for(sample = 0L; sample != num_samples; sample++)
    {
    for(src_channel = 0; src_channel != num_src_channels; src_channel++)
      eeg[src_channel] = read_Intel_short(in);
    for(obj_channel = 0; obj_channel != num_obj_channels; obj_channel++)
      {
      float new_eeg;
      new_eeg = 0.0;
      for(src_channel = 0; src_channel != num_src_channels; src_channel++)
	new_eeg += trans_matrix[src_channel][obj_channel]*eeg[src_channel];
      write_Intel_short((short)(rint(new_eeg)), out);
      }
    }
/*Adjust event table offsets.*/
  fseek(in, evtpos, SEEK_SET);
  teeg_type = getc(in);
  if((teeg_type != 1) && (teeg_type != 2))
    fprintf(stderr, "warning: unknown TEEG type %d - assuming Type 1\n", teeg_type);
  putc(teeg_type, out);
  num_events = read_Intel_long(in);
  write_Intel_long(num_events, out);
  num_events /= (teeg_type>1? packed_sizeof_EVENT2: packed_sizeof_EVENT1);
  for(acc = 0; acc != 4; acc++)
    {
    getc(in);
    putc(0, out);
    }
  for(event = 0L; event != num_events; event++)
    {
    for(acc = 0; acc != packed_sizeof_EVENT1-sizeof(long); acc++)
      putc(getc(in), out);
    write_Intel_long(((read_Intel_long(in)-packed_sizeof_SETUP-num_src_channels*packed_sizeof_ELECTLOC)/num_src_channels+packed_sizeof_ELECTLOC)*num_obj_channels+packed_sizeof_SETUP, out);
    if(teeg_type > 1)
      for(acc = packed_sizeof_EVENT1; acc != packed_sizeof_EVENT2; acc++)
	putc(getc(in), out);
    }
  fclose(out);
  fclose(in);
  exit(0);
  }
