package nxm.dsp.prim;

import nxm.sys.inc.Units;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.DataOp;
import nxm.sys.lib.MidasException;
import nxm.sys.lib.Primitive;
import nxm.sys.libm.Fft;

/**
  Perform Inverse FFT <br><br>
  <p style="text-align: center">
  **** ******************** *** <br>
  **    EXPERIMENTAL VERSION  *<br>
  **** ******************** *** <br>
  </p>
  @author Jeff Schoen
  @since NeXtMidas 3.3.2
*/
public class invfft extends Primitive {
   /** windowing function name list */
  private static final String windowList = DataOp.windowList;
  /** PARTIAL mode (process or discard partial input data frames) - NONE, ALL or START */
  private static final String partialList= "None,All,Start";
  /** PARTIAL mode: discard any partial input data insufficient to form a whole frame */
  private static final int P_NONE        = 1;
  /** PARTIAL mode: zero-pad any partial input data to produce an output frame */
  private static final int P_ALL         = 2;
  /** PARTIAL mode: zero-pad partial input data only at the start of processing;
   *  discard partial data if we have already output at least one frame */
  private static final int P_START       = 3;
  /** Packing mode (how the nyquist bin is handled) for real data (P_UNPACK ,P_DC or P_ZERO) */
  private static final String packList = "Unpack,DC,Zero";
  /** Packing mode that unpacks the nyquist into its own bin which may cause unoptimal transfer lengths */
  private static final int P_UNPACK    = 1;
  /** Packing mode that packs the real part of the nyquist bin into the always zero imaginary component of the P_DC bin */
  private static final int P_DC        = 2;
  /** Packing mode that drops the nyquist bin */
  private static final int P_ZERO      = 3;
  /** DataFile for inputs/outputs/cross output file */
  private DataFile hi1, ho1;
  /** Data buffers for inputs/outputs/computations */
  private Data dbi1,dbf1,dbw,dba1,db, dbr1, dbr2, dbo;
  private int nfft;      // transform size
  private int nin;       // transfer length
  private int nrot;      // number rotate for complex
  private int nout;      // number out
  private int stride;    // stride (num elements to skip)
  private int pack;      // packing mode
  private int partial;   // padding option for incomplete input
  private double over;   // fraction of {@link #nfft} to overlap blocks
  private int keepSamples;
  private double size;   // size of input files (smaller of IN1 or IN2)
  private boolean cxm;   // input data is in complex mode
  private boolean out1d; // output a type 1000 file
  private boolean outReal; // indicates a real output file
  private boolean in2d;  // input 1 is a type 2000 file
  private boolean dp;    // do double precision buffers & output
  private boolean deRotate; // deRotate input data to [0->Fs)
  private Fft     fft;   // FFT library
  private String  wname; // Name of windowing function

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

  public int open () {
    hi1 = MA.getDataFile ("IN","1000,2000","S#,C#",0);
    hi1.open(DataFile.WARNEMPTY);
    in2d = (hi1.typeClass==2);
    outReal = MA.getState( "/REAL", false );
    //  FFT explain file: default packing single frame=UNPACK, multi=ZERO
    if( in2d ) {
      pack = MA.getSelectionIndex("/PACK",packList, P_ZERO);
    }
    else {
      pack = MA.getSelectionIndex("/PACK",packList,P_UNPACK );
    }
    nfft  = MA.getL("NFFT");
    if (nfft == -1) nfft = (int)hi1.getSize(); //Matches X-Midas for -1 to get file size (2GB int limit)
    nin   = MA.getL("/TL",nfft);
    if( nin != nfft ) {
      throw new MidasException("Transfer lengths other than FFT size not yet implemented.");
    }
    if( outReal ) {
      //  Check the packing to determine the number of input samples
      if( pack == P_UNPACK ) {
        nin = nfft/2 + 1;
      } else {
        nin = nfft/2;
      }
    }
    out1d = MA.getState("/1D");
    partial = MA.getSelectionIndex("/PARTIAL",partialList,P_START);
    if (state!=RESTART) {                             // retain realtime changes
      over = MA.getD("OVER");
      wname = MA.getS("WIN");
    }
    double delf = hi1.getXDelta();
    double delx = 0;
    double dely;
    double stry; // input's Start (X Start for 1D, Y Start for 2D)
    cxm = (hi1.spa==2);
    if (in2d) {
      nin  = hi1.getSubSize();
      dely = hi1.getYDelta();
      stry = hi1.getYStart();
      //  The forward FFT sets the time to the center of thedata used to compute the FFT, adjust
      stry -= ( hi1.getYDelta() / 2.0 );
      hi1.setFS(0);
    } else {
      dely = hi1.getYDelta();
      stry = hi1.getYStart();
      dely = delx*stride;
      stry = hi1.getXStart();
    }
    hi1.setXferLength(nin);
    hi1.setConsLength(nin);
    setOverlap( over );
    size = hi1.size;                                  // estimate size
    int flags = Fft.INVERSE;
    byte dtype = (dp)? DOUBLE:FLOAT;
    String cform = (dp)? "CD":"CF";
    String sform = (dp)? "SD":"SF";
    deRotate = MA.getState("/ROTATE", true );
    if ( !outReal ) {
      nout = nfft;
      nrot = nfft/2-1;
      flags |= Fft.COMPLEX;
      pack = 0;                      // no packing for complex data
      if (!deRotate) nrot=0;
    } else {
      nout = nfft;
      nrot = 0;
      flags |= Fft.COMPLEX;
    }
    fft = new Fft(nfft,flags);
    String oform =   ( outReal )? sform : cform;  // output file format
    double osize = (out1d)? 1.0   : Math.max(1.0,Math.floor((size-nfft))); // output file size
    ho1 = MA.getDataFile("OUT","2000",oform,0);
    ho1.propagate(hi1,"TIMECODE|KEYWORDS"); // propagate timecode and keywords from IN1 to OUT1 (Bug 2101)
    ho1.setSubSize (keepSamples);
    ho1.setXStart(0.0);
    ho1.setXDelta(1.0/( delf * nfft ));
    ho1.setXUnits(Units.TIME);
    ho1.setYStart( stry );
    ho1.setYDelta( dely );
    if (!out1d) {
      ho1.setYUnits(Units.TIME);  // if not /1D then set YU=time (since 3.1.0 for BUG 2495)
    }
    ho1.setSize(osize);
    ho1.open(DataFile.WARNEMPTY|DataFile.OPTIONAL);
    setWindow (wname);
    dbi1 = hi1.getDataBuffer(nin,dtype);
    dbf1 = new Data (cform,nout);
    if( over > 0 ) {
      dbo = new Data( cform, nout );
    }
    dbr1 = new Data (cform,nout);
    dbr2 = new Data (sform,nout);
    dba1 = ho1.getDataBuffer(1,dtype);
    return NORMAL;
  }

  public int process () {
    int stat = NORMAL;
    int n1 = (hi1.isOpen)? hi1.read(dbi1) : 0;
    //  For outReal, the amount of data read will not match the requested FFT size because only half the spectrum exists
    if ( n1<0 ) {                               // EOF
      stat = FINISH;
      return stat;
    }
    if( outReal ) {
      n1 = nfft;
    }
    if( outReal ) {
      if (( n1 == nin*2 || partial != P_NONE ) ) {
        if (hi1.isOpen) processInput (n1,ho1,dbi1,dbf1,dba1);
      }
    } else {
      if (( n1 == nin || partial != P_NONE ) ) {
        if (hi1.isOpen) processInput (n1,ho1,dbi1,dbf1,dba1);
      }
    }
    if (out1d) stat = FINISH;   // our work here is done
    if (partial == P_START) partial = P_NONE; // don't want any more partials
    return stat;
  }

  private void processInput (int n, DataFile ho, Data dbi, Data dbf, Data dba) {
    try {
    dbf.setSPA(dbi.spa);                              // see dbf as input
    DataOp.copy( dbi, dbf );
    if( outReal ) {
      if( pack == P_DC ) {
        DataOp.copy( dbi, dbi.bpe, dbi.bpe, dbf, ( nfft - 1)*dbf.bpe, -dbf.bpe, dbf.bpe, nin-1 );
        for( int i= (nfft/2); i < nfft; i++ ) {
          dbf.setF(2*i+1, -dbf.getF(2*i+1));
        }
        dbf.setF(nfft, dbf.getF(1));
        dbf.setF( 1, (float)0.0 );
        //  Nyquist is packed into the imaginary part of (always real) DC
      } else if ( pack == P_ZERO ) { //  P_ZERO = no nyquist bin
        DataOp.copy( dbi, dbi.bpe, dbi.bpe, dbf, ( nfft - 1)*dbf.bpe, -dbf.bpe, dbf.bpe, nin-1 );
        for( int i= (nfft/2); i < nfft; i++ ) {
          dbf.setF(2*i+1, -dbf.getF(2*i+1));
        }
      } else { // P_UNPACK
        //  FFT on real input outputs only the positive half of the spectrum
        //  In "UNPACK" form, spectrum runs from dc to nyquist mirror/phase shift the positive spectrum
        DataOp.copy( dbi, dbi.bpe, dbi.bpe, dbf, ( nfft - 1)*dbf.bpe, -dbf.bpe, dbf.bpe, nin-1 );
        for( int i= (nfft/2) + 1; i < nfft; i++ ) {
          dbf.setF(2*i+1, -dbf.getF(2*i+1));
        }
      }
    }
    if (n<nfft) DataOp.zero(dbf,n,nfft-n);            // zero fill
    if (cxm&&dbi.spa==1) DataOp.cmplx(dbf,dbf,nfft,DataOp.SC);  // force complex FFT
    //  Rotate data back to 0->Fs if it has been rotated
    if( nrot != 0 ) {
      fft.unrotate( dbf, nrot );
    }
    fft.work (dbf);
    if( !( wname.equalsIgnoreCase( "none" )) ) {
      DataOp.copy( dbf, dbr1 );
      DataOp.div(dbr1,dbw,dbf);
    }
    dbf.setSPA(2); db=dbf;                            // see dbf as output
    //  Discard overlapped samples
    if( keepSamples < nin ) {
      //  Using overlap in the forward FFT will result in duplicate samples
      //  Discard the duplicates - not sure how this will work with pipes,
      //  but Jimbo Bauerbo said to undo as many of the forward options as possible.
      //  For a file, could just always write the entire fft output
      //  then backup by the number of samples to discard.  This would allow
      //  for an exact undoing - however, pipes (of this kind) cannot be
      //  backed up...
      //  So, just always keep the keepSamples.
      //  This will result in a loss of data on the final FFT
      DataOp.copy( dbf, dbo );
      dbo.setSize( keepSamples );
      db = dbo;
    }
    if( outReal ) {
      //  Keep the real portion of the buffer
      DataOp.real(db, dbr2);
      db=dbr2;
    }
    ho.write (db);             // write output
    }
    catch ( Exception e ) {
      e.printStackTrace();
    }
  }

  public int close () {
    fft.free();
    if (out1d && ho1.isOpen) ho1.setFS(0);
    hi1.close();
    ho1.close();
    return NORMAL;
  }

  public int processReady () {
    if (state!=PROCESS) return 0;
    int n = hi1.processReady(nin);
    if (ho1.isOpen) n = Math.min(n, ho1.processReady(1) );
    return n;
  }

  /** set the fractional overlap, zero for no overlap, negative to skip.
      @param value fractional overlap
   */
  public void setOverlap (double value) {
    over = value;
    keepSamples = Math.max(0,(int)((1.0-over)*nfft));
  }

  /** get the FFT size.
      @return FFT size
   */
  public int getNFft () { return nfft; }

  /** set the FFT size.
      @param value new FFT size
   */
  public void setNFft (int value) {
    if (value==nfft) return;
    MA.put("NFFT",""+value);
    setState(RESTART);
  }

  /** set the windowing function.
      @param value new windowing fuction to use
   */
  public void setWindow (String value) {
    // TODO adjust this when implementing zero padding for INVFFT
    // int nwind = Math.min(nin,nfft);
    int nwind = nfft;
    wname = value;
    dbw = DataOp.getWindow(wname,nwind,1.0,M);
    // Remove any zeros (doing a divide by)
    for( int i=0; i < dbw.getSize();i=i+1 ) {
      if( Math.abs(dbw.getF(i)) < 1.0e-4 ) {
        dbw.setF(i, (float)1.0 );
      }
    }
  }

}
