package nxm.dsp.prim;

import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.Primitive;

/** Port of the Fortran host primitive convolve from XMidas.

    <pre>
  ##############
  #            #
  #  CONVOLVE  #
  #            #
  ##############

  NAME:  CONVOLVE

  PURPOSE:
     To convolve the data in file &lt;in&gt; with that of file
     &lt;impulse&gt; to effect a linear time-invariant filter.

  DESCRIPTION:
     Get command line parameters
     Set-up necessary MIDAS files
     Read buffer of data
     Perform convolution
     Write buffer of data
     Post headers, close files and return to MIDAS

  CALLING SEQUENCE:

     CONVOLVE &lt;in&gt; &lt;impulse&gt; &lt;out&gt;

  INPUT FILES:
     MIDAS input file containing data
     MIDAS input file containing impulse response

  OUTPUT FILES:
     MIDAS file of results

  RESTRICTIONS:
     1.  There is no length limitation on &lt;impulse&gt;.
     2.  Both &lt;in&gt; and &lt;impulse&gt; must be array data or scalar data files.
         The output file &lt;out&gt; is always written as an array data file.
     3.  Data formats in &lt;in&gt; and &lt;impulse&gt; can be byte, integer, long, real,
         or double.  There is no requirement that these files be of the same
         format.  The format of &lt;out&gt; follows that of &lt;in&gt;.
     4.  &lt;out&gt; can overwrite &lt;in&gt;, but &lt;out&gt; cannot overwrite &lt;impulse&gt;.
     5.  There is no length limitation on &lt;in&gt;.

  MODIFICATIONS:
     Fixed Java port                              MRD     08/05
     Port to Java for NeXtMidas                   HJH     06/05
  </pre>

  @author Hoyt Heaton
  @author Marcus Deyoung
 */
@SuppressWarnings("unused")
public class convolve extends Primitive {

    private DataFile inputFile;
    private DataFile outputFile;
    private DataFile impulseFile;
    private Data inputData;
    private Data impulseData;
    private Data outputData;

    /** number of points to output this time */
    private int bufsize;
    /** number of taps in the impulse response */
    private int impulseTaps;
    /** this should be the same as impulse taps, but is the actual number read in */
    private int numImpulse;
    /** input work starting sample for this pass - the last of the last pass' input values (or zeros) are held in the lower locations */
    private int inpst;
    private float[] impulse;
    private float[] impulseWork;
    private float[] inputWork;

    /** default constructor */
    public convolve () {}

    public int open () {

      inputFile = MA.getDataFile ("IN", "1000", "S#",0);
      inputFile.open();

      impulseFile = MA.getDataFile ("IMPULSE", "1000", "S#", 0);
      impulseFile.open();
      impulseTaps = bufsize = (int)impulseFile.getSize();

      outputFile = MA.getDataFile ("OUT",inputFile,0);
      outputFile.setXStart( inputFile.getXStart() + inputFile.getXDelta()
              * (impulseFile.getXStart() / impulseFile.getXDelta()) );
      outputFile.setXDelta( inputFile.getXDelta() );
      outputFile.setXUnits( inputFile.getXUnits() );
      outputFile.setSize( inputFile.getSize() + impulseFile.getSize() - 1.0 );
      outputFile.open();

      if (impulseFile.filename.equals(outputFile.filename)) {
        M.error("Output file may not override impulse response");
      }

      //Get buffers for input and output
      inputData = inputFile.getDataBuffer(impulseTaps*4, FLOAT);
      outputData = outputFile.getDataBuffer(impulseTaps*4, FLOAT);

      impulseData = impulseFile.getDataBuffer(impulseTaps, FLOAT);

      //read in impulse file. If empty, exit.
      numImpulse = impulseFile.read(impulseData);
      inpst = numImpulse - 1;
      if (numImpulse <= 0) {
        return FINISH;
      }

      //Cast byte buffer to floats
      impulse = impulseData.castF(true);
      impulseWork = new float[impulseTaps];

      // reverse elements in impulse buffer
      // note: in the strange case we did not read in all the taps, the end of the impulse work buffer will have zeros
      for (int i=0; i<numImpulse; i++) {
          impulseWork[i] = impulse[(numImpulse - 1) - i];
      }

      impulseData.uncast(impulse, false);

      // Set up working input array (stores data from previous iteration
      //  of convolve).
      inputWork = new float[impulseTaps*4];
      return NORMAL;
    }

    public int process () {
      int numInput = inputFile.read(inputData, numImpulse);

      if (numInput == 0) { return NOOP; }

      // mrd -- need to run through conv last time if finished with the inputfile
      //if ( numInput < 0) {
      //    return FINISH;
      //}

      // mrd -- conv expects minimum of 0,
      // uses  ( bufsize = numInput + inpst; ) and will be off by 1 if
      // numInput is -1
      if(numInput < 0 ) { numInput = 0; }

      conv(inputData, outputData, numInput, numImpulse);

      outputFile.write(outputData, bufsize);

      // mrd -- done with inputFile, finish up
      if ( numInput < numImpulse) {return FINISH;}
      else                        {return (NORMAL);}

    }

    public int close () {
      inputFile.close();
      impulseFile.close();
      outputFile.close();
      return (NORMAL);
    }

    /** perform convolution using the reversed impulse array ({@link #impulseWork}). */
    private void conv (Data inputData, Data outputData,  int numInput, int numImpulse) {
      // Get arrays of floats from input/output data buffers to
      //   use during calculations.
      float[] input = inputData.castF(true);
      float[] output = outputData.castF(false);

      // Load data read from input file into the proper location in the
      // working input array
      if (numImpulse + numInput > inputWork.length) {
        M.error("Working input array not large enough.\n" +
                "numImpulse = " + numImpulse + "\n" +
                "numImput = " + numInput + "\n" +
                "inputWork.length = " + inputWork.length);
      }
      if (numInput > input.length) {
        M.error("input array(" + input.length + ") returned by cast is"
                + "smaller than numRead(" + numInput + ")" );
      }

      int j=0;
      for (int i=inpst; i<inpst + numInput; i++) {
        inputWork[i] = input[j++];
      }


      //if numInput = bufsize, read a complete chunk, so we are in the middle
      if ( numInput == bufsize  ) {
        //convolve data
        for (int i=0; i<bufsize; i++) {
          output[i] = 0.0f;
          for (j=0; j<numImpulse; j++) {
            output[i] += inputWork[i+j] * impulseWork[j];
          }
        }

        //move last numImpulse-1 points down to beginning of inputWork
        for (int i=0; i < (numImpulse-1); i++) {
          inputWork[i] = inputWork[i+bufsize];
        }
      } else {
        //otherwise, (numInput<bufsize), incomplete chunk, reached the end of the file
        //since final block, zero out top of inputwork array.
        bufsize = numInput + inpst;
        for (int i=bufsize; i< (bufsize + numImpulse - 1); i++) {
          inputWork[i] = 0.0f;
        }

        //convolve data
        for (int i=0; i < bufsize; i++) {
          output[i] = 0.0f;
          for (j=0; j<numImpulse; j++) {
            output[i] += inputWork[i+j] * impulseWork[j];
          }
        }
      }

      inputData.uncast(input, false);
      outputData.uncast(output, true);
    }
}
