package nxm.dsp.prim;

import static java.lang.Math.cos;
import static java.lang.Math.sin;
import static nxm.sys.libm.Tolerance.equalsWithTolerance;

import nxm.dsp.lib.FilterFIR;
import nxm.sys.inc.InternalUseOnly;
import nxm.sys.inc.Units;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.DataOp;
import nxm.sys.lib.Keywords;
import nxm.sys.lib.MidasException;
import nxm.sys.lib.Primitive;
import nxm.sys.lib.Time;
import nxm.sys.libm.Fft;

/** Tune and/or Filter and/or Decimate a signal.
    @author  Ken Smith / Jason Broussard / J. Czechowski IV
 */
public class tuner extends Primitive {
  //BUG 1757  Tuner - IndexOutOfBoundsException from certain input file sizes
  /** PARTIAL mode (process or discard partial input data frames) - NONE, ALL or START. */
  public static final String partialList= "None,All,Start";
  /** PARTIAL mode: discard any partial input data insufficient to form a whole frame. */
  public static final int P_NONE        = 1;
  /** PARTIAL mode: zero-pad any partial input data to produce an output frame. */
  public 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.
   */
  public static final int P_START       = 3;

  /** Keyword update mode:
      <pre>
        AUTO     - Same as UPDATE if /SBT= is set on command
                   line or NONE if /SBT= is not set
        INSERT   - Updates SBT keyword if present, inserts it if
                   not present (insertion assumes INPUT at RF)
        NONE     - Do not update/insert, but warn if present
        IGNORE   - Ignore all keywords and do not issue any
                   warnings related to them
        UPDATE   - Updates SBT keyword if present
      </pre>
      @since NeXtMidas 3.3.1
   */
  @InternalUseOnly // Only intended for use within DSP option tree
  public static enum KwdUpdate {
    /** Same as UPDATE if /SBT= is set on command line or NONE if /SBT= is not set */
    AUTO,
    /** Updates SBT keyword if present, inserts it if not present (insertion assumes INPUT at RF) */
    INSERT,
    /** Do not update/insert, but warn if present*/
    NONE,
    /** Ignore all keywords and do not issue any warnings related to them */
    IGNORE,
    /** Updates SBT keyword if present */
    UPDATE
  };

  private DataFile hin;                // Input file
  private DataFile hout;               // Output file
  private DataFile hdebug;             // Debug file (output)
  private Data     buf;                // Data buffer
  private Data     filter;             // Data buffer containing the filter
  private String   outFormat;          // Output format
  private String   windowName;         // Window name (always "HANN" prior to NXM 3.3.1)
  private Double   inputIF;            // Input IF  (Hz)
  private Double   outputIF;           // Output IF (Hz)
  private Double   bandwidth;          // Filter BW (Hz) (null if n/a)
  private String   resOutputIF;        // Label for actual Output IF result (since NXM 3.3.1)
  private double   timeDelta;          // Time delta imposed by filter/dec  (since NXM 3.3.1)
  private boolean  setTimeDelta;       // Update based on timeDelta?        (since NXM 3.3.1)
  private double   phasor_real;        // Current phasor (real)
  private double   phasor_imag;        // Current phasor (imag)
  private double   dphasor_real;       // Phasor delta (real)
  private double   dphasor_imag;       // Phasor delta (imag)
  private int      nfft;               // FFT size
  private int      decimation;         // Decimation facrot (1=NoDecimation)
  private int      fftFlagsFwd;        // FFT flags (forward)
  private int      fftFlagsInv;        // FFT flags (inverse)
  private int      nbad;               // Number of filter transient points to throw away
  private int      partial;            // Padding option for incomplete input (P_NONE, P_ALL, P_START)
  private boolean  doTune;             // Do tuning (false if inputIF ~= outputIF)
  private boolean  doFilter;           // Use a filter?
  private boolean  extFilter;          // Use external filter file?
  private boolean  filterOut;          // Write filter out?
  private boolean  complexIn;          // Input data format
  private boolean  complexOut;         // Output data format
  private boolean  complex;            // Internal data format

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

  @Override
  public int open () {
    // Get command line arguments first (some are used in setting up the files)
    nfft           = MA.getL("/FFT",        8192);   //  8192  matches X-Midas SBTUNER
    windowName     = MA.getS("/WIND",       "HANN"); // "HANN" matches X-Midas SBTUNER
    inputIF        = MA.getD("CFREQ",       MA.getD("INPUTIF",   null));
    bandwidth      = MA.getD("BW",          MA.getD("BANDWIDTH", null));
    decimation     = MA.getL("DECIMATION",  1);
    resOutputIF    = MA.getU("ACTOUTPUTIF", null);
    setTimeDelta   = MA.getZ("/SETDELAY",   true);
    hin            = (MA.isPresent("INFILE"))? MA.getDataFile("INFILE")
                                             : MA.getDataFile("INPUT");
    partial        = MA.getSelectionIndex("/PARTIAL",partialList,P_NONE);
    buf            = new Data("CF", nfft*2);
    outputIF       = null;  // default value (set below as needed)
    extFilter      = false; // default value (set below as needed)
    filterOut      = false; // default value (set below as needed)
    doFilter       = false; // default value (set below as needed)
    complex        = false; // default value (set below as needed)
    doTune         = false; // default value (set below as needed)
    timeDelta      = 0.0;   // default value (set below as needed)
    nbad           = 0;     // default value (set below as needed)
    fftFlagsFwd    = Fft.FORWARD|Fft.UNPACKED;  // Must add COMPLEX or REAL to this!
    fftFlagsInv    = Fft.INVERSE|Fft.COMPLEX;   // Always COMPLEX for inverse

    if ((resOutputIF != null) && resOutputIF.isEmpty()) resOutputIF = null;

    DataFile  hfilter   = null;            // Filter file (IMPRESP or FILTOUT)
    double    tolr      = 1e-8;
    String    kwdSBT    = MA.getS("/SBT", "");
    KwdUpdate kwdUpdate = MA.getSelection("/KWDUPDATE", KwdUpdate.AUTO);
    int       ntaps     = MA.getL("/TAPS", MA.getL("/NTAP", nfft/8));

    if ((this instanceof sbtuner) && MA.isPresent("BWIMPRESP") && (bandwidth == null)) {
      // X-Midas SBTUNER allows <BWIMPRESP> to be <BW> or <IMPRESP>
      try { bandwidth = MA.getD("BWIMPRESP"); }
      catch (MidasException e) { hfilter = MA.getDataFile("BWIMPRESP"); }
    }

    if (MA.isPresent("IMPRESP")) {
      if (bandwidth != null)       M.warning("Don't specify both <BANDWIDTH> (or <BW>) and <IMPRESP>, <IMPRESP> overrides");
      if (hfilter   != null)       M.warning("Don't specify both <BANDWIDTH> (or <BW>) and <IMPRESP>, <IMPRESP> overrides");
      if (MA.isPresent("FILTOUT")) M.warning("Don't specify both <IMPRESP> and <FILTOUT>, <IMPRESP> overrides");
      hfilter = MA.getDataFile("IMPRESP");
    }

    // ==== CHECK WARNING/ERROR CONDITIONS =====================================
    if (MA.isPresent("INFILE")) {
      M.deprecate("Since 2.9.1, use of 'INFILE=' is deprecated. Please use 'INPUT='");
    }
    if (MA.isPresent("OUTFILE")) {
      M.deprecate("Since 2.9.1, use of 'OUTFILE=' is deprecated. Please use 'OUTPUT='");
    }
    if (MA.isPresent("IMP_RESPONSE")) {
      throw new MidasException("Use of IMP_RESPONSE= was deprecated in NeXtMidas 2.9.1 "
                             + "and has been removed, use IMPRESP=");
    }
    if (MA.isPresent("/PH")) {
      // We have not implemented this feature from X-Midas SBTUNER, but we want
      // to make sure to error out in case someone is porting code from X-Midas
      // and has it set.
      throw new MidasException("Use of /PH= is not currently supported.");
    }
    if (!setTimeDelta) {
      // Enable backwards-compatible behavior which is known to be incorrect.
      M.deprecate("Since 3.3.1, use of '/SETDELAY=FALSE' is deprecated.");
    }

    // ==== OPEN INPUT FILE ====================================================
    hin.open();
    if (hin.spa == 2)  {
      fftFlagsFwd |= Fft.COMPLEX;
      complex   = true; // internal data must be complex when input is complex
      complexIn = true;
    }
    else {
      fftFlagsFwd |= Fft.REAL;
      complexIn = false;
    }

    // ==== INIT FILTER (1 of 2) ===============================================
    if ((bandwidth != null) || (hfilter != null)) {
      if ((hfilter != null)) {
        extFilter = true;
        hfilter.open();
      }
      doFilter = true;
      complex  = true;  // The internal data is always complex when a filter is used

      // Now see if we're supposed to write the impulse response out
      if (!extFilter && MA.isPresent("FILTOUT")) {
        hfilter = MA.getDataFile("FILTOUT", hin, DataFile.OUTPUT);
        hfilter.setFormat("CF");
        hfilter.open();
        filterOut = hfilter.isOpen();
      }
    }

    // ==== OPEN OUTPUT FILE ===================================================
    outFormat = MA.getS("OUTFORMAT", hin.getFormat());
    // TODO: comment out next 2 lines as X-Midas now supports D
    if (outFormat.equals("SD")) outFormat = "SF"; // <-- Matches X-Midas SBTUNER
    if (outFormat.equals("CD")) outFormat = "CF"; // <-- Matches X-Midas SBTUNER

    if (MA.isPresent("OUTPUT") || MA.isPresent("OUTFILE")) {
      hout = (MA.isPresent("OUTFILE"))? MA.getDataFile("OUTFILE", hin, DataFile.OUTPUT)
                                      : MA.getDataFile("OUTPUT",  hin, DataFile.OUTPUT);
      hout.setFormat(outFormat);
      hout.setXDelta(hin.getXDelta()*decimation);
      hout.open();
      complexOut = (hout.spa == 2);
    }
    else {
      hout       = null;
      complexOut = false;
    }

    // ==== DETERMINE PROCESSING REQUIRED ======================================
    if (inputIF != null) {
      String outIF = MA.getS("OUTPUTIF", "");

      outputIF = (outIF.isEmpty())? inputIF : MA.getD("OUTPUTIF");
      if (equalsWithTolerance(inputIF, outputIF, tolr)) {
        doTune   = false;
        outputIF = inputIF;
      }
      else {
        // An output IF of -1 implies center it in the output band,
        // either 0 for complex or fs/4 for real
        if (equalsWithTolerance(outputIF, -1.0, tolr)) {
          outputIF = (complexOut)? 0.0 : 0.25 / hin.getXDelta(); // NXM-4055 fix OUTPUTIF for real out
        }
        double dphase = (outputIF - inputIF) * TWOPI * hin.getXDelta()*decimation;
        dphasor_real = cos(dphase);
        dphasor_imag = sin(dphase);
        phasor_real  = 1.0;
        phasor_imag  = 0.0;
        doTune       = true;
      }
    }

    // ==== INIT FILTER (2 of 2) ===============================================
    if (doFilter) {
      if (extFilter) {
        filter    = hfilter.getDataBuffer(nfft*2, FLOAT);
        ntaps     = hfilter.read(filter);
        timeDelta = -(hfilter.getXStart() + hin.getXDelta())
                   * (hin.getXDelta() / hfilter.getXDelta());
      }
      else {
      // TODO: Add support for BANDWIDTH=0 to match X-Midas
//        if (bandwidth<=0) { // X-Midas added this
//          bandwidth = (complexIn) ? 2 * hin.getNyquist() : hin.getNyquist();
//        }
        FilterFIR f  = FilterFIR.initBandPass(M, ntaps, inputIF, bandwidth, hin.getNyquist(),
                                              decimation, complexOut, true, true, windowName, 1.0);
        filter    = new Data(f.makeFilter(nfft*2), 2);
        timeDelta = f.getDelay(hin.getXDelta());

        // Write the filter out, if requested
        if (filterOut) {
          if (setTimeDelta) {
            hfilter.setXStart(-timeDelta-hin.getXDelta()); // <-- since NeXtMidas 3.3.1
          }
          hfilter.setXDelta(hin.getXDelta());
          hfilter.setXUnits(Units.TIME_S); // <-- since NeXtMidas 3.3.1
          hfilter.write(filter, ntaps);
        }
      }

      // Zero-pad, scale and FFT the impulse response (remember it's complex)
      // Note: We do this here rather than in the process() as is done by the
      //       X-Midas SBTUNER since this saves us a call to DataOp.mul(..) with
      //       the scale factor. If we omit the scale factor the output result
      //       will be way, way bigger than it should be.
      DataOp.mul(filter, new Data(1.0/nfft), filter, ntaps, DataOp.CSC);
      DataOp.fill(filter, ntaps, 0.0, nfft-ntaps);
      DataOp.fft(filter, nfft, Fft.FORWARD|Fft.COMPLEX);

      nbad    = ntaps - 1;  // Calculate number of filter transient points to throw away
      complex = true;       // The internal data is always complex when a filter is used
    }

    // ==== HANDLE DELAY THROUGH FILTER (since NeXtMidas 3.3.1) ================
    if (setTimeDelta && (hout != null)) {
      hout.setXStart(hout.getXStart() + timeDelta);
    }

    // ==== UPDATE KEYWORDS (since NeXtMidas 3.3.1) ============================
    // The computation for newSBT (when oldSBT!=null) is taken directly from
    // SBTUNER in X-Midas.
    if ((hout != null) && doTune && (kwdUpdate != KwdUpdate.IGNORE)) {
      if (kwdUpdate == KwdUpdate.AUTO) {
        kwdUpdate = (kwdSBT.isEmpty())? KwdUpdate.NONE : KwdUpdate.UPDATE;
      }
      if (kwdSBT.isEmpty()) {
        kwdSBT = "SBT"; // This must come AFTER the check for KwdUpdate.AUTO
      }
      Keywords kwds    = hout.getKeywordsObject();
      Double   oldSBT  = kwds.getD(kwdSBT, null);
      boolean  hasSBT  = (oldSBT != null);
      double   fsIn    = 1 / hin.getXDelta();
      double   fsOut   = fsIn / decimation;
      double   fMinIn  = (complexIn )? -fsIn/2  : 0.0;
      double   fMinOut = (complexOut)? -fsOut/2 : 0.0;
      double   fbwIn   = (fsIn/2)  - fMinIn;
      double   fbwOut  = (fsOut/2) - fMinOut;
      double   ftune   = outputIF - inputIF;
      double   newSBT  = (oldSBT != null)? (oldSBT + (fMinOut + fbwOut/2) - (fMinIn + fbwIn/2) - ftune)
                                         : (0.0    + (fMinOut + fbwOut/2) - (fMinIn + fbwIn/2) - ftune);

      switch (kwdUpdate) {
        case INSERT:             kwds.put(kwdSBT, new Data(newSBT)); break;
        case UPDATE: if (hasSBT) kwds.put(kwdSBT, new Data(newSBT)); break;
        case NONE:   if (hasSBT) M.warning("Skipping update for "+kwdSBT); break;
        default:     throw new AssertionError("Unexpected KwdUpdate."+kwdUpdate);
      }

      if (kwds.find("SBT_OFFSET_CORRECT") != null) {
        // We look for SBT_OFFSET_CORRECT since it is an obscure keyword to find
        // outside of a Platinum file with SBT keywords.
        M.warning("Skipping update for (Platinum) SBT_* keywords.");
      }
    }

    // ==== SET XFER/CONS LENGTH (must come after filter init) =================
    hin.xferLength = nfft;
    hin.consLength = ((hin.xferLength - nbad) / decimation) * decimation;

    // ==== FINAL ERROR CHECKING ===============================================
    if (doTune && !doFilter && !complex) {
      // We can't tune a real input unless we filter it (since the Hilbert
      // transform is part of the filter).
      throw new MidasException("You must filter a real input to tune it!");
    }

    // ==== SET RESULTS & CLOSE FILTER FILE ====================================
    if (resOutputIF != null) {
      // Since 3.3.1 write actual outputIF to result (matches X-Midas SBTUNER)
      MR.put(resOutputIF, outputIF);
    }

    if (hfilter != null) {
      hfilter.close();
    }

    return (hout != null)? NORMAL : FINISH;
  }

  @Override
  public int process () {
    // ==== READ INPUT FILE ====================================================
    int ngot      = hin.read(buf);
    int ngot_orig = ngot;   //BUG 1757 - Need the orig. ngot because ngot gets
                            //           recalculated by subtracting nbad from it
    if (ngot == 0) return NOOP;
    if (ngot <  0) return FINISH;

    // ==== FILTER DATA (if necessary) =========================================
    if (doFilter) {
      // Zero-pad out to the FFT size
      DataOp.fill(buf, ngot, 0.0, nfft-ngot);
      DataOp.fft(buf, nfft, fftFlagsFwd);

      // Conjugate reflection to create negative frequency part of complex FFT,
      //  as done in X-Midas primitive SBTUNER.
      if (!complexIn) reflect(buf);

      DataOp.mul(buf, filter, buf, nfft, DataOp.CCC);
      DataOp.fft(buf, nfft, fftFlagsInv);
    }

    // ==== TUNE & DECIMATE DATA (as necessary) ================================
    // We need to trim off the "bad" points at the start of the array, then
    // decimate and/or tune. This isn't easily vectorizable, so do it on the array.
    float[] fbuf    = buf.castF(true);
    int     nout    = (ngot - nbad)/decimation;
    int     nscalar = (complexOut)? nout*2 : nout;

    if (doTune) {                 // doTune -> complex=true
      for (int i=nbad*2, j=0; j < nscalar; ) {
        fbuf[j++] = (float)(fbuf[i] * phasor_real - fbuf[i+1] * phasor_imag);
        if (complexOut) {
          fbuf[j++] = (float)(fbuf[i+1] * phasor_real + fbuf[i] * phasor_imag);
        }
        double rtemp = phasor_real*dphasor_real - phasor_imag*dphasor_imag;
        double itemp = phasor_imag*dphasor_real + phasor_real*dphasor_imag;
        phasor_real = rtemp;
        phasor_imag = itemp;
        i += 2;

        // Now perform the decimation
        if (complex) i += 2*(decimation-1);
        else         i += decimation-1;
      }
//double phasorAbs = Math.sqrt(phasor_real*phasor_real + phasor_imag*phasor_imag);
//phasor_real = phasor_real / phasorAbs;
//phasor_imag = phasor_imag / phasorAbs;
    }
    else {
      for (int i=nbad*2, j=0; j < nscalar; ) {
        fbuf[j++] = fbuf[i++];
        if (complexOut) {
          fbuf[j++] = fbuf[i];
        }
        // Skip the imaginary part (if it exists) regardless of
        // whether we copied it over above
        if (complex) i++;

        // Now perform the decimation
        if (complex) i += 2*(decimation-1);
        else         i += decimation-1;
      }
    }
    buf.uncast(fbuf, true);
    ngot -= nbad;

    // ==== PAD OUTPUT (as necessary, see Bug 1757) ============================
    // BUG 1757  Tuner - IndexOutOfBoundsException from certain input file sizes
    //  If the amount of data to output is < hin.consLength, then we have to pad the
    //  output data buffer with 0's if the user specified as such with the /partial switch
    if (nout < hin.consLength) {
      if (partial == P_ALL || partial == P_START) {
        DataOp.fill(buf, ngot_orig, 0.0, buf.size - ngot_orig);
        // populated the frame out to consLength, so we can reset hout to consLength
        // to output a full frame of data (padded out)
        nout = hin.consLength;
      }
    }

    // BUG 1757 - P_START only applies to the first frame of data.  If the first
    // frame is not padded, then P_START should be set to P_NONE, because
    // /partial=start == /partial=none after the first frame of data is processed.
    if (partial == P_START) {
      partial = P_NONE;
    }

    // ==== ACCOUNT FOR FILTER DELAY (see Bug 2698) & OUTPUT DATA ==============
    if (nout > 0) {
      // Update time stamp in output header
      Object hdr = buf.getHeader();
      if (!setTimeDelta) {
        buf.setHeader(null); // <-- Matches NeXtMidas 3.2 (though not very useful)
      }
      else if (hdr instanceof Time) {
        buf.setHeader(((Time)hdr).addSec(timeDelta));
      }

      // Write out data
      hout.write(buf, nout);
    }
    return NORMAL;
  }

  @Override
  public int close () {
    hin.close();
    if (hout    != null) hout.close();
    if (hdebug  != null) hdebug.close();

    if (resOutputIF != null) {
      // Since 3.3.1 write actual outputIF to result (matches X-Midas SBTUNER)
      MR.put(resOutputIF, outputIF);
    }
    return NORMAL;
  }

  private void reflect (Data buffer) {
    float[] fbuf = buffer.castF(true);

    // We're looking a a complex array as a real one, so
    //  our indices are doubled
    for (int i = 1; i < nfft/2; i++) {
      int j1 = (nfft - i)*2;
      int j2 = i*2;
      fbuf[j1++] = fbuf[j2++];  // Real part
      fbuf[j1] = - fbuf[j2];    // Imag part
    }
    buffer.uncast(fbuf,true);
  }
}
