// file: $(NEDC_NFC)/class/cpp/Edf/edf_06.cc
//
// This file contains signal processing and interpolation methods.
//

// Revision History:
//
// 20240307 (JP): refactored code to expand set/get methods
//

// local include files
//
#include "Edf.h"

//*****************************************************************************
//
// public methods: signal processing and interpolation methods
//
//*****************************************************************************

// method: copy_signal
//
// arguments:
//  VVectorDouble& sigo: selected signal (output)
//  VVectorDouble& sigi: signal from file (input)
//
// return: a boolean value indicating status
//
// This method copies one signal data structure to another.
//
bool Edf::copy_signal(VVectorDouble& sigo_a, VVectorDouble& sigi_a) {

  // resize the output
  //
  long nchan = sigi_a.size();
  Edf::resize(sigo_a, nchan, false);

  // loop over each channel
  //
  for (long i = 0; i < nchan; i++) {
    Edf::resize(sigo_a[i], sigi_a[i].size(), false);
    sigo_a[i].assign(sigi_a[i]);
  }

  // exit gracefully
  //
  return true;
}

// method: copy_signal
//
// arguments:
//  VVectorDouble& sigo: selected signal (output)
//  long& nlo: the number of labels (output)
//  cha** labelso: the label array (output)
//  VVectorDouble& sigi:  signal from file (input)
//  long nli: the number of labels (input)
//  char** labelsi: the label array (input)
//
// return: a boolean value indicating status
//
// This method copies the signal data from input to output,
// and also copies the channel labels. The latter is done because
// this same routine is used for channel selection and montage processing.
// The arrays are passed as arguments because the source arrays are
// different for channel selection and montage processing.
//
bool Edf::copy_signal(VVectorDouble& sigo_a, long& nlo_a, char** labelso_a,
		      VVectorDouble& sigi_a, long nli_a, char** labelsi_a) {

  // display debug information
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: starting to copy signal\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }
  
  // create output space
  //
  nlo_a = sigi_a.size();
  Edf::resize(sigo_a, nlo_a, false);

  // save all the channels
  //
  for (long i = 0; i < nlo_a; i++) {

    // create space
    //
    long j_end = sigi_a[i].size();
    Edf::resize(sigo_a[i], j_end, false);

    // copy the data
    //
    for (long j = 0; j < j_end; j++) {
      sigo_a[i][j] = sigi_a[i][j];
    }
  }

  // display debug information
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: done copying signal\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
    fprintf(stdout, "%s (line: %d) %s: starting copying labels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }
  
  // copy the labels
  //  note that it is assumed that the labels array has
  //  had space allocated for it.
  //
  for (long i = 0; i < nlo_a; i++) {
    Edf::resize(labelso_a[i], strlen(labelsi_a[i]) + 1);
    strcpy(labelso_a[i], labelsi_a[i]);
  }

  // display debug information
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: done copying labels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }
  
  // exit gracefully
  //
  return true;
}

// method: set_interp_chans
//
// arguments:
//  long new_channels_a: the number of new channels to interpolate (input)
//
// return: a boolean value indicating status
//
// This method sets the number of channels to be interpolated.
//      
bool Edf::set_interp_chans(long new_channels_a) {
  num_clabels_d = new_channels_a;
  if (new_channels_a == 0) {
    return false;
  }
  else {
    return true;
  }
}

// method: add_interp_channel
//
// arguments:
//  VVectorDouble& sigo: new signal (output)
//  VVectorDouble& sigi: original signal (input)
//  VectorDouble& chani: new interpolated channels (input)
//  INTERPOLATE_OMODE int_omode: output mode (input)
//  
// return: a boolean value indicating status
//
// This method adds the interpolated channels in the original set of
// channels and updates the header information.
//
// Note that all information regarding the new channels and their adjacent
// channels are stored in protected variables once the parameter file is
// parsed. Therefore, those variables are used in this method to add the new
// channels.
//
// Note also that this method will preserve the original set of channels,
// or copy just the new channels depending on the output mode.
//
bool Edf::add_interp_channel(VVectorDouble& sigo_a, VVectorDouble& sigi_a, 
			     VVectorDouble& chani_a,
			     INTERPOLATE_OMODE int_omode_a) {
  
  // the vector containing the positions of the adjacent_channels:
  //  as the new channel is being interpolated regarding its adjacent
  //  channels, their position in the original signal should be
  //  known. These positions will be used to create the header
  //  specific information related to new channels. 
  //
  long* pos_adj_channels[num_clabels_d];

  // variable containing the old number of channels:
  //  as the interpolation mode outputs can change either to 
  //  new channel, adjacent channels, or all channels, the old 
  //  number of channels might change.
  // 
  long old_num_channels = 0;

  // variable conatining the new number of channels
  //  as the interpolation mode outputs can change either to 
  //  new channel, adjacent channels, or all channels, the new 
  //  number of channels might change as well.
  //
  long new_num_channels = 0;

  // create output space to store the new set of channels:
  //  these vectors will work as auxiliary variables to 
  //  update the header.
  //  
  char* new_chan_labels[num_clabels_d];
  char* new_chan_trans_type[num_clabels_d];
  char* new_chan_phys_dim[num_clabels_d];
  double new_chan_phys_min[num_clabels_d];
  double new_chan_phys_max[num_clabels_d];
  long new_chan_dig_min[num_clabels_d];
  long new_chan_dig_max[num_clabels_d];
  char* new_chan_prefilt[num_clabels_d];
  long new_chan_rec_size[num_clabels_d];
  //long new_sample_freqs[num_clabels_d];  

  // display debug information
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: beginning addition of new channels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }

  // loop over new channels numbers and take the position of its
  // adjacent channels one by one
  //
  for (long i = 0; i < num_clabels_d; i++) {

    // number of adjacent channels related to
    // new channel i
    //
    long nadj_chans = num_adj_chan_d[i];
    
    // resize columns of position to the number of adjacent channels
    // related to new_channel i
    //
    pos_adj_channels[i] = new long[nadj_chans];
    
    // loop over adjacent channels of new channel i
    //
    for (long j = 0; j < nadj_chans; j++) {

      // find the corresponding channel in the input signal
      //
      long pos;
      
      // match mode:
      //  smmode_d has the last value assigned to it. That means when
      //  interpolation is performed, the actual match_mode might
      //  be saved. 
      //  
      if ((pos = Edf::find_match(adj_chan_labels_d[i][j],
			         hdr_ghdi_nsig_rec_d, hdr_chan_labels_d,
			         smmode_d)) < 0) {
        fprintf(stdout, "Error: %s (line: %d) %s: no match for [%s]\n",
		__FILE__, __LINE__, __PRETTY_FUNCTION__,
	        adj_chan_labels_d[i][j]);

        return false;
      }

      // take position of adjacent channel j in the original signal 
      //
      pos_adj_channels[i][j] = pos;
    } 
  }  

  // case 1: preserve all channels including the interpolated channels:
  //  in this case, all the information regarding the original signal
  //  should be preserved, including its channels discrete points and
  //  the header specific information.
  //
  if (int_omode_a == INTOMODE_CONCAT) {

    // save number of channels preserved
    //
    old_num_channels = Edf::get_num_channels_file();
  }

  // case 2: preserve just the interpolated channel:
  //  in this case, just the new channel will be preserved and the specific 
  //  information stored will be from any of its adjacent channels.
  //
  else if (int_omode_a == INTOMODE_REPLACE) {

    // save number of channels preserved
    //
    old_num_channels = 0;
  }
  
  // case 3: unknown
  //
  else {
    fprintf(stdout, "%s (line: %d) %s: unknown output mode\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }  

  // new number of channels
  //
  new_num_channels = old_num_channels + num_clabels_d;
    
  // resize output signal to the final size
  //  the output signal rows might be resized to the number of new channels
  //  since the last row will represent the new channel. 
  //
  Edf::resize(sigo_a, new_num_channels, false);   

  // add new channel's information to new variables:  
  //  update header information related to the new channel.
  //  it copies the header information of one of the adjacent channels
  //  since the new channel is the interpolation of them.
  //
  for (long i = 0; i < num_clabels_d; i++) {    
    new_chan_labels[i] = new char[strlen(new_chan_labels_d[i]) + 1];
    strcpy(new_chan_labels[i], new_chan_labels_d[i]);

    new_chan_trans_type[i] =
    new char[strlen(hdr_chan_trans_types_d[pos_adj_channels[i][0]]) + 1];
    strcpy(new_chan_trans_type[i],
	   hdr_chan_trans_types_d[pos_adj_channels[i][0]]); 

    new_chan_phys_dim[i] =
    new char[strlen(hdr_chan_phys_dim_d[pos_adj_channels[i][0]]) + 1];
    strcpy(new_chan_phys_dim[i], hdr_chan_phys_dim_d[pos_adj_channels[i][0]]); 

    new_chan_phys_min[i] = hdr_chan_phys_min_d[pos_adj_channels[i][0]];
    new_chan_phys_max[i] = hdr_chan_phys_max_d[pos_adj_channels[i][0]];
    new_chan_dig_min[i] = hdr_chan_dig_min_d[pos_adj_channels[i][0]];
    new_chan_dig_max[i] = hdr_chan_dig_max_d[pos_adj_channels[i][0]];
    new_chan_dig_max[i] = hdr_chan_dig_max_d[pos_adj_channels[i][0]];

    new_chan_prefilt[i] =
    new char[strlen(hdr_chan_prefilt_d[pos_adj_channels[i][0]]) + 1];
    strcpy(new_chan_prefilt[i], hdr_chan_prefilt_d[pos_adj_channels[i][0]]);

    new_chan_rec_size[i] = hdr_chan_rec_size_d[pos_adj_channels[i][0]];
  }
    
  // save the signal data from the original channels
  //
  for (long i = 0; i < old_num_channels; i++) {
    // create output space
    //
    long j_end = sigi_a[i].size();
    Edf::resize(sigo_a[i], j_end, false);

    // copy the data
    //
    for (long j = 0; j < j_end; j++) {
      sigo_a[i][j] = sigi_a[i][j];
    }       
  }

  // save new channel's information:   
  //  update header information related to the new channel.
  //  it copies the header information of one of its  adjacent
  //  channels since the new channel is the interpolation of them.
  //
  for (long i = old_num_channels; i < new_num_channels; i++) {
    Edf::resize(hdr_chan_labels_d[i],
    	        strlen(new_chan_labels[i - old_num_channels]) + 1, false);
    strcpy(hdr_chan_labels_d[i],new_chan_labels[i - old_num_channels]);

    Edf::resize(hdr_chan_trans_types_d[i],
		strlen(new_chan_trans_type[i - old_num_channels]) + 1, false);
    strcpy(hdr_chan_trans_types_d[i],
	   new_chan_trans_type[i - old_num_channels]);

    Edf::resize(hdr_chan_phys_dim_d[i],
    	        strlen(new_chan_phys_dim[i - old_num_channels]) + 1, false);
    strcpy(hdr_chan_phys_dim_d[i],new_chan_phys_dim[i - old_num_channels]);

    hdr_chan_phys_min_d[i] = new_chan_phys_min[i - old_num_channels]; 
    hdr_chan_phys_max_d[i] = new_chan_phys_max[i - old_num_channels]; 
    hdr_chan_dig_min_d[i] = new_chan_dig_min[i - old_num_channels]; 
    hdr_chan_dig_max_d[i] = new_chan_dig_max[i - old_num_channels]; 

    Edf::resize(hdr_chan_prefilt_d[i],
    	        strlen(new_chan_prefilt[i - old_num_channels]) + 1, false);
    strcpy(hdr_chan_prefilt_d[i],new_chan_prefilt[i - old_num_channels]);

    hdr_chan_rec_size_d[i] = new_chan_rec_size[i - old_num_channels];

    // add the new channel's signal data:
    // create output space to hold size of each new channel
    //
    long j_end = chani_a[i - old_num_channels].size();
    Edf::resize(sigo_a[i], j_end, false);
    
    // copy the data
    //
    for (long j = 0; j < j_end; j++) {
      sigo_a[i][j] = chani_a[i - old_num_channels][j];
    }
  }

  // clean up memory
  //
  Edf::cleanup(new_chan_labels, num_clabels_d);
  Edf::cleanup(new_chan_trans_type, num_clabels_d);
  Edf::cleanup(new_chan_phys_dim, num_clabels_d);
  Edf::cleanup(new_chan_prefilt, num_clabels_d);

  // loop over number of new channels and clean memory allocated
  // to position of adjacent channels 
  //
  for (long i = 0; i < num_clabels_d; i++) {
    delete [] pos_adj_channels[i];
  }

  // update the number of channels and the header size 
  //
  hdr_ghdi_nsig_rec_d = new_num_channels;
  hdr_ghdi_hsize_d = compute_header_size(hdr_ghdi_nsig_rec_d);

  // display debug information
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: done with addition of new channels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }
  
  // exit gracefully
  //
  return true;
}

// method: interpolate
//
// arguments:
//  VVectorDouble& sigo: new set of channels (input/output)
//  MATCH_MODE match_mode: match mode (input)
//  INTERPOLATE_MODE mode: interpolate mode (input)
//  INTERPOLATE_OMODE omode: output mode (input)
//  
// return: a boolean value indicating status
//
// This method branches to the interpolation mode selected by the
// parameter file and add the interpolated channels signal.
//
// Note that this method will update the header specific information.
//
bool Edf::interpolate(VVectorDouble& sigo_a,
		      MATCH_MODE match_mode_a,
		      INTERPOLATE_MODE mode_a,
		      INTERPOLATE_OMODE omode_a) {

  // declare local variables
  //
  VVectorDouble sig_t;
  bool status = false;

  // save match_mode
  //
  smmode_d = match_mode_a;
  
  // variable to store the new_channels discrete points
  //
  VVectorDouble new_channels;

  // display a debug message
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: beginning interpolation\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }

  // branch on the interpolate mode
  //
  // case 1: interpolation by average
  //
  if (mode_a == INTMODE_AVERAGE) {

    // display debug information
    //
    if (dbgl_d > Dbgl::BRIEF) {
      fprintf(stdout, "%s (line: %d) %s: average mode (%lu)\n",
	      __FILE__, __LINE__, __PRETTY_FUNCTION__,
	      (long)mode_a);
    }
    if (!(status = interpolate_average(new_channels, sigo_a))) {
      return status;
    }
  }
    
  // case 2: unknown
  //
  else {
    fprintf(stdout, "%s (line: %d) %s: unknown mode \n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }

  // return the new signal
  //
  // copy the original signal:
  //  make a copy of the original signal to add the
  //  interpolated channels into the new signal.
  //
  if (!(status = Edf::copy_signal(sig_t, sigo_a))) {
    return status;
  }
  
  // add interpolated channels into the original set of channels
  //
  if (!(status = Edf::add_interp_channel(sigo_a, sig_t,
					 new_channels, omode_a))) {
    return status;
  }    

  // display a debug message
  //
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: original signal: %ld channels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__,
	    (long)sig_t.size());
    fprintf(stdout, "%s (line: %d) %s: new signal: %ld channels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__,
	    (long)sigo_a.size());
  }
  if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: done interpolating channels\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }
  
  // exit gracefully
  //
  return status;
}

// method: test_signal
//
// arguments:
//  VVectorDouble& sig: the EEG signal data (input/output)
//  const char* mode: specifies the type of signal (input)
//
// return: a boolean value indicating status
//
// This method is provided as a simple way to generate a signal for
// debugging. It must be called after read() because the signal
// an all its parameters must exist. It is not intended for general
// use. It was originally developed to debug the front end code. It
// replaces each channel with the exact same signal (e.g., sinewaves).
//
bool Edf::test_signal(VVectorDouble& sig_a, const char* mode_a) {
 
  // declare local variables
  //
  bool status = true;

  // create test signal only if parameter is not (null)
  //
  if (strcmp(mode_a, "(null)") != 0) {

   // display debugging information
   //
   if (dbgl_d > Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: creating signal [%s]\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__,
	    mode_a);
   }

   // loop for all channels
   //
   long nc = sig_a.size();
   double fs_inv = 1.0 / hdr_sample_frequency_d;

   for (long i = 0; i < nc; i++) {
     
     // loop over all samples
     //
     long ns = sig_a[i].size();

    for (long j = 0; j < ns; j++) {

      // case 1: single sinewave
      //
      if (strcmp(mode_a, "single_sinewave") == 0) {

	// compute a single sinewave
	//
	double t = (double)j * fs_inv;
	double arg = 2 * M_PI * EDF_TST_F4 * t;
	sig_a[i][j] = EDF_TST_A4 * sin(arg);
      }

      // case 2: three sinewaves
      //
      else if (strcmp(mode_a, "three_sinewaves") == 0) {

	// compute the arguments of the sine functions
	//
	double t = (double)j * fs_inv;
	double arg1 = 2 * M_PI * EDF_TST_F4 * t;
	double arg2 = 2 * M_PI * EDF_TST_F6 * t;
	double arg3 = 2 * M_PI * EDF_TST_F7 * t;

	// compute the sine and sum
	//
	sig_a[i][j] = EDF_TST_A4 * sin(arg1) +
	  EDF_TST_A6 * sin(arg2) + EDF_TST_A7 * sin(arg3);
      }

      // case 3: six sinewaves
      //
      else if (strcmp(mode_a, "six_sinewaves") == 0) {

	// compute the arguments of the sine functions
	//
	double t = (double)j * fs_inv;
	double arg1 = 2 * M_PI * EDF_TST_F1 * t;
	double arg2 = 2 * M_PI * EDF_TST_F2 * t;
	double arg3 = 2 * M_PI * EDF_TST_F3 * t;
	double arg4 = 2 * M_PI * EDF_TST_F4 * t;
	double arg5 = 2 * M_PI * EDF_TST_F5 * t;
	double arg6 = 2 * M_PI * EDF_TST_F6 * t;

	// compute the sine and sum
	//
	sig_a[i][j] = EDF_TST_A1 * sin(arg1) + EDF_TST_A2 * sin(arg2) +
	  EDF_TST_A3 * sin(arg3) + EDF_TST_A4 * sin(arg4) +
	  EDF_TST_A5 * sin(arg5) + EDF_TST_A6 * sin(arg6);
      }

      // case 4: rectified_sinewave
      //
      else if (strcmp(mode_a, "rectified_sinewave") == 0) {

	// compute the arguments of the sine function
	//
	double t = (double)j * fs_inv;
	double arg1 = 2 * M_PI * EDF_TST_F1 * t;
	double sum = sin(arg1);

	// compute the sine and sum
	//
	if (sum > 0) {
	  sig_a[i][j] = EDF_TST_A1 * sum;
	}
	else {
	  sig_a[i][j] = 0.0;
	}
      }

      // case 5: sine_plus_noise
      //
      else if (strcmp(mode_a, "sine_plus_noise") == 0) {

	// compute a single sinewave plus noise
	//
	double t = (double)j * fs_inv;
	double arg = 2 * M_PI * EDF_TST_F4 * t;
	sig_a[i][j] = EDF_TST_A4 * (sin(arg) + drand48());
      }

      // else: error
      //
      else {
	status = false;
      }

      // display debugging information
      //
      if ((dbgl_d > Dbgl::DETAILED) && (j < DEF_DBG_NS)) {
	fprintf(stdout, "%s (line: %d) %s: sig[%ld][%ld] = %f\n",
		__FILE__, __LINE__, __PRETTY_FUNCTION__,
		i, j, sig_a[i][j]);
      }
    }
  }

  // display debugging information
  //
  if (dbgl_d >= Dbgl::BRIEF) {
    fprintf(stdout, "%s (line: %d) %s: done creating signal\n",
	    __FILE__, __LINE__, __PRETTY_FUNCTION__);
  }

  // exit gracefully
  //
  return status;
 }

 else {
   
   //exit gracefully
   //
   return status;
 }
}

//
// end of file