package nxm.dsp.prim;

import static java.lang.Double.NEGATIVE_INFINITY;
import static java.lang.Double.POSITIVE_INFINITY;
import static java.lang.Math.abs;
import static java.lang.Math.ceil;
import static java.lang.Math.floor;
import static java.lang.Math.log10;
import static java.lang.Math.max;
import static java.lang.Math.min;
import static java.lang.Math.pow;
import static java.lang.Math.rint;
import static java.lang.Math.sqrt;
import static nxm.sys.lib.StringUtil.padLeft;

import java.util.LinkedList;
import java.util.List;
import java.util.Objects;

import nxm.sys.inc.InternalUseOnly;
import nxm.sys.inc.Tablizable;
import nxm.sys.inc.Units;
import nxm.sys.lib.Convert;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.DataOp;
import nxm.sys.lib.Midas;
import nxm.sys.lib.MidasException;
import nxm.sys.lib.Primitive;
import nxm.sys.lib.Table;
import nxm.sys.lib.Time;
import nxm.sys.libm.Fft;

/** Computes statistics about a signal.
 *  @author J. Czechowski IV (signal statistics version)
 *  @author Jeff Schoen      (initial FFT version)
 *  @since NeXtMidas 3.1.2
 */
public class sigstatistics extends Primitive  {
  // Implementation Notes:
  //   1. This currently assumes that the output from the FFT is REAL (i.e. not COMPLEX). There
  //      are some portions of the primitive that will work with COMPLEX output (generally
  //      portions copied from STATISTICS); however, others will clearly need updating if COMPLEX
  //      output from the FFT is ever used (which is highly unlikely since the PSD algorithm only
  //      works with a REAL output).
  //   2. Prior to NeXtMidas 4.1.4 this did not support:
  //        a. Complex input
  //        b. Sub-band tuned files (this was unintentional and was due minor issues
  //           with the way some of the keywords were handled)
  //        c. Modern keywords (i.e. JICD 4.2)
  //        d. Extra output statistics
  //   3. The version prior to NeXtMidas 4.1.4 has a number of documentation errors, particularly
  //      related to the way the algorithm was documented. This was the result of an oversight when
  //      this primitive was first imported into NeXtMidas (by those unfamiliar with the code) and
  //      has been subsequently corrected.
  //   4. A number of the computations in the pre-NeXtMidas 4.1.4 version used float rather than
  //      double values as intermediates. This was changed to using double values since, although
  //      the extra precision is not needed, keeping all of the computations using double values
  //      is more efficient (given the FPU on the x86/x86_64 chipset).

  /** windowing function name list */
  public static final String windoList     = nxm.sys.prim.fft.windowList;

  /** PARTIAL mode (process or discard partial input data frames) - NONE, ALL or START */
  public static final String partialList   = nxm.sys.prim.fft.partialList;

  /** PARTIAL mode: discard any partial input data insufficient to form a whole frame */
  public static final int P_NONE           = nxm.sys.prim.fft.P_NONE;

  /** PARTIAL mode: zero-pad any partial input data to produce an output frame */
  public static final int P_ALL            = nxm.sys.prim.fft.P_ALL;

  /** 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          = nxm.sys.prim.fft.P_START;

  /** Packing mode (how the nyquist bin is handled) for real data (P_UNPACK, P_DC or
   *  P_ZERO)
   */
  public static final String packList      = nxm.sys.prim.fft.packList;

  /** Packing mode that unpacks the nyquist into its own bin which may cause
   *  sub-optimal transfer lengths
   */
  public static final int P_UNPACK         = nxm.sys.prim.fft.P_UNPACK;

  /** Packing mode that packs the real part of the nyquist bin into the always zero
   *  imaginary component of the P_DC bin
   */
  public static final int P_DC             = nxm.sys.prim.fft.P_DC;

  /** Packing mode that drops the nyquist bin */
  public static final int P_ZERO           = nxm.sys.prim.fft.P_ZERO;

  /** Indicator that there is no applicable result. */
  private static final SNRResult NO_RESULT = new SNRResult();

  private       DataFile inFile;        // Input file     (required)
  private       DataFile outFile;       // Output file    (optional)
  private       DataFile fftFile;       // Output file    (optional) [FFT output for debugging]
  private       String   outResult;     // Output result  (optional)
  private       Data     inData;        //
  private       Data     fftData;       //
  private       Data     outData;       //
  private       Data     winData;       //
  private       Data     addData;       //
  private       Data     acData;        //
  private       int      nfft;          // transform size
  private       int      nin;           // transfer length
  private       int      nrot;          // number rotate for complex
  private       int      nout;          // number out
  private       int      navg;          // num transforms to average
  private       int      iavg;          // averaging counter
  private       int      stride;        // stride (num elements to skip)
  private       int      nexp;          // num exponential averages
  private       int      pack;          // packing mode
  private       int      partial;       // padding option for incomplete input
  private       double   factor;        // windowing scale factor
  private       double   over;          // fraction of {@link #nfft} to overlap blocks
  private       double   size;          // size of input files (smaller of IN1 or IN2)
  private       double   rnavg;         // 1/navg or 1/navg/nexp 1/{@link #navg} or 1/{@link/navg}/{@link #nexp}
  private       double   rnexp;         // ~(1-1/nexp)/rnavg (1-1/@link #nexp})*{@link #nexp}*{@link #navg}
  private       boolean  mag;           // output magnitude squared (MAG)
  private final boolean  psd   = true;  // output power spectral density (PSD) [matches fft.java -- required]
  private final boolean  log   = false; // output MAG or PSD in 20log format [matches fft.java -- not used]
  private final boolean  fdb   = true;  // use fast dB log computation [matches fft.java -- required]
  private       boolean  cxm;           // input data is in complex mode
  private final boolean  out1d = false; // output a type 1000 file [matches fft.java -- required]
  private       boolean  in2d;          // input 1 is a type 2000 file
  private       boolean  dp;            // do double precision buffers & output
  private       boolean  ac;            // Take out DC bias before window and FFT
  private       String   wname;         // Name of windowing function
  private       Fft      fft;           // FFT library
  private       double   fftFreqMin;    // FFT output minimum frequency
  private       double   fftFreqMax;    // FFT output maximum frequency
  private       double   fftFreqCentre; // FFT centre frequency
  private       double   fftFreqDelta;  // FFT output frequency delta
  private       double   fftBandwidth;  // FFT output bandwidth
  private       double   signalFreq;    // Frequency of the signal of interest (-1 if n/a)
  private       double   signalBW;      // Bandwidth of the signal of interest  (-1 if n/a)
  private       double   rcvrCF;        // Receiver center freq
  private       int      rcvrTF;        // Receiver translation freq
  private       double   rcvrBW;        // Receiver bandwidth
  private       Double   startTime;     // The start time in J1950 seconds (null if n/a).
  private       Double   endTime;       // The end   time in J1950 seconds (null if n/a).
  private       String   kwMode;        // keyword mode
  private       int      signalBinMin;  // First FFT bin for the signal of interest (-1 if n/a)
  private       int      signalBinMax;  // Last FFT bin for the signal of interest (-1 if n/a)
  private       int      rcvrBinMin;    // First FFT bin for the receiver's bandwidth (0 if n/a)
  private       int      rcvrBinMax;    // Last FFT bin for the receiver's bandwidth (N if n/a)
  private       int      maxHistPoints; // Maximum number of historical points to average over
  private       Table    statTable;     // The most recent statistics

  private final List<SNRResult> histSignal  = new LinkedList<>(); // Historical signal values
  private final List<SNRResult> histPeak    = new LinkedList<>(); // Historical peak values

  //////////////////////////////////////////////////////////////////////////////
  // THIS GOES IN COMMAND & COMMANDABLE
  //////////////////////////////////////////////////////////////////////////////
  // in Commandable add:
  //    private static final int COMPLETED = 12;
  //
  // in Command add:
  //    private boolean waitForPipe;  // is the /WAITFORPIPE switch set?
  //
  // in Command.initSwitches() add:
  //    waitForPipe = MA.getState("/WAITFORPIPE");
  //
  //
  // in Command.runSingle() change:
  //    status = process();
  //    if (status!=NORMAL) break;
  // to
  //    status = process();
  //    if ((status==FINISH) && waitForPipe) status = COMPLETED;
  //    if (status!=NORMAL) break;
  //
  // The code immediately following, mimics this, but is limited to functions that can be overridden.
  private static final int COMPLETED = 12;
  private boolean waitForPipe;

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

  @Override
  public synchronized int checkin (Midas M) {
    int status  = super.checkin(M);
    waitForPipe = MA.getState("/WAITFORPIPE");
    return status;
  }

  @Override
  public int process () {
    int status = _process();
    if (waitForPipe && (status == FINISH)) status = COMPLETED;
    return status;
  }
  //////////////////////////////////////////////////////////////////////////////

  @Override
  @SuppressWarnings("fallthrough")
  public int open () {
    inFile        = MA.getDataFile("IN", "1000","S#,C#",0); inFile.open(DataFile.WARNEMPTY);
    in2d          = (inFile.typeClass == 2);
    nfft          = MA.getL("NFFT");                    if(nfft == -1) nfft = (int)inFile.getSize();
    nin           = MA.getL("/TL", nfft);
    mag           = fdb||log||psd;
    ac            = MA.getState("/AC");
    dp            = MA.getState("/DP");
    partial       = MA.getSelectionIndex("/PARTIAL", partialList, P_START);
    signalFreq    = MA.getD("FREQ", -1);
    signalBW      = MA.getD("BW", -1);
    signalBinMin  = -1; // set below
    signalBinMax  = -1; // set below
    maxHistPoints = MA.getL("/AVGPTS", 1024);
    statTable     = new Table();
    outResult     = MA.getU("/RES", null);
    rcvrCF        = MA.getD("/RCVR_CF");
    rcvrBW        = MA.getD("/RVCR_BW");
    kwMode        = MA.getS("/KWMODE", "AUTO");
    startTime     = MA.getD("STARTTIME", null);
    endTime       = MA.getD("ENDTIME",   null);

    if ((startTime != null) && (startTime <= 0)) startTime = null;
    if ((endTime   != null) && (endTime   <= 0)) endTime   = null;

    if (outResult != null) {
      MR.put(outResult, getStatTable());
    }

    if (state != RESTART) {                           // retain realtime changes
      navg  = MA.getL("/NAVG", 5);
      over  = MA.getD("/OVER", 0.25);
      nexp  = MA.getL("/NEXP", 1);
      wname = MA.getS("/WIN", "HANN");
    }

    double delx = inFile.getXDelta();
    double delf = 1.0/(nfft*delx);
    double dely;
    double stry; // input's Start (X Start for 1D, Y Start for 2D)
    cxm = (inFile.spa == 2);

    if (in2d) {
      nin  = inFile.getSubSize();
      dely = inFile.getYDelta();
      stry = inFile.getYStart();
      inFile.setFS(0);
      setOverlap(over);                                     // initialize stride
    } else  {
      setOverlap(over);                                     // initialize stride
      dely = delx*stride;
      stry = inFile.getXStart();
    }

    size = inFile.size;                                    // estimate size
    if (navg < 0) navg = (int)((size - nfft) / stride) +1; // average the entire file
    setNExp (nexp);                                        // needs navg
    setNAvg (navg);                                        // needs nexp

    int flags    = Fft.FORWARD;
    byte dtype   = (dp) ? DOUBLE : FLOAT;
    String cform = (dp) ? "CD" : "CF";
    String sform = (dp) ? "SD" : "SF";
    pack = out1d ? P_UNPACK : P_ZERO;       // unpack only if single output frame (i.e. /1D)
    pack = MA.getSelectionIndex("/PACK", packList, pack);

    if (cxm) {
      nout = nfft;
      nrot = nfft/2-1;
      flags |= Fft.COMPLEX;
      factor = 1.0/nfft;
      pack = 0;               // no packing for complex data
    }
    else {
      nout = nfft/2;
      if (nfft != nout*2) {
        // This can happen if nfft is odd
        throw new MidasException("Invalid FFT size parameter: "+nfft);
      }
      nrot = 0;
      flags |= Fft.REAL;
      if (pack == P_UNPACK) { flags |= Fft.UNPACKED; nout += 1; }
      else                  { flags |= Fft.PACKED;              }
      factor = 1.0/(2*nfft);
    }

    fft = new Fft(nfft, flags);
    if (psd) factor /= sqrt(delf);
    if (!MA.getState("/ROTATE", cxm)) nrot = 0;

    fftFreqDelta  = delf;
    fftFreqMin    = -nrot*delf;
    fftBandwidth  = delf * nout;
    fftFreqMax    = fftFreqMin + (fftBandwidth);
    fftFreqCentre = fftFreqMin + (fftBandwidth * 0.50);

    String oform = (mag) ? sform : cform;   //output file format
    double osize = (out1d) ? 1.0 : max(1.0, floor((size - nfft)/stride/navg));  //output file size
    fftFile = MA.getDataFile("FFT", "2000", oform, 0);
    fftFile.propagate(inFile, "TIMECODE|KEYWORDS");
    fftFile.setSubSize(nout);
    fftFile.setXStart(fftFreqMin);
    fftFile.setXDelta(delf);
    fftFile.setXUnits(Units.FREQUENCY);
    fftFile.setYStart(stry + (nfft + (navg - 1) * stride) / 2 * delx);
    fftFile.setYDelta(dely * navg);
    fftFile.setYUnits(Units.TIME);
    fftFile.setSize(osize);
    fftFile.open(DataFile.WARNEMPTY | DataFile.OPTIONAL);

    outFile = MA.getDataFile("OUT");
    outFile.propagate(inFile, "TIMECODE|KEYWORDS");
    outFile.setXStart(stry + (nfft + (navg - 1) * stride) /2 * delx);
    outFile.setXDelta(dely * navg);
    outFile.setYUnits(Units.TIME);
    Table quals = outFile.getQualifiers();
    int   type  = (quals == null) ? 6000 : quals.getL("TYPE", 6000);
    if ((type / 1000) < 6) {
      outFile.setType(type);
      outFile.addSubRec("NFLR", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("NTHR", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIG",  "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SSNR", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("FREQ", "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("PEAK", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("PSNR", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("ANF",  "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("ANT",  "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("ASIG", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("ASNR", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AFRQ", "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("APSG", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("APSN", "SL", -1, 0, Units.REL_POWER_DB);
    }
    else {
      outFile.setType(type);
      outFile.addSubRec("NOISE_FLOOR",              "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("NOISE_THRESHOLD",          "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIGNAL_STRENGTH",          "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIGNAL_SNR",               "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("PEAK_SIGNAL_FREQ",         "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("PEAK_SIGNAL_STRENGTH",     "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("PEAK_SIGNAL_SNR",          "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AVG_NOISE_FLOOR",          "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AVG_NOISE_THRESHOLD",      "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AVG_SIGNAL_STRENGTH",      "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AVG_SIGNAL_SNR",           "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AVG_PEAK_SIGNAL_FREQ",     "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("AVG_PEAK_SIGNAL_STRENGTH", "SL", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("AVG_PEAK_SIGNAL_SNR",      "SL", -1, 0, Units.REL_POWER_DB);

      // SIGNAL_* entries added in NeXtMidas 4.1.4...
      outFile.addSubRec("SIGNAL_MEAN",              "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIGNAL_SDEV",              "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIGNAL_ELEM",              "SL", -1, 0, Units.NONE        );
      outFile.addSubRec("SIGNAL_MIN",               "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIGNAL_MAX",               "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SIGNAL_AMIN",              "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("SIGNAL_AMAX",              "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("SIGNAL_OMIN",              "SL", -1, 0, Units.NONE        );
      outFile.addSubRec("SIGNAL_OMAX",              "SL", -1, 0, Units.NONE        );
      outFile.addSubRec("SIGNAL_SUM",               "SF", -1, 0, Units.REL_POWER_DB);

      // SPECTRUM_* entries added in NeXtMidas 4.1.4...
      outFile.addSubRec("SPECTRUM_MEAN",            "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SPECTRUM_SDEV",            "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SPECTRUM_ELEM",            "SL", -1, 0, Units.NONE        );
      outFile.addSubRec("SPECTRUM_MIN",             "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SPECTRUM_MAX",             "SF", -1, 0, Units.REL_POWER_DB);
      outFile.addSubRec("SPECTRUM_AMIN",            "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("SPECTRUM_AMAX",            "SL", -1, 0, Units.FREQUENCY_HZ);
      outFile.addSubRec("SPECTRUM_OMIN",            "SL", -1, 0, Units.NONE        );
      outFile.addSubRec("SPECTRUM_OMAX",            "SL", -1, 0, Units.NONE        );
      outFile.addSubRec("SPECTRUM_SUM",             "SF", -1, 0, Units.REL_POWER_DB);
    }
    outFile.setSize(osize);
    outFile.open(DataFile.WARNEMPTY | DataFile.OPTIONAL);

    setWindow(wname);
    inData  = inFile.getDataBuffer(nin, dtype);
    fftData = new Data(cform, nout);
    outData = fftFile.getDataBuffer(1, dtype);
    addData = fftFile.getDataBuffer(1, dtype);

    // LOOK UP RECEIVER CF AND BW KEYWORDS
    switch (kwMode.toUpperCase()) {
      case "AUTO":
        // This deprecation warning added in NeXtMidas 4.1.4
        M.deprecate("/KWMODE=AUTO currently matches /KWMODE=LEGACY to be backwards-compatible "
                   +"with NeXtMidas 4.1.3; this will change in a future release to give preference to the "
                   +"/KWMODE=JICD keywords");
        // FALLTHROUGH

      case "LEGACY":
        // This is the set of "commonly-used keyword names," that were supported back when
        // this primitive was first written. Considered "legacy" since NeXtMidas 4.1.4.
        if (rcvrCF <= 0) rcvrCF = Convert.o2d(inFile.getKeyWord("TLA_VRF"        )) * 1e6;
        if (rcvrCF <= 0) rcvrCF = Convert.o2d(inFile.getKeyWord("TLA_CENTER_FREQ")) * 1e6;
        if (rcvrCF <= 0) rcvrCF = Convert.o2d(inFile.getKeyWord("VRF"            )) * 1e6;
        if (rcvrCF <= 0) rcvrCF = Convert.o2d(inFile.getKeyWord("ETF"            )) * 1e6;

        if (rcvrBW <= 0) rcvrBW = Convert.o2d(inFile.getKeyWord("TLA_VRBW"       )) * 1e6;
        if (rcvrBW <= 0) rcvrBW = Convert.o2d(inFile.getKeyWord("TLA_BANDWIDTH"  )) * 1e6;
        if (rcvrBW <= 0) rcvrBW = Convert.o2d(inFile.getKeyWord("VRBW"           )) * 1e6;
        if (rcvrBW <= 0) rcvrBW = Convert.o2d(inFile.getKeyWord("EBW"            )) * 1e6;
        break;

      case "JICD":
      case "JICD4.2":
      case "JICD 4.2":
        // These are the keywords supported by JICD 4.2 and many of the more "modern" tools
        // which defer to base-units for all values. Take care to note the difference in units
        // for VRF between LEGACY and JICD 4.2 mode.
        if (rcvrCF <= 0) rcvrCF = Convert.o2d(inFile.getKeyWord("VRF"      ));
        if (rcvrBW <= 0) rcvrBW = Convert.o2d(inFile.getKeyWord("BANDWIDTH"));
        break;

      default:
        throw new MidasException("Unsupported keyword mode /KWMODE="+kwMode);
    }

    // SETUP THE IF->RF TRANSLATIONS & RECEIVER BW FILTERS
    if (MA.getState("/IF", false)) {
      rcvrCF = fftFreqCentre;
      rcvrTF = 0;
    }
    else if (rcvrCF <= 0) {
      M.warning("No receiver RVCR_CF given, assuming " + signalFreq + " is IF.");
      rcvrCF = fftFreqCentre;
      rcvrTF = 0;
    }
    else {
      rcvrTF = (int)rint(rcvrCF - fftFreqCentre);
      fftFile.setXStart(fftFreqMin + rcvrTF);
    }

    if (rcvrBW <= 0) {
      M.warning("No receiver RVCR_BW given, assuming full-band.");
      rcvrBW = fftBandwidth;
    }
    else if (rcvrBW > fftBandwidth) {
      M.warning("Receiver RVCR_BW exceeds full-band, assuming full-band.");
      rcvrBW = fftBandwidth;
    }

    if (rcvrBW <= 0) {
      rcvrBinMin = 0;
      rcvrBinMax = nout;
    }
    else {
      double minFreq = fftFreqCentre - rcvrBW/2;
      double maxFreq = fftFreqCentre + rcvrBW/2;
      rcvrBinMin = (int)floor((minFreq - fftFreqMin) / fftFreqDelta);  // inclusive
      rcvrBinMax = (int)ceil( (maxFreq - fftFreqMin) / fftFreqDelta);  // exclusive
    }

    // SETUP THE SIGNAL BW FILTERS
    double minFreq = (signalFreq - rcvrTF) - signalBW/2;  // signal min
    double maxFreq = (signalFreq - rcvrTF) + signalBW/2;  // signal max

    if ((minFreq >= fftFreqMin) && (maxFreq <= fftFreqMax)) {
      signalBinMin = (int)floor((minFreq - fftFreqMin) / fftFreqDelta);  // inclusive
      signalBinMax = (int)ceil( (maxFreq - fftFreqMin) / fftFreqDelta);  // exclusive

      if ((signalBW / fftFreqDelta) < 0.25) {
        // Anything sufficiently smaller than 1.00 bin will start having
        // problems, usually the real problems are seen around 0.10 and lower.
        // At that point the SNR will typically be under reported by 2-4 dB
        // for every doubling of the bin size (i.e. halving of the number of FFT
        // bins).
        //
        // The most common cases where a user will see this are:
        //  1. Choosing a NFFT that is just way too small
        //  2. The signal is (or approximates) a tone.
        M.warning("Signal bandwidth ("+signalBW+" Hz) is sufficiently below the width of the "
                 +"FFT bins ("+fftFreqDelta+" Hz) that the SNR results will be unreliable.");
      }
    }
    else {
      signalBinMin = -1; // <-- disable
      signalBinMax = -1; // <-- disable
    }

    inFile.setTransLength(nin);
    inFile.setConsLength(nin);
    iavg = 1;
    rnavg = 1.0 / navg;                            // initial block average
    return NORMAL;
  }

  /** underlying process operation for process()
   *  @return Command Status
   */
  public int _process () {
    int    stat    = NORMAL;
    double time0   = inFile.getTimeAt();
    int    numRead = inFile.read(inData);
    double time1   = inFile.getTimeAt();

    if (numRead < 0) {                             // EOF
      stat = FINISH;
      if (iavg == 1) return stat;                  // already a clean finish
      navg = iavg; rnavg = 1.0 / navg;             // prep for last partial average
      numRead = max(0, numRead);                   // finish off last frame
    }
    if (numRead == nin || partial != P_NONE) {
      Data statData = processInput(nfft, fftFile, inData, fftData, addData);

      if (statData != null) {
        if ((startTime != null) && (time1 < startTime)) {
          // too early
        }
        else if ((endTime   != null) && (time0 > endTime  )) {
          // too late
        }
        else {
          computeSigStatistics(statData);
          if (outResult != null) {
            MR.put(outResult, getStatTable());
          }
        }
      }
    }
    if (++iavg > navg) {                           // have we done enough for the next (averaged) output frame?
      iavg = 1;                                    // reset counter
      rnavg = 1.0/navg/nexp;                       // now start exp averaging
      if (out1d) stat = FINISH;                    // our work here is done
      if (partial == P_START) partial = P_NONE;    // don't want any more partials
    }
    return stat;
  }

  @Override
  public int processReady () {
    if (state != PROCESS) return 0;
    int n = inFile.processReady(nin);
    if (fftFile.isOpen) n = min(n, fftFile.processReady(1));
    return n;
  }

  @Override
  public int close () {
    if (verbose || ((outResult == null) && !outFile.isOpen())) {
      // since 3.7.1 using M.println instead of MT.writeln to allow logging - NXM-1813
      M.println("  NOISE_FLOOR              : "+roundTo3(statTable.get("NOISE_FLOOR"             ))+" dB\t"
              + "  AVG_NOISE_FLOOR          : "+roundTo3(statTable.get("AVG_NOISE_FLOOR"         ))+" dB");
      M.println("  NOISE_THRESHOLD          : "+roundTo3(statTable.get("NOISE_THRESHOLD"         ))+" dB\t"
              + "  AVG_NOISE_THRESHOLD      : "+roundTo3(statTable.get("AVG_NOISE_THRESHOLD"     ))+" dB");
      M.println("  ");
      M.println("  SIGNAL_STRENGTH          : "+roundTo3(statTable.get("SIGNAL_STRENGTH"         ))+" dB\t"
              + "  AVG_SIGNAL_STRENGTH      : "+roundTo3(statTable.get("AVG_SIGNAL_STRENGTH"     ))+" dB");
      M.println("  SIGNAL_SNR               : "+roundTo3(statTable.get("SIGNAL_SNR"              ))+" dB\t"
              + "  AVG_SIGNAL_SNR           : "+roundTo3(statTable.get("AVG_SIGNAL_SNR"          ))+" dB");
      M.println("  ");
      M.println("  PEAK_SIGNAL_FREQ         : "+roundMhz(statTable.get("PEAK_SIGNAL_FREQ"        ))+" MHz\t"
              + "  AVG_PEAK_SIGNAL_FREQ     : "+roundMhz(statTable.get("AVG_PEAK_SIGNAL_FREQ"    ))+" MHz");
      M.println("  PEAK_SIGNAL_STRENGTH     : "+roundTo3(statTable.get("PEAK_SIGNAL_STRENGTH"    ))+" dB\t"
              + "  AVG_PEAK_SIGNAL_STRENGTH : "+roundTo3(statTable.get("AVG_PEAK_SIGNAL_STRENGTH"))+" dB");
      M.println("  PEAK_SIGNAL_SNR          : "+roundTo3(statTable.get("PEAK_SIGNAL_SNR"         ))+" dB\t"
              + "  AVG_PEAK_SIGNAL_SNR      : "+roundTo3(statTable.get("AVG_PEAK_SIGNAL_SNR"     ))+" dB");
      M.println("  ");
      M.println("  SIGNAL_MEAN              : "+roundTo3(statTable.get("SIGNAL_MEAN"             ))+" dB\t"
              + "  SPECTRUM_MEAN            : "+roundTo3(statTable.get("SPECTRUM_MEAN"           ))+" dB");
      M.println("  SIGNAL_SDEV              : "+roundTo3(statTable.get("SIGNAL_SDEV"             ))+" dB\t"
              + "  SPECTRUM_SDEV            : "+roundTo3(statTable.get("SPECTRUM_SDEV"           ))+" dB");
      M.println("  SIGNAL_ELEM              : "+roundTo3(statTable.get("SIGNAL_ELEM"             ))+"\t"
              + "  SPECTRUM_ELEM            : "+roundTo3(statTable.get("SPECTRUM_ELEM"           ))+"");
      M.println("  SIGNAL_MIN               : "+roundTo3(statTable.get("SIGNAL_MIN"              ))+" dB\t"
              + "  SPECTRUM_MIN             : "+roundTo3(statTable.get("SPECTRUM_MIN"            ))+" dB");
      M.println("  SIGNAL_MAX               : "+roundTo3(statTable.get("SIGNAL_MAX"              ))+" dB\t"
              + "  SPECTRUM_MAX             : "+roundTo3(statTable.get("SPECTRUM_MAX"            ))+" dB");
      M.println("  SIGNAL_AMIN              : "+roundMhz(statTable.get("SIGNAL_AMIN"             ))+" MHz\t"
              + "  SPECTRUM_AMIN            : "+roundMhz(statTable.get("SPECTRUM_AMIN"           ))+" MHz");
      M.println("  SIGNAL_AMAX              : "+roundMhz(statTable.get("SIGNAL_AMAX"             ))+" MHz\t"
              + "  SPECTRUM_AMAX            : "+roundMhz(statTable.get("SPECTRUM_AMAX"           ))+" MHz");
      M.println("  SIGNAL_OMIN              : "+roundTo3(statTable.get("SIGNAL_OMIN"             ))+"\t"
              + "  SPECTRUM_OMIN            : "+roundTo3(statTable.get("SPECTRUM_OMIN"           ))+"");
      M.println("  SIGNAL_OMAX              : "+roundTo3(statTable.get("SIGNAL_OMAX"             ))+"\t"
              + "  SPECTRUM_OMAX            : "+roundTo3(statTable.get("SPECTRUM_OMAX"           ))+"");
      M.println("  SIGNAL_SUM               : "+roundTo3(statTable.get("SIGNAL_SUM"              ))+" dB\t"
              + "  SPECTRUM_SUM             : "+roundTo3(statTable.get("SPECTRUM_SUM"            ))+" dB");
    }
    if (outResult != null) {
      MR.put(outResult, getStatTable());
    }
    fft.free();
    if (out1d && fftFile.isOpen) fftFile.setFrameSize(0);
    inFile.close();
    outFile.close();
    fftFile.close();
    histSignal.clear();
    histPeak.clear();
    return NORMAL;
  }

  /** Round to 3 decimal places if numeric. */
  private static Object roundTo3 (Object obj) {
    if (obj instanceof Number) {
      obj = new Data(((Number)obj).doubleValue());
    }
    if (obj instanceof Data) {
      return padLeft(((Data)obj).toString("0.000"), 12);
    }
    return padLeft(Objects.toString(obj), 12);
  }

  /** Convert Hz to MHz and round to 3 decimal places if numeric. */
  private static Object roundMhz (Object obj) {
    if ((obj instanceof Number) || (obj instanceof Data)) {
      return roundTo3(Convert.o2d(obj) / 1e6);
    }
    return padLeft(Objects.toString(obj), 12);
  }

  /** Processes the input. This function is nearly-identical to the one in the
   *  FFT primitive.
   */
  private Data processInput (int numRead, DataFile testFile, Data inData, Data fftData, Data addData) {
    Data db;
    fftData.setSPA(inData.spa);                                                 // see dbf as input

    if (ac) {                                                                   // remove DC bias
      DataOp.sum(inData, acData, numRead);
      DataOp.scale(acData, acData, 1, (-1.0/numRead));
      DataOp.add(inData, acData, inData, numRead);
    }

    DataOp.mul(winData, inData, fftData);                                       // apply window

    if (numRead < nfft) DataOp.zero(fftData, numRead, nfft-numRead);            // zero fill

    if (cxm && inData.spa==1) DataOp.cmplx(fftData, fftData, nfft, DataOp.SC);  //force complex FFT

    fft.work(fftData);                                                          // perform FFTs
    fftData.setSPA(2); db=fftData;                                              // see dbf as output
    if (pack==P_ZERO) fftData.setF(1, 0.0F);

    if (mag) { DataOp.mag2(db, outData); db=outData; }                          // compute PSD

    if (nexp > 1) {
      DataOp.add(db, addData, addData);
      if (iavg < navg) return null;
      DataOp.scale(addData, outData, nout, rnavg);  // scale by 1/navg
      DataOp.scale(outData, addData, nout, rnexp);  // prescale by (1-1/nxp)/rnavg
      db=outData;
    }
    else if (navg > 1) {                            // average
      if (iavg == 1) DataOp.copy(db, addData);
      else DataOp.add(db, addData, addData);
      if (iavg < navg) return null;
      DataOp.scale(addData, addData, nout, rnavg);  // scale by 1/navg
      db=addData;
    }

    if      (fdb) DataOp.dbf(db, nout, 0);          // compute fast dB
    else if (log) DataOp.dblog(db, db);             // compute log

    if (nrot != 0) fft.rotate(db);                  // rotate complex fft bins
    db.setHeader( inData.getHeader() );             // propagate time line
    testFile.write(db);                             // write output
    return db;
  }

  /** Computes the signal statistics, only called by process(..) */
  private void computeSigStatistics (Data statData) {
    double[] dbuf = ( dp)? statData.castD(true) : null;
    float[]  fbuf = (!dp)? statData.castF(true) : null;
    SNRResult peakSNR;
    SNRResult tempSNR;
    SNRResult signalSNR;

    if (signalBinMin < 0) {
      peakSNR   = findSNR(dbuf, fbuf, 1, 0, nout);
      signalSNR = NO_RESULT;
    }
    else {
      peakSNR = findSNR(dbuf, fbuf, 1, rcvrBinMin,   rcvrBinMax);
      tempSNR = findSNR(dbuf, fbuf, 1, signalBinMin, signalBinMax); // The "floor" here is carrier+noise

      int noise  = peakSNR.noise;
      int signal = tempSNR.noise;

      signalSNR = new SNRResult(tempSNR.freq, signal, noise, peakSNR.thresh, "dB", tempSNR.stats);
    }

    if (dp) statData.uncast(dbuf, false);
    else    statData.uncast(fbuf, false);

    synchronized (this) {
      histSignal.add(signalSNR);
      histPeak.add(peakSNR);

      if (histSignal.size() > maxHistPoints) {
        histSignal.remove(0);
        histPeak.remove(0);
      }
    }

    SNRResult avgPeak   = avgResult(histPeak);
    SNRResult avgSignal = avgResult(histSignal);

    statTable_put("NOISE_FLOOR",              "NFLR", peakSNR.noise   );
    statTable_put("NOISE_THRESHOLD",          "NTHR", peakSNR.thresh  );
    statTable_put("SIGNAL_STRENGTH",          "SIG",  signalSNR.signal);
    statTable_put("SIGNAL_SNR",               "SSNR", signalSNR.snr   );
    statTable_put("PEAK_SIGNAL_FREQ",         "FREQ", peakSNR.freq    );
    statTable_put("PEAK_SIGNAL_STRENGTH",     "PEAK", peakSNR.signal  );
    statTable_put("PEAK_SIGNAL_SNR",          "PSNR", peakSNR.snr     );
    statTable_put("AVG_NOISE_FLOOR",          "ANF",  avgPeak.noise   );
    statTable_put("AVG_NOISE_THRESHOLD",      "ANT",  avgPeak.thresh  );
    statTable_put("AVG_SIGNAL_STRENGTH",      "ASIG", avgSignal.signal);
    statTable_put("AVG_SIGNAL_SNR",           "ASNR", avgSignal.snr   );
    statTable_put("AVG_PEAK_SIGNAL_FREQ",     "AFRQ", avgPeak.freq    );
    statTable_put("AVG_PEAK_SIGNAL_STRENGTH", "APSG", avgPeak.signal  );
    statTable_put("AVG_PEAK_SIGNAL_SNR",      "APSN", avgPeak.snr     );

    if (avgSignal.stats != null) {
      statTable_put("SIGNAL_MEAN",  null, avgSignal.stats.mean);
      statTable_put("SIGNAL_SDEV",  null, avgSignal.stats.stdDev);
      statTable_put("SIGNAL_ELEM",  null, avgSignal.stats.nelem);
      statTable_put("SIGNAL_MIN",   null, avgSignal.stats.min);
      statTable_put("SIGNAL_MAX",   null, avgSignal.stats.max);
      statTable_put("SIGNAL_AMIN",  null, avgSignal.stats.minAbs);
      statTable_put("SIGNAL_AMAX",  null, avgSignal.stats.maxAbs);
      statTable_put("SIGNAL_OMIN",  null, avgSignal.stats.minIndex);
      statTable_put("SIGNAL_OMAX",  null, avgSignal.stats.maxIndex);
      statTable_put("SIGNAL_SUM",   null, avgSignal.stats.sum);
    }

    if (avgPeak.stats != null) {
      statTable_put("SPECTRUM_MEAN",  null, avgPeak.stats.mean);
      statTable_put("SPECTRUM_SDEV",  null, avgPeak.stats.stdDev);
      statTable_put("SPECTRUM_ELEM",  null, avgPeak.stats.nelem);
      statTable_put("SPECTRUM_MIN",   null, avgPeak.stats.min);
      statTable_put("SPECTRUM_MAX",   null, avgPeak.stats.max);
      statTable_put("SPECTRUM_AMIN",  null, avgPeak.stats.minAbs);
      statTable_put("SPECTRUM_AMAX",  null, avgPeak.stats.maxAbs);
      statTable_put("SPECTRUM_OMIN",  null, avgPeak.stats.minIndex);
      statTable_put("SPECTRUM_OMAX",  null, avgPeak.stats.maxIndex);
      statTable_put("SPECTRUM_SUM",   null, avgPeak.stats.sum);
    }
    if (outFile.isOpen()) {
      if (statData.getHeader() instanceof Time) {
        outFile.setTimeAt((Time)statData.getHeader());
      }
      outFile.writeDataTable(statTable);
    }
  }

  /** Computes the average value of the set of results. */
  private SNRResult avgResult (List<SNRResult> res) {
    switch (res.size()) {
      case 0: return NO_RESULT;
      case 1: return res.get(0);
    }

    double  freq     = 0;
    double  signal   = 0;
    double  noise    = 0;
    double  thresh   = 0;
    double  rmean    = 0;
    double  imean    = 0;
    double  stdDev   = 0;
    double  min      = 0;
    double  max      = 0;
    int     minIndex = 0;
    int     maxIndex = 0;
    int     count    = 0;
    Double  rsum     = null;
    Double  isum     = null;
    Integer nelem    = null;
    Boolean complex  = null;
    boolean hasStats = true;
    Stats   stats    = null;

    for (SNRResult r : res) {
      freq   += r.freq;
      signal += r.signal;
      noise  += r.noise;
      thresh += r.thresh;
      count++;

      if (!hasStats || (r.stats == null)) {
        hasStats = false;
      }
      else if (count == 1) {
        rmean    = r.stats.rmean;
        imean    = r.stats.imean;
        rsum     = r.stats.rsum;
        isum     = r.stats.isum;
        stdDev   = r.stats.stdDev;
        nelem    = r.stats.nelem;
        complex  = r.stats.complex;
        min      = r.stats.min;
        max      = r.stats.max;
        minIndex = r.stats.minIndex;
        maxIndex = r.stats.maxIndex;
      }
      else {
        rmean   += r.stats.rmean;
        imean   += r.stats.imean;
        stdDev  += r.stats.stdDev;
        rsum     = ((rsum == null) || (r.stats.rsum == null))? null : rsum + r.stats.rsum;
        isum     = ((isum == null) || (r.stats.isum == null))? null : isum + r.stats.isum;

        if (((nelem == null) && (r.stats.nelem != null)) ||
            ((nelem != null) && ((r.stats.nelem == null) ||
                                 (nelem.intValue() != r.stats.nelem.intValue())))) {
          hasStats = false; // illegal mix of frame sizes
        }
        if (((complex == null) && (r.stats.complex != null)) ||
            ((complex != null) && ((r.stats.complex == null) ||
                                   (complex.booleanValue() != r.stats.complex.booleanValue())))) {
          hasStats = false; // illegal mix of real and complex values
        }
        if (min > r.stats.min) {
          min      = r.stats.min;
          minIndex = r.stats.minIndex;
        }
        if (max < r.stats.max) {
          max      = r.stats.max;
          maxIndex = r.stats.maxIndex;
        }
      }
    }

    if (hasStats) {
      double minAbs = fftFreqMin + rcvrTF + (minIndex * fftFreqDelta);
      double maxAbs = fftFreqMin + rcvrTF + (maxIndex * fftFreqDelta);

      stats = new Stats(rmean/count, imean/count, rsum, isum, stdDev/count, min, max, minAbs,
                        maxAbs, minIndex, maxIndex, nelem, complex);
    }

    return new SNRResult((int)rint(freq   / count),
                         (int)rint(signal / count),
                         (int)rint(noise  / count),
                         (int)rint(thresh / count),
                         "dB", stats);
  }

  /** Used by computeSigStatistics(..) to insert values into the statTable. */
  private void statTable_put (String t6key, String t3key, Object val) {
    if (val != null) {
      if (t6key != null) statTable.put(t6key, val);
      if (t3key != null) statTable.put(t3key, val);
    }
    else {
      if (t6key != null) statTable.remove(t6key);
      if (t3key != null) statTable.remove(t3key);
    }
  }

  /** Finds the SNR in the spectrum using statistical methods.
   *  @param dbuf  Buffer if double-precision (either dbuf or fbuf must be set).
   *  @param fbuf  Buff if single-precision.
   *  @param spa   Complex mode (scalar=1, complex=2).
   *  @param first The index of the first atom to read in (inclusive).
   *  @param last  The index of the last atom to read in (exclusive).
   *  @return The SNR results.
   */
  private SNRResult findSNR (double[] dbuf, float[] fbuf, int spa, int first, int last) {
    final boolean FUNC_DEBUG  = false;
    final double  CLIP_WEIGHT = 1.0;
    final double  CLIP_DELTA  = 0.5;

    Stats   initStats   = (dbuf != null) ? computeStatistics(dbuf, spa, NEGATIVE_INFINITY, POSITIVE_INFINITY, CLIP_WEIGHT, first, last)
                                         : computeStatistics(fbuf, spa, NEGATIVE_INFINITY, POSITIVE_INFINITY, CLIP_WEIGHT, first, last);
    Stats   lastStats   = initStats;
    double  peakFreq    = (initStats.maxIndex * fftFreqDelta) + fftFreqMin + rcvrTF;
    double  peak        = initStats.max;
    double  upperClip   = initStats.max;
    double  lowerClip   = initStats.min;
    double  upperFloor  = initStats.max;
    double  lowerFloor  = initStats.min;

    // FIND THE NOISE THRESHOLD (this is intended to discard any effects of receiver/A2D roll-off)
    while (true) {
      Stats thisStats = (dbuf != null) ? computeStatistics(dbuf, spa, lowerClip, upperClip, CLIP_WEIGHT, first, last)
                                       : computeStatistics(fbuf, spa, lowerClip, upperClip, CLIP_WEIGHT, first, last);
      if (thisStats == null) {
        break;  // This means that all values were above/below the clip
      }
      else {
        double dStdDev  = thisStats.stdDev  - lastStats.stdDev;
        double dMean    = thisStats.rmean   - lastStats.rmean;

        if (FUNC_DEBUG) {
          System.out.format("%n");
          System.out.format("v  clip  = [%8.3f, %8.3f]%n", lowerClip, upperClip);
          System.out.format("v  rmean = %8.3f  stdDev = %8.3f%n ", thisStats.rmean, thisStats.stdDev);
          System.out.format("v  dmean = %8.3f  dStdDev  = %8.3f%n ", dMean, dStdDev);
        }

        if ((abs(dMean) > abs(dStdDev)) || (thisStats.min >= thisStats.max)) {
          break;
        }
      }
      lowerFloor  = lowerClip;
      lowerClip   = lowerClip + CLIP_DELTA;
      lastStats   = thisStats;
    }

    if (FUNC_DEBUG) {
      System.out.format("%n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%n");
    }

    // FIND THE NOISE FLOOR
    lastStats = initStats;

    while (true) {
      Stats thisStats = (dbuf != null) ? computeStatistics(dbuf, spa, lowerClip, upperClip, CLIP_WEIGHT, first, last)
                                       : computeStatistics(fbuf, spa, lowerClip, upperClip, CLIP_WEIGHT, first, last);
      if (thisStats == null) {
        break; // This means that all values were above/below the clip
      }
      else {
        double dStdDev = thisStats.stdDev - lastStats.stdDev;
        double dMean   = thisStats.rmean  - lastStats.rmean;

        if (FUNC_DEBUG) {
          System.out.format("%n");
          System.out.format("^  clip  = [%8.3f, %8.3f]%n", lowerClip, upperClip);
          System.out.format("^  rmean = %8.3f  stdDev   = %8.3f%n ", thisStats.rmean, thisStats.stdDev);
          System.out.format("^  dmean = %8.3f  dStdDev  = %8.3f%n ", dMean,
              dStdDev);
        }

        if ((abs(dMean) > abs(dStdDev)) || (thisStats.min >= thisStats.max)) {
          break;
        }
      }
      upperFloor  = upperClip;
      upperClip   = upperClip - CLIP_DELTA;
      lastStats   = thisStats;
    }

    if (FUNC_DEBUG) {
      System.out.format("%n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%n");
      System.out.format("lowerFloor = %8.3f%n", lowerFloor);
      System.out.format("upperFloor = %8.3f%n", upperFloor);
    }

    return _snr(peak, upperFloor, lowerFloor, peakFreq, initStats);
  }

  /** Get the number of input frames to block average for each output frame.
   *  @return number of frames to average
   */
  public int getNAvg () {
    return navg;
  }

  /** Get the number of output frames for exponentially average (sliding avg) for each output frame.
   *  @return number of frames to exponentially average
   */
  public int getNExp () {
    return nexp;
  }

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

  /** get the fractional overlap, zero for no overlap, negative to skip.
      @return overlap
   */
  public double getOverlap () {
    return over;
  }

  /** Get the windowing function name.
   *  @return name of window function used
   */
  public String getWindow () {
    return wname;
  }

  /** Set the number of input frames to block average for each output frame.
   *  @param value number of frames to average
   */
  public void setNAvg (int value) {
    navg  = max(1, value);
    rnavg = 1.0/navg/nexp;
  }

  /** Set the number of output frames to exponentially average (sliding avg) for each output frame.
   *  @param value number of frames to exponentially average
   */
  public void setNExp (int value) {
    nexp  = max(1, value);
    rnexp = (1.0-(1.0/nexp)) * nexp * navg;
  }

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

  /** Set the fractional overlap, zero for no overlap, negative to skip.
   *  @param value overlap to use
   */
  public void setOverlap (double value) {
    if (in2d && stride%nin != 0) M.warning("Type 2000 inputs don't have integral overlap="+over);
    over   = value;       // note stride of 0 can be useful
    stride = max(0,(int)((1.0-over)*nin));
    inFile.setXferLength(nin);
    inFile.setConsLength(stride);
  }

  /** Set the windowing function.
   *  @param value windowing function to use
   */
  public void setWindow (String value) {
    int nwind = min(nin, nfft);
    wname     = value;
    winData   = DataOp.getWindow(wname, nwind, factor, M);
  }

  /** Gets a copy of the current statistics.
   *  @return table with current statistics
   */
  public synchronized Table getStatTable () {
    return statTable.copy();
  }

  /** Gets the most recent noise floor (NOISE_FLOOR) in dB.
   *  @return current noise floor
   */
  public int getNoiseFloor () {
    return statTable.getL("NOISE_FLOOR");
  }

  /** Gets the average noise floor (AVG_NOISE_FLOOR) in dB.
   *  @return average noise floor
   */
  public int getAvgNoiseFloor () {
    return statTable.getL("AVG_NOISE_FLOOR");
  }

  /** Gets the most recent noise threshold (NOISE_THRESHOLD) in dB.
   *  @return current noise threshold
   */
  public int getNoiseThreshold () {
    return statTable.getL("NOISE_THRESHOLD");
  }

  /** Gets the average noise threshold (AVG_NOISE_THRESHOLD) in dB.
   *  @return average noise threshold
   */
  public int getAvgNoiseThreshold () {
    return statTable.getL("AVG_NOISE_THRESHOLD");
  }

  /** Gets the most recent signal strength (SIGNAL_STRENGTH) in dB.
   *  @return current signal strength
   */
  public int getSignalStrength () {
    return statTable.getL("SIGNAL_STRENGTH");
  }

  /** Gets the average signal strength (AVG_SIGNAL_STRENGTH) in dB.
   *  @return average signal strength
   */
  public int getAvgSignalStrength () {
    return statTable.getL("AVG_SIGNAL_STRENGTH");
  }

  /** Gets the most recent signal SNR (SIGNAL_SNR) in dB.
   *  @return current signal signal to noise ratio
   */
  public int getSignalSNR () {
    return statTable.getL("SIGNAL_SNR");
  }

  /** Gets the average signal SNR (AVG_SIGNAL_SNR) in dB.
   *  @return average signal signal to noise ratio
   */
  public int getAvgSignalSNR () {
    return statTable.getL("AVG_SIGNAL_SNR");
  }

  /** Gets the most recent peak signal frequency (PEAK_SIGNAL_FREQ) in Hz.
   *  @return current peak signal frequency
   */
  public int getPeakSignalFreq () {
    return statTable.getL("PEAK_SIGNAL_FREQ");
  }

  /** Gets the average signal frequency (AVG_PEAK_SIGNAL_FREQ) in Hz (returns 0 if it jumps around).
   *  @return average peak signal frequency
   */
  public int getAvgPeakSignalFreq () {
    return statTable.getL("AVG_PEAK_SIGNAL_FREQ");
  }

  /** Gets the most recent peak signal strength (PEAK_SIGNAL_STRENGTH) in dB.
   *  @return current peak signal strength
   */
  public int getPeakSignalStrength () {
    return statTable.getL("PEAK_SIGNAL_STRENGTH");
  }

  /** Gets the average peak signal strength (AVG_PEAK_SIGNAL_STRENGTH) in dB.
   *  @return average peak signal strength
   */
  public int getAvgPeakSignalStrength () {
    return statTable.getL("AVG_PEAK_SIGNAL_STRENGTH");
  }

  /** Gets the most recent peak signal SNR (PEAK_SIGNAL_SNR) in dB.
   *  @return current peak signal SNR
   */
  public int getPeakSignalSNR () {
    return statTable.getL("PEAK_SIGNAL_SNR");
  }

  /** Gets the average signal SNR (AVG_PEAK_SIGNAL_SNR) in dB.
   *  @return average peak signal SNR
   */
  public int getAvgPeakSignalSNR () {
    return statTable.getL("AVG_PEAK_SIGNAL_SNR");
  }

  /** Computes the SNR of a signal. When doing this computation based on a PSD,
   *  the assumption here is that noise is uniform across the bins, the signal is
   *  contained within a single bin and the signal bandwidth (in Hz) is equal to
   *  the width of the bin (in Hz). In effect, this is identical to a special
   *  case of the equation used in {@link #_snr(double,double,double,double)}
   *  where <code>Sb=1</code> and <code>Cb=1</code>.
   *  <pre>
   *  What can you measure in the spectrum:
   *    Sp = Cp + Np  = (carrier + noise)         in dB
   *    Np = (noise)                              in dB
   *
   *  What you can't measure in the spectrum:
   *    Cp = (carrier)                            in dB
   *
   *  What you know (omits the reference signal since it cancels out):
   *              (Sp / 10)
   *    Sp' = 10                                  not in dB
   *              (Np / 10)
   *    Np' = 10                                  not in dB
   *              (Cp / 10)
   *    Cp' = 10                                  not in dB
   *
   *  What you need:
   *    SNR = 10 log (Cp' / Np') = Cp - Np        in dB
   *
   *  Solving:
   *                    (SpNp/10)     (Np/10)
   *                  10          - 10
   *    SNR = 10 log (-----------------)          in dB [This is the version used here]
   *                        (Np/10)
   *                      10
   *
   *                    (Sp - Np)/10
   *    SNR = 10 log (10            - 1)        in dB [This version is in common use]
   *  </pre>
   *  @param Sp  The power of the carrier plus the noise floor in dB. (This is
   *             what is shown in the PSd within the BW of the signal.)
   *  @param Np  The noise floor in dB.
   *  @return The SNR in dB.
   */
  public static int snr (int Sp, int Np) {
    return _snr(Sp, Np, 0.0, 0.0, null).snr;
  }

  /** Internal use only: Version of the snr(..) function that takes in noise threshold (Nt) and freq.
   *  @param Sp  The power of the carrier plus the noise floor in dB. (This is
   *             what is shown in the PSd within the BW of the signal.)
   *  @param Np  The noise floor in dB.
   *  @param Nt  The noise lower threshold in dB.
   *  @param freq frequency
   *  @return The SNR in dB.
   *  @deprecated Since NeXtMidas 4.1.4 -- Use the version that takes in a {@code stats} argument.
   */
  @InternalUseOnly
  @Deprecated
  public static SNRResult _snr (double Sp, double Np, double Nt, double freq) {
    return _snr(Sp, Np, Nt, freq, null);
  }

  /** Internal use only: Version of the snr(..) function that takes in noise threshold (Nt) and freq.
   *  @param Sp    The power of the carrier plus the noise floor in dB. (This is
   *               what is shown in the PSd within the BW of the signal.)
   *  @param Np    The noise floor in dB.
   *  @param Nt    The noise lower threshold in dB.
   *  @param freq  Frequency in Hz.
   *  @param stats Statistics to include with the result (can be null).
   *  @return The SNR in dB.
   *  @since NeXtMidas 4.1.4
   */
  @InternalUseOnly
  public static SNRResult _snr (double Sp, double Np, double Nt, double freq, Stats stats) {
    if (Sp > Np) {
      double cn = pow(10, (Sp)/10.0);   // (carrier + noise)
      double n  = pow(10, (Np)/10.0);   // (noise)
      double t  = pow(10, (Nt)/10.0);   // (noise - lower threshold)
      return new SNRResult(freq, cn-n, n, t, stats);
    }
    else {
      // If Sp = Np then there is no signal here, just noise.   => SNR = 0 dB
      //
      // If Sp < Np our "noise is uniform" assumption has failed
      // and we should just consider anything here to be noise. => SNR = 0 dB
      //
      // [J. Czechowski] One could argue that Sp < Np represents destructive interference
      // between the signal and the noise; but the odds of that being the case
      // (as opposed to noise not being uniform) are near nil, and anyone who
      // needs that functionality should probably use a custom tool to measure
      // that, I see no good reason to add support for that here.
      double n = pow(10, (Np)/10.0);    // (noise)
      double t = pow(10, (Nt)/10.0);    // (noise - lower threshold)
      return new SNRResult(freq, n, n, t, stats);
    }
  }

  /** Compute the statistics for a given set of data.
   *  @param buf        The data read in (cast to double).
   *  @param spa        The number of scalars per atom (1=Scalar, 2=Complex, ...).
   *  @param upperClip  The upper clip value (+Inf to disable). Does not affect maximum value.
   *  @param lowerClip  The lower clip value (-Inf to disable). Does not affect minimum value.
   *  @param clipWeight The weight given to a value outside the clip region when calculating
   *                    mean and standard deviation (1=Include, 0=Omit).
   *  @param first      The index of the first atom to read in (inclusive).
   *  @param last       The index of the last atom to read in (exclusive).
   *  @return The statistics.
   */
  private Stats computeStatistics (double[] buf, int spa, double lowerClip, double upperClip,
                                   double clipWeight, int first, int last) {
    // THIS IS 100% IDENTICAL TO THE float[] VERSION EXCEPT FOR THE TYPE OF BUFFER PASSED IN
    double  rsum     = 0;
    double  isum     = 0;
    double  stdDev   = 0;
    double  min      = POSITIVE_INFINITY;
    double  max      = NEGATIVE_INFINITY;
    double  count    = 0;
    int     nelem    = last - first;
    int     minIndex = -1;
    int     maxIndex = -1;
    boolean complex = (spa == 2);
    double  upperClipSum = (upperClip * clipWeight);
    double  lowerClipSum = (lowerClip * clipWeight);
    double  upperClipSD  = upperClipSum * upperClipSum;  // => (upperClip * clipWeight)^2
    double  lowerClipSD  = lowerClipSum * lowerClipSum;  // => (lowerClip * clipWeight)^2

    if (spa == 2) {
      // This isn't the correct implementation for clipping....
      //   for (int i = first*spa; i< last*spa; i+=2) {
      //     double rval = buf[ i ];
      //     double ival = buf[i+1];
      //     double abs  = rval*rval + ival*ival;
      //     if (abs > clip) continue; // skip it
      //     rsum   += rval;
      //     isum   += ival;
      //     stdDev += abs;
      //     count  += 1;
      //     // leave dmin/dmax as squared until end for speed
      //     if (abs > max) { max = abs; maxIndex = i; }
      //     if (abs < min) { min = abs; minIndex = i; }
      //   }
    }
    else {
      for (int i = first*spa; i < last*spa; i++) {
        double val = buf[i];

        if (val > upperClip) {
          rsum   += upperClipSum;
          stdDev += upperClipSD;
          count  += clipWeight;
        }
        else if (val < lowerClip) {
          rsum   += lowerClipSum;
          stdDev += lowerClipSD;
          count  += clipWeight;
        }
        else {
          double abs = val*val;
          rsum   += val;
          stdDev += abs;
          count  += 1;
        }
        if (val > max) { max = val; maxIndex = i; }
        if (val < min) { min = val; minIndex = i; }
      }
    }

    // finish computation
    double rmean = rsum / count;
    double imean = isum / count;
    stdDev = stdDev/count - rmean*rmean - imean*imean;
    stdDev = (stdDev < 0.0) ? 0.0 : sqrt(stdDev);

    // revert complex squared back to magnitude
    if (complex) {
      min = sqrt(min);
      max = sqrt(max);
    }
    double minAbs = fftFreqMin + rcvrTF + (minIndex * fftFreqDelta);
    double maxAbs = fftFreqMin + rcvrTF + (maxIndex * fftFreqDelta);

    return (count == 0) ? null : new Stats(rmean, imean, rsum, isum,  stdDev, min, max,
                                           minAbs, maxAbs, minIndex, maxIndex, nelem, complex);
  }

  /** Compute the statistics for a given set of data.
   *  @param buf        The data read in (cast to float).
   *  @param spa        The number of scalars per atom (1=Scalar, 2=Complex, ...).
   *  @param upperClip  The upper clip value (+Inf to disable). Does not affect maximum value.
   *  @param lowerClip  The lower clip value (-Inf to disable). Does not affect minimum value.
   *  @param clipWeight The weight given to a value outside the clip region when calculating
   *                    mean and standard deviation (1=Include, 0=Omit).
   *  @param first      The index of the first atom to read in (inclusive).
   *  @param last       The index of the last atom to read in (exclusive).
   *  @return The statistics.
   */
  private Stats computeStatistics (float[] buf, int spa, double lowerClip, double upperClip,
                                   double clipWeight, int first, int last) {
    // THIS IS 100% IDENTICAL TO THE double[] VERSION EXCEPT FOR THE TYPE OF BUFFER PASSED IN
    double  rsum     = 0;
    double  isum     = 0;
    double  stdDev   = 0;
    double  min      = POSITIVE_INFINITY;
    double  max      = NEGATIVE_INFINITY;
    double  count    = 0;
    int     nelem    = last - first;
    int     minIndex = -1;
    int     maxIndex = -1;
    boolean complex = (spa == 2);
    double  upperClipSum = (upperClip * clipWeight);
    double  lowerClipSum = (lowerClip * clipWeight);
    double  upperClipSD  = upperClipSum * upperClipSum;  // => (upperClip * clipWeight)^2
    double  lowerClipSD  = lowerClipSum * lowerClipSum;  // => (lowerClip * clipWeight)^2

    if (spa == 2) {
      // This isn't the correct implementation for clipping....
      //   for (int i = first*spa; i< last*spa; i+=2) {
      //     double rval = buf[ i ];
      //     double ival = buf[i+1];
      //     double abs  = rval*rval + ival*ival;
      //     if (abs > clip) continue; // skip it
      //     rsum   += rval;
      //     isum   += ival;
      //     stdDev += abs;
      //     count  += 1;
      //     // leave dmin/dmax as squared until end for speed
      //     if (abs > max) { max = abs; maxIndex = i; }
      //     if (abs < min) { min = abs; minIndex = i; }
      //   }
    }
    else {
      for (int i = first*spa; i < last*spa; i++) {
        double val = buf[i];

        if (val > upperClip) {
          rsum   += upperClipSum;
          stdDev += upperClipSD;
          count  += clipWeight;
        }
        else if (val < lowerClip) {
          rsum   += lowerClipSum;
          stdDev += lowerClipSD;
          count  += clipWeight;
        }
        else {
          double abs = val*val;
          rsum   += val;
          stdDev += abs;
          count  += 1;
        }
        if (val > max) { max = val; maxIndex = i; }
        if (val < min) { min = val; minIndex = i; }
      }
    }

    // finish computation
    double rmean = rsum / count;
    double imean = isum / count;
    stdDev = stdDev/count - rmean*rmean - imean*imean;
    stdDev = (stdDev < 0.0)? 0.0 : sqrt(stdDev);

    // revert complex squared back to magnitude
    if (complex) {
      min = sqrt(min);
      max = sqrt(max);
    }
    double minAbs = fftFreqMin + rcvrTF + (minIndex * fftFreqDelta);
    double maxAbs = fftFreqMin + rcvrTF + (maxIndex * fftFreqDelta);

    return (count == 0)? null : new Stats(rmean, imean, rsum, isum, stdDev, min, max,
                                          minAbs, maxAbs, minIndex, maxIndex, nelem, complex);
  }

  //////////////////////////////////////////////////////////////////////////////
  // NESTED CLASSES
  //////////////////////////////////////////////////////////////////////////////

  /** The statistical information about the data. */
  @InternalUseOnly
  public static final class Stats implements Tablizable {
    final double  rmean;     // real mean
    final double  imean;     // complex mean
    final double  stdDev;    // standard deviation
    final double  min;       // min value
    final double  max;       // max value
    final int     minIndex;  // index of min value (-1 if n/a)
    final int     maxIndex;  // index of max value (-1 if n/a)
    // Added in NeXtMidas 4.1.4...
    final Double  rsum;      // real sum
    final Double  isum;      // complex sum
    final Double  minAbs;    // min abscissa value
    final Double  maxAbs;    // max abscissa value
    final Integer nelem;     // number of elements
    final Boolean complex;   // complex (true) or real (false)
    final Data    sum;       // sum of values
    final Data    mean;      // mean of values

    /** Creates a new instance.
     *  @param rmean    Real mean
     *  @param imean    Complex mean
     *  @param stdDev   Standard deviation
     *  @param min      Min value
     *  @param max      Max value
     *  @param minIndex Index of min value (-1 if n/a)
     *  @param maxIndex Index of max value (-1 if n/a)
     */
    public Stats (double rmean, double imean, double stdDev, double min, double max,
                  int minIndex, int maxIndex) {
      this(rmean, imean, null, null, stdDev, min, max, null, null, minIndex, maxIndex, null, null);
    }

    /** Creates a new instance.
     *  @param rmean    Real mean
     *  @param imean    Complex mean
     *  @param rsum     Real sum
     *  @param isum     Complex sum
     *  @param stdDev   Standard deviation
     *  @param min      Min value
     *  @param max      Max value
     *  @param minAbs   Min abscissa value
     *  @param maxAbs   Max abscissa value
     *  @param minIndex Index of min value (-1 if n/a)
     *  @param maxIndex Index of max value (-1 if n/a)
     *  @param nelem    Number of elements
     *  @param complex  Is the data complex (true) or real (false)?
     */
    Stats (double rmean, double imean, Double rsum, Double isum, double stdDev,
           double min, double max, Double minAbs, Double maxAbs, int minIndex,
           int maxIndex, Integer nelem, Boolean complex) {
      this.rmean    = rmean;
      this.imean    = imean;
      this.rsum     = rsum;
      this.isum     = isum;
      this.stdDev   = stdDev;
      this.min      = min;
      this.max      = max;
      this.minAbs   = minAbs;
      this.maxAbs   = maxAbs;
      this.minIndex = minIndex;
      this.maxIndex = maxIndex;
      this.nelem    = nelem;
      this.complex  = complex;

      if (complex == null) {
        sum  = null;
        mean = null;
      }
      else if (complex) {
        sum  = new Data("CD", 1);
        mean = new Data("CD", 1);

        sum.setD(0,rsum);   if (complex) sum.setD(1,isum);
        mean.setD(0,rmean); if (complex) mean.setD(1,imean);
      }
      else {
        sum  = new Data("SD", 1);
        mean = new Data("SD", 1);
        sum.setD(0,rsum);
        mean.setD(0,rmean);
      }
    }

    @Override
    public String toString () {
      return toTable().toString();
    }

    @Override
    public Table toTable () { // Since NeXtMidas 4.1.4
      Table table = new Table();

      // The ordering here matches STATISTICS (though it doesn't need to)
      if (mean   != null) table.put("MEAN", mean);
      if (true          ) table.put("SDEV", stdDev);
      if (nelem  != null) table.put("ELEM", nelem);
      if (true          ) table.put("MIN",  min);
      if (true          ) table.put("MAX",  max);
      if (minAbs != null) table.put("AMIN", minAbs);
      if (maxAbs != null) table.put("AMAX", maxAbs);
      if (true          ) table.put("OMIN", minIndex);
      if (true          ) table.put("OMAX", maxIndex);
      if (sum    != null) table.put("SUM",  sum);

      return table;
    }
  }

  /** The results of an SNR computation. */
  @InternalUseOnly
  public static final class SNRResult implements Tablizable {
    final int   freq;   // in Hz
    final int   signal; // in dB
    final int   noise;  // in dB
    final int   thresh; // in dB
    final int   snr;    // in dB
    final Stats stats;  // Statistics (since NeXtMidas 4.1.4)

    /** Creates a new instance with all values set to 0.
     *  @deprecated Since NeXtMidas 4.1.4 -- This is for internal use and should not be public.
     */
    @Deprecated // Change to package-private (this is special-purpose and only used by NO_RESULT).
    public SNRResult () {
      this(0, 0, 0, 0, "dB", null);
    }

    /** Creates a new instance.
     *  @param freq   Frequency in Hz.
     *  @param signal Signal in dB/Hz.
     *  @param noise  Noise in dB/Hz.
     *  @param thresh Threshold in dB/Hz.
     *  @deprecated Since NeXtMidas 4.1.4 -- This was for internal use and should not have been public.
     */
    @Deprecated
    public SNRResult (double freq, double signal, double noise, double thresh) {
      this(freq, signal, noise, thresh, null);
    }

    /** Creates a new instance.
     *  @param freq   Frequency in Hz.
     *  @param signal Signal in dB/Hz.
     *  @param noise  Noise in dB/Hz.
     *  @param thresh Threshold in dB/Hz.
     *  @param stats  Additional statistics (can be null).
     */
    SNRResult (double freq, double signal, double noise, double thresh, Stats stats) {
      this.freq   = (int)rint(freq);                        // in Hz
      this.signal = (int)rint(10 * log10(signal));          // in dB
      this.noise  = (int)rint(10 * log10(noise));           // in dB
      this.thresh = (int)rint(10 * log10(thresh));          // in dB
      this.snr    = (int)rint(10 * log10(signal / noise));  // in dB
      this.stats  = stats;
    }

    /** Creates a new instance.
     *  @param freq   Frequency in Hz.
     *  @param signal Signal in dB.
     *  @param noise  Noise in dB.
     *  @param thresh Threshold in dB.
     *  @param units  This must be equal to "dB".
     *  @param stats  Additional statistics (can be null).
     */
    SNRResult (int freq, int signal, int noise, int thresh, String units, Stats stats) {
      if (!units.equals("dB")) {
        throw new IllegalArgumentException("Expected units=dB but found units="+units);
      }
      this.freq   = freq;
      this.signal = signal;
      this.noise  = noise;
      this.thresh = thresh;
      this.snr    = signal - noise;
      this.stats  = stats;
    }

    @Override
    public String toString () {
      return toTable().toString();
    }

    @Override
    public Table toTable () { // Since NeXtMidas 4.1.4
      Table table = (stats != null)? stats.toTable() : new Table();

      table.put("FREQ",   freq);
      table.put("SIGNAL", signal);
      table.put("NOISE",  noise);
      table.put("THRESH", thresh);
      table.put("SNR",    snr);

      return table;
    }
  }
}
