package nxm.dsp.prim;

import java.text.NumberFormat;
import java.util.StringTokenizer;

import nxm.sys.lib.BaseFile;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.DataOp;
import nxm.sys.lib.ICEPacket;
import nxm.sys.lib.Primitive;
import nxm.sys.lib.Table;
import nxm.sys.lib.Time;
import nxm.sys.libm.Fft;
import nxm.sys.libm.Window;

/** Primitive which takes a time domain input signal
    and produces the cyclic spectrum at specified
    cycle frequencies. The primitive has been optimized
    using phase rotations of the frequency bins to avoid
    multiple FFT's and LO computation.
    <br><br>

    Unlike FASTCYCLO, this primitive uses implements a
    2-Stage Accumulator. Stage-1 uses Coherent Integration
    and Stage-2 uses Normalized Integration. Stage-1 feeds
    Stage-2 in serial fashion. The number of integrations
    for each stage is specified using the "|" separator in
    the &lt;INTEGRATION&gt; command line parameter.

    @author RGH (X-Midas host primitive)
    @author  J. Czechowski IV (NeXtMidas Java port)
    @version $Id: cyclofunc.java,v 1.6 2014/10/17 15:03:32 ntn Exp $
    @since NeXtMidas 3.2.0
 */
/*
data arrays:
  grab        - from file, nfft at a time
  psd         - autospectrum of input, optionally normalized by 1/sqrt(max(1e-20,sumsq/nsub))
  sxx         - autospectrum of input, optionally normalized by square of above
  win         - window values
  wddat       - windowed data before and after FFT

per cycle freq:
  syy         - magnitude of wddat convolved with cxlo if multiple coherent accumulations
  csxy        - coherent accumulation area
  cnmsxy      - non-coherent accumulation area
  cxlo        - centered response for oscillator of this freq

other data per cycle freq:
  blkphas     - maintains phase from grab to grab
  freq        - this cycle frequency

data flow:

                         grab (r,c)      win (r)
                            \            /
                             win_stat_rot()
                                   |
                        wddat (r,c)     - complex data is rotated (odd samples sign
                                   |            reversed)
                                   |
                                  fft()
                                   |
                                   |
                  ---------------wddat (c)      - rotate by one bin if the input is complex
                 /                 |
                /           accum_sxx_psd()
               /           /              \
              /          psd (r)          sxx (r)    both get autospectrum
             /  (opt) scale by const   (opt) scale by square of that const
            |             |                   \
            |          psd file                \
            |                                   \
            |   for each cycle:                  \-----------------------\
            |                                                             \
            \                  angle = cyc[i].blkphase*2pi                 \
             \                 phase = (cos(angle),isin(angle)              \
              \                    |                                         \
               \                   |         cyc[i].cxlo (c)                  \
                \--------------\   |        /                                  \
                 \              accum_coher()                                   \
                  \            /            \                                    \
                   \      cyc[i].csxy (c)   cyc[i].syy (r) (only if multiple      |
                    \                       /               coherent integration) |
                     \                     /                                     /
                      \                   /                                     /
                       \                 /                                     /
                        \               /    ---------------------------------/
                         \             /    /
                          \----\      /    /
                                accum_norm()
                                     |
                                cyc[i].cnmsxy (c)
                                     |
                                clean_cycspec() - optional filter
                               /              \
                       cyc[i].cnmsxy (c)   cyc[i].mag (r)
                                 \           /
                                  cyc[i] file  (choose one of the above)
 */
public class cyclofunc extends Primitive {

  private static final int DEBUG_BASIC    = 0x01;
  private static final int DEBUG_EXTRA    = 0x02;
  private static final int DEBUG_STATROT  = 0x04;
  private static final int DEBUG_FFT      = 0x08;
  private static final int DEBUG_SXXPSD   = 0x10;
  private static final int DEBUG_COHER    = 0x20;
  private static final int DEBUG_NCOHER   = 0x40;

  private static final int max_ncyclo     = 10;
  private static final int max_fft_size   = 4194394;
  private static final int med_fft_size   = 1024;
  private static final int min_fft_size   = 16;
  private static final int max_xlate_taps = 11;
  private static final float floorval     = 1.0e-15f;
  //  private static final int p_max_modes    = 3; // commented out unused constant in 3.3.3
  private static final double sqrt2o2     = Math.sqrt(2.0)/2.0;
  private static final String rptSubrecordFormat =
      "FILE/XA,NCYC/SD,FREQ/9D,FTSZ/SD,AVGS/SD,OVLP/SD,TTAG/SD," +
      "TDUR/SD,TINC/SD,SIZE/SD,SAMP/SD,MEAN/SD,MMAX/SD,MMIN/SD," +
      "VARI/SD,VMAX/SD,VMIN/SD";

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

  private static final class RealTime {
    public static final int DEFAULT = 0;
    public static final int NOW = 1;
    public static final int SPECIFIED = 2;
  }

  class CycRecord {
    Data syy;                        //auto spectrum
    Data csxy;                       //coherent accumulator
    Data cnmsxy;                     //non-coherent accumulator
    Data mag;                        //Magnitude if they want that
    double      cycfreq;             //Cycle frequency
    double      blkphas;             //Block phase correction
    int         nfft;                //Size of fft
    int         nfft2;               //Half size of fft
    int         ffttype;             //0=Real, 1=Complex
    int         nshift;              //Number of bins to shift
    int         nsmth;               //Number of taps in smoothing function
    int         ntaps;               //Number of active taps in cxlo
    float[]     cxlo;                //Complex translation coeff.
    boolean     isrslt;              //Is cycfrq a result?
    String      nmrslt;              //Result name

    CycRecord () {
      //
      // Reserve room for the translation coefficients.  num_xlate_taps
      // describes the number of coefficients, but the entries are complex
      // so we need twice the space.
      cxlo = new float[2*num_xlate_taps];
    }
  }

  private DataFile hin,hpsd,hrpt;
  private ICEPacket pkt;
  private int debug;
  private int thinfact,thincnt; // output every thinfact frame
  private boolean log_smooth;
  private boolean log_magout;
  private boolean log_normalize;
  private double dminpwr;
  private String c_wind;
  private int num_xlate_taps;     // Number of translation taps
  private float rwhiten, rwhiten_avg;
  private double dsmooth;
  private int izerobins;
  private double dfrqbgn;
  private int ifrqbgn;
  private double dfrqend;
  private int ifrqend;
  private String name_realtime;
  private int mode_realtime;
  private int nfft;
  private double doverlap;
  private boolean log_navg;
  private int navg;               // Number of FFTs to average in PSD
  private String c_navg;
  private int nacc_coher, nacc_coher_max;
  private int nacc_norm, nacc_norm_max;
  private int ncyclo;             // Number of cycle frequencies
  private double dxstart;
  private int isubsize;
  private CycRecord[] Cyc;
  private DataFile[] hCyc;
  private Data grab;              // grab buffer
  private float[] win;            // time domain window
  private Data sxx;               // sxx spectrum estimate
  private Data psd;               // output power spectrum
  private Data wddat;             // reformatted/windowed data
  private double dcons_len;       // nominal consumption length
  private double dfsamp;
  private Table Rpt;              // stub entries for report file
  private Fft fft;                // FFT engine
  String name_navg;
  private int gnum = 0;
  private double ttag = 0.0;
  //  private int recOffset = 0; // comment out unused field in 3.3.3

  private double dst_rsum;
  private double dst_isum;
  private double dst_sqsum;

  private class Statistics {
    public double mean; // mean over accumulation window
    public double mmax; // maximum mean over all FFT windows
    public double mmin; // minimum mean over all FFT windows
    public double vari; // variance over accumulation window
    public double vmax; // maximum variance over all FFT windows
    public double vmin; // minimum variance over all FFT windows
    public int    ntot; // total samples in accumulation window
    public int    nsub; // total samples in FFT window

    Statistics () {
      mean = 0.0;
      mmax = -1.0e20;
      mmin = 1.0e20;
      vari = 0.0;
      vmax = -1.0e20;
      vmin = 1.0e20;
      ntot = 0;
      nsub = nfft;
    }
  }

  Statistics Stat;

  private boolean log_gate;
  private boolean log_init_psd;
  private boolean log_init_sxx;
  private boolean log_init_sxy;
  private boolean log_init_nmsxy;

  public int open() {

    //
    // Command line switches
    //
    debug = MA.getL("/DEBUG");
    // No debug flag specified => we get -1
    if (-1 == debug)
      debug = 0;
    // Force basic if extra is on
    if (0 != (debug & DEBUG_EXTRA))
      debug |= DEBUG_BASIC;

    if (MA.find("/SMOOTH")) {
      log_smooth = true;
    }
    else {
      log_smooth = false;
    }

    if (MA.find("/MAGOUT")) {
      log_magout = true;
    }
    else {
      log_magout = false;
    }

    if (MA.find("/NORM")) {
      log_normalize = true;
    }
    else {
      log_normalize = false;
    }

    dminpwr = MA.getD("/MINPWR", 1.0e-20);

    c_wind = MA.getString ("/WIND","BH61",0);

    num_xlate_taps = MA.getL ("/NTAPS", 5);
    if ((num_xlate_taps % 2) == 0) {
      num_xlate_taps++;
      num_xlate_taps = Math.min(Math.max(3,num_xlate_taps),max_xlate_taps);
    }

    rwhiten = (float) MA.getD("/WHITE",10.0*(double)floorval);
    rwhiten = (float) Math.min(Math.max(0.0,rwhiten),1.0);

    dsmooth = MA.getD ("/SMOOTH", 0.1);
    if (1.0 == dsmooth)
      dsmooth = 0.1;
    dsmooth = Math.min (Math.max (0.0, dsmooth), 1.0);

    izerobins = MA.getL ("/ZEROBINS",3);
    izerobins = Math.min(Math.max(0, izerobins), min_fft_size);

    if (MA.find("/FRQBGN")) {
      dfrqbgn = MA.getD("/FRQBGN");
      ifrqbgn = 0;
    }
    else {
      dfrqbgn = 0.0d;
      ifrqbgn = -1;
    }

    if (MA.find("/FRQEND")) {
      dfrqend = MA.getD ("/FRQEND");
      ifrqend = 0;
    }
    else {
      dfrqend = 0.0d;
      ifrqend = -1;
    }

    if (MA.find("/REALTIME")) {
      name_realtime = MA.getU("/REALTIME");
      if (name_realtime.length() > 0)
        mode_realtime = RealTime.SPECIFIED;
      else
        mode_realtime = RealTime.NOW;
    }
    else {
      mode_realtime = RealTime.DEFAULT;
    }

    thinfact = MA.getL("/THINFACT",1);
    thinfact = Math.max(1,thinfact);
    thincnt = 0;

    //
    // The input file
    //

    hin = MA.getDataFile ("IN","1000|2000","S#|C#|-SP|-CP",0);
    hin.open();
    if (hin.typeClass > 1) hin.setType(1000);
    pkt = (ICEPacket)hin.getPacketHandler();
    if (pkt==null) M.info("Input file is not packetized: "+hin.getTag());

    //
    // Other arguments
    //

    //
    // FFT size
    //
    nfft = MA.getL("NFFT");
    // TODO: Need to ensure it will factor - orig code does this with m$genradix
    // libnmfft has fft_checkradix, but it's C code
    if (nfft <= 1)
      nfft = med_fft_size;
    nfft = Math.min(max_fft_size,Math.max(min_fft_size, nfft));
    //	if (false) {
    //	    nfft = MMath.power2GE (nfft);
    //	    M.warning ("Unable to factor FFT size.  Using " + nfft);
    //	}
    if (0 != (debug & DEBUG_BASIC))
      M.info ("FFT size: " + nfft);

    grab = hin.getDataBuffer (nfft);

    //
    // Time Domain overlap
    //
    doverlap = MA.getD ("OVERLAP", 0.0);
    doverlap = Math.min(Math.max (0.0, doverlap), 0.999999);
    if (0 != (debug & DEBUG_BASIC))
      M.info ("Overlap: " + doverlap);

    //
    // Number of averages
    // It may be a result parameter
    //
    name_navg = MA.getU ("INTEGRATE");
    c_navg = MR.getS (name_navg);
    if (null != c_navg) {
      log_navg = true;
    }
    else {
      c_navg = name_navg;
      log_navg = false;
    }
    String noBars = c_navg.replace ('|', ' ');
    String noQuotes = noBars.replace ('"', ' ');
    nacc_coher_max = 0;
    nacc_norm_max = 0;
    try {
      StringTokenizer tokens = new StringTokenizer (noQuotes);
      nacc_coher_max = Integer.parseInt(tokens.nextToken());
      nacc_norm_max  = Integer.parseInt(tokens.nextToken());
    }
    catch (Exception e) {
      M.warning("Invalid integration parameter");
    }
    nacc_coher_max = Math.min(Math.max(1, nacc_coher_max),2147483647);
    nacc_norm_max  = Math.min(Math.max(1, nacc_norm_max), 2147483647);
    navg = (int) Math.min((double)nacc_coher_max * (double)nacc_norm_max, 2147483647);
    rwhiten_avg = rwhiten * (float)nacc_coher_max * (float)nacc_coher_max;
    if (0 != (debug & DEBUG_EXTRA))
      M.info ("navg,ncoher,nnorm: " + navg + " " + nacc_coher_max +" "+nacc_norm_max);

    //
    // Cycle frequency
    // Might be a value, a result parameter or a list
    //
    ncyclo = 0;
    String cycloParam = MA.getU("CFREQ");
    try {
      StringTokenizer tokens = new StringTokenizer(cycloParam,"|");
      while (tokens.hasMoreTokens()) {
        ncyclo++;
        tokens.nextToken();
      }
    }
    catch (Exception e) {
    }
    if (ncyclo > max_ncyclo) {
      ncyclo = max_ncyclo;
      M.warning ("Too many cycle frequencies specified.");
    }
    Cyc = new CycRecord[ncyclo];
    try {
      StringTokenizer tokens = new StringTokenizer(cycloParam,"|");
      for (int cycloNum = 0; cycloNum < ncyclo; cycloNum++) {
        Cyc[cycloNum] = new CycRecord();
        Cyc[cycloNum].nmrslt = tokens.nextToken();
        if (null == MR.get(Cyc[cycloNum].nmrslt)) {
          Cyc[cycloNum].isrslt = false;
          Cyc[cycloNum].cycfreq = Double.parseDouble(Cyc[cycloNum].nmrslt);
        }
        else {
          Cyc[cycloNum].isrslt = true;
          Cyc[cycloNum].cycfreq = MR.getD (Cyc[cycloNum].nmrslt);
        }
      }
    }
    catch (Exception e) {
    }
    if (0 != (debug & DEBUG_EXTRA)) {
      M.info ("ncyclo: "+ncyclo);
      for (int i = 0; i < ncyclo; i++) {
        M.info (Cyc[i].nmrslt + "\t" + Cyc[i].isrslt + "\t" + Cyc[i].cycfreq);
      }
    }

    // Determine the output subsize and xstart based on the input type
    // Also set up for the FFT based on the data type.
    int fftFlags = Fft.FORWARD;
    if (hin.getFormatMode() == 'C') {
      dxstart = -((double)(nfft/2 - 1)) / (hin.getXDelta()*((double)nfft));
      isubsize = nfft;
      fftFlags |= Fft.COMPLEX;
    }
    else {
      dxstart = 0.0;
      isubsize = nfft/2 + 1;
      fftFlags |= Fft.REAL;
    }
    fft = new Fft();
    fft.init (nfft, fftFlags);

    //
    // Output files
    //

    //
    // Set up the spectra output files
    //

    hCyc = new DataFile[ncyclo];
    for (int i = 0; i < ncyclo; i++) {
      String fileName = MA.getS ("OUT");
      if (ncyclo > 1) {
        // Build output file names suffixed by _##, where ## is the two digit
        // spectrum number.  An alternate approach would be NxM's PrintfFormat class.
        NumberFormat f = NumberFormat.getInstance();
        f.setMinimumIntegerDigits(2);
        fileName = fileName + "_" + f.format (i+1);
      }
      if (log_magout)
        hCyc[i] = new DataFile (M, fileName, "2000",
            "SF", BaseFile.OUTPUT|BaseFile.OPTIONAL, 1.0, null);
      else
        hCyc[i] = new DataFile (M, fileName, "2000",
            "CF", BaseFile.OUTPUT|BaseFile.OPTIONAL, 1.0, null);

      hCyc[i].setXDelta (1.0/(hin.getXDelta()*(double)nfft));
      hCyc[i].setXStart (dxstart);
      hCyc[i].setXUnits (3);

      hCyc[i].setYDelta (hin.getXDelta()*(1.0-doverlap)*((double)(nfft*navg)));
      hCyc[i].setYStart (hin.getXStart() + 0.5*hCyc[i].getYDelta());
      hCyc[i].setYUnits (1);

      hCyc[i].setSubSize (isubsize);
      hCyc[i].setTimeCode (hin.getTimeCode ());
      hCyc[i].open();
    }

    // PSD file/pipe

    hpsd = MA.getDataFile ("PSDOUT", "2000", "SF", BaseFile.OUTPUT|BaseFile.OPTIONAL,
        1.0, null);
    hpsd.setXDelta (1.0/(hin.getXDelta()*(double)nfft));
    hpsd.setXStart (dxstart);
    hpsd.setXUnits (3);

    hpsd.setYDelta (hin.getXDelta() * (1.0 - doverlap) * (double)(nfft * navg));
    hpsd.setYStart (hin.getXStart() + 0.5 * hpsd.getYDelta ());
    hpsd.setYUnits (1);

    hpsd.setSubSize (isubsize);
    hpsd.setTimeCode (hin.getTimeCode ());
    hpsd.open();

    // RPT file/pipe

    hrpt = MA.getDataFile ("RPTOUT", null, BaseFile.OUTPUT|BaseFile.OPTIONAL);

    hrpt.setXDelta (hin.getXDelta() * (1.0 - doverlap) * (double)(nfft * navg));
    // TODO: Should't this be 0.5*hrpt.getXDelta() ?  The FORTRAN says ydelta, but
    // at this point ydelta will be zero.
    hrpt.setXStart (hin.getXStart() + 0.5*hrpt.getYDelta());
    hrpt.setType (3000);
    hrpt.setSubRecords (rptSubrecordFormat);
    hrpt.setTimeCode (hin.getTimeCode ());
    hrpt.open();

    //
    // Set up buffers
    //

    sxx   = new Data ("SF", nfft);
    psd   = hpsd.getDataBuffer (hpsd.getSubSize());
    if (hin.getFormatMode() == 'C')
      wddat = new Data ("CF", nfft);
    else
      wddat = new Data ("SF", nfft+2);

    for (int i = 0; i < ncyclo; i++) {
      Cyc[i].cnmsxy = new Data ("CF", hCyc[i].getSubSize());
      Cyc[i].mag    = new Data ("SF", hCyc[i].getSubSize());
      Cyc[i].csxy   = new Data ("CF", hCyc[i].getSubSize());
      Cyc[i].syy    = new Data ("SF", hCyc[i].getSubSize());
    }

    //
    // Create a window for time domain data
    //
    win = Window.get (c_wind, nfft, 1.0/((float)(nfft * navg)));

    //
    // Compute frequency start/stop bins for cyclic spectra
    //
    if (ifrqbgn == -1) dfrqbgn = hpsd.getXStart();
    if (ifrqend == -1) dfrqend = hpsd.getXStart() + ((double)(hpsd.getFrameSize()-1))*hpsd.getXDelta();
    dfrqbgn = (dfrqbgn-hpsd.getXStart())/Math.max(1.0e-20,hpsd.getXDelta());
    dfrqend = (dfrqend-hpsd.getXStart())/Math.max(1.0e-20,hpsd.getXDelta());
    ifrqbgn = (int) (Math.min(Math.max(0.0,dfrqbgn),2147483647.0));
    ifrqend = (int) (Math.min(Math.max(0.0,dfrqend),2147483647.0));
    ifrqbgn = Math.min(Math.max(0,ifrqbgn),hpsd.getFrameSize()-1);
    ifrqend = Math.min(Math.max(0,ifrqend),hpsd.getFrameSize()-1);
    ifrqbgn = Math.min(ifrqbgn,ifrqend);

    if (0 != (debug & DEBUG_EXTRA))
      M.info ("Freq Start/Stop Bins: " + ifrqbgn + "\t" + ifrqend);

    //
    // Compute dataflow transfer stuff
    //

    hin.setXferLength(nfft);

    if (doverlap > 0.0) {
      hin.setConsLength((int) ((1.0-doverlap)*(double)(nfft)));
      dcons_len      = (double)hin.getConsLength();
    }
    else {
      hin.setConsLength(-1);
      dcons_len      =  (int) ((1.0-doverlap)*(double)(nfft));
    }

    //
    // Compute the bin offset and translation coefficients
    //
    for (int i = 0; i < ncyclo; i++) {
      Cyc[i].ntaps = num_xlate_taps;
      Cyc[i].nshift = gen_xlate_taps (Cyc[i].cxlo,
          Cyc[i].ntaps,
          Cyc[i].cycfreq,
          hin.getXDelta());
    }

    //
    // Compute smoothing window sizes (assume at least 20% excess bandwidth)
    //
    dfsamp = 1.0/hin.getXDelta();;
    for (int i = 0; i < ncyclo; i++) {
      Cyc[i].nsmth   = (int) ((dsmooth*Cyc[i].cycfreq)/
          (dfsamp/(double)nfft));
      Cyc[i].nsmth   = Math.min(Math.max(1,Cyc[i].nsmth), Cyc[i].nshift);
      Cyc[i].blkphas = 0.0;
      if (0 != (debug & DEBUG_BASIC))
        M.info ("Bin offset: " + Cyc[i].nshift + "\t" + Cyc[i].nsmth);
    }

    //
    // Initialize report structure.  Build the report structure from the
    // rptSubrecordFormat string so it will match the structure of the
    // report file created with that format string.
    //
    Rpt = new Table(Table.KV);

    try {
      StringTokenizer tokens = new StringTokenizer(rptSubrecordFormat,",");
      while (tokens.hasMoreTokens()) {
        StringTokenizer key = new StringTokenizer (tokens.nextToken(), "/");
        Rpt.setKey (key.nextToken (), key.nextToken ());
      }
    }
    catch (Exception e) {
      M.info ("Error parsing " + rptSubrecordFormat);
    }
    Rpt.put ("FILE", (hin.getURL().length() < 80) ? hin.getURL() : hin.getURL().substring (0, 80));
    Rpt.put ("NCYC", ncyclo);
    //double[] frq = new double[9];
    Rpt.put ("FREQ", 0);
    Rpt.put ("FTSZ", (double) nfft);
    Rpt.put ("AVGS", (double) navg);
    Rpt.put ("OVLP", doverlap);
    Rpt.put ("TTAG", 0.0);
    ttag = 0.0;
    Rpt.put ("TDUR", 0.0);
    Rpt.put ("TINC", 0.0);
    Rpt.put ("SIZE", 0.0);
    Rpt.put ("SAMP", 0.0);
    Rpt.put ("MEAN", 0.0);
    Rpt.put ("MMAX", 0.0);
    Rpt.put ("MMIN", 0.0);
    Rpt.put ("VARI", 0.0);
    Rpt.put ("VMAX", 0.0);
    Rpt.put ("VMIN", 0.0);

    if (mode_realtime == RealTime.SPECIFIED) {
      Rpt.put ("SIZE",-1.0);
      if (null == MR.get (name_realtime))
        M.warning ("Time result not found: "+name_realtime);
      else
        Rpt.put("TTAG", MR.getD (name_realtime));
    }
    else if (mode_realtime == RealTime.NOW) {
      Rpt.put ("SIZE", -1.0);
      Rpt.put ("TTAG", Time.current());
    }
    else {
      Rpt.put ("SIZE", hin.getSize());
      ttag = hin.getXStart() +
          (hin.getOffset() -(1.0-doverlap) * (double) (nfft*navg)) * hin.getXDelta();
      Rpt.put ("TTAG", ttag);
    }
    Rpt.put ("TDUR", hpsd.getYDelta());
    Rpt.put ("TINC", hin.getXDelta());
    Rpt.put ("SAMP", hin.getOffset());

    //
    // Initialize statistics structure
    //
    Stat = new Statistics();

    dst_rsum = 0.0;
    dst_isum = 0.0;
    dst_sqsum = 0.0;

    //
    // Other variable initialization
    //
    log_gate       = false;
    log_init_psd   = true;
    log_init_sxx   = true;
    log_init_sxy   = true;
    log_init_nmsxy = true;

    nacc_coher = 0;
    nacc_norm  = 0;

    return NORMAL;
  }

  public int process() {

    int n = hin.read(grab);
    thincnt++;
    if (thincnt<thinfact) return NOOP;
    thincnt=0;

    gnum = gnum + 1;

    if (0 != (debug & (DEBUG_STATROT|DEBUG_FFT|DEBUG_SXXPSD|DEBUG_COHER|DEBUG_NCOHER))) {
      M.info ("**************************************************************");
      M.info ("grab "+gnum+" got "+n);
    }

    if (n>0) {

      // Window (zero-pad if required) and take statistics
      // from the input time data

      int istat = win_stat_rot (win, grab, wddat, n);

      if (0 != (debug & DEBUG_STATROT)) {
        M.info ("wddat after win_stat_rot");
        float[] x = wddat.castF(true);
        if (hin.getFormatMode() == 'C') {
          for (int i = 0; i < 10; i++)
            M.info ("  "+i+", ("+x[2*i]+","+x[2*i+1]+")");
          if (n < nfft) {
            for (int i = n-5; i < n+5; i++)
              M.info ("  "+i+", ("+x[2*i]+","+x[2*i+1]+")");
          }
          for (int i = nfft-10; i < nfft; i++)
            M.info ("  "+i+", ("+x[2*i]+","+x[2*i+1]+")");
        }
        else {
          for (int i = 0; i < 10; i++)
            M.info ("  "+i+", "+x[i]);
          if (n < nfft) {
            for (int i = n-5; i < n+5; i++)
              M.info ("  "+i+", "+x[i]);
          }
          for (int i = nfft-10; i<nfft; i++)
            M.info ("  "+i+", "+x[i]);
        }
        wddat.uncast (x, true);
      }

      // Update the time domain statistics
      Stat.mean += dst_rsum;
      Stat.vari += dst_sqsum;
      Stat.ntot += nfft;
      Stat.nsub  = nfft;

      if (0 != (debug & DEBUG_STATROT)) {
        M.info ("Stat.mean "+Stat.mean);
        M.info ("Stat.vari "+Stat.vari);
        M.info ("Stat.ntot "+Stat.ntot);
        M.info ("Stat.nsub "+Stat.nsub);
      }

      if (dst_rsum > Stat.mmax)  Stat.mmax = dst_rsum;
      if (dst_rsum < Stat.mmin)  Stat.mmin = dst_rsum;
      if (dst_sqsum > Stat.vmax) Stat.vmax = dst_sqsum;
      if (dst_sqsum < Stat.vmin) Stat.vmin = dst_sqsum;

      if (istat < 0)
        M.warning("Unknown input data type: "+hin.getFormat());

      // Compute gating logical
      if (dst_sqsum > ((double) nfft)*dminpwr)
        log_gate = false;
      else
        log_gate = true;

      // FFT the data (rotate complex by 1 bin)

      if (!log_gate) {
        fft.work (wddat);
        if (hin.getFormatMode() == 'C') {
          // would like to:
          //fft.rotate (wddat);
          // but that doesn't appear to do what we want
          // manually rotate instead
          float b0r = wddat.getF(0);
          float b0i = wddat.getF(1);
          DataOp.copy (wddat, 1, wddat, 0, wddat.getSize()-1);
          wddat.setF(2*(wddat.getSize()-1), b0r);
          wddat.setF(2*(wddat.getSize()-1)+1, b0i);
        }
        else {
          // Unpack Nyquist bin
          // TODO - This will be more efficient with individual access rather than
          // the full cast - look into that
          float[] x = wddat.castF(true);
          x[nfft] = x[1];
          x[1] = 0.0f;
          wddat.uncast (x, true);
        }

        if (0 != (debug & DEBUG_FFT)) {
          M.info ("after fft");
          float[] x = wddat.castF(true);
          for (int j = 0; j < 10; j++)
            System.out.println (j+"\t("+x[2*j]+","+x[2*j+1]+")");
          for (int j = nfft-10; j < nfft; j++)
            System.out.println (j+"\t("+x[2*j]+","+x[2*j+1]+")");
          wddat.uncast(x,true);
        }

        // Compute Sxx/PSD update
        float ftemp = 1.0f;
        if (log_normalize)
          ftemp = (float)
          (1.0 / Math.sqrt(Math.max(1.0e-20,dst_sqsum/((double)Stat.nsub))));

        // Do SXX and Power spectrum updates
        istat = accum_sxx_psd (wddat, sxx,
            psd, ftemp, hpsd.getSubSize(),
            izerobins);
        if (0 != (debug & DEBUG_SXXPSD)) {
          M.info ("after accum_sxx_psd");
          M.info ("sxx");
          float[] x = sxx.castF (true);
          for (int i = 0; i < 20; i++)
            M.info ("\t"+i+"\t"+x[i]);
          for (int i = hpsd.getSubSize()-10; i < hpsd.getSubSize(); i++)
            M.info ("\t"+i+"\t"+x[i]);
          sxx.uncast (x, true);
          M.info ("psd");
          x = psd.castF (true);
          for (int i = 0; i < 10; i++)
            M.info ("\t"+i+"\t"+x[i]);
          for (int i = hpsd.getSubSize()-10; i < hpsd.getSubSize(); i++)
            M.info ("\t"+i+"\t"+x[i]);
          psd.uncast (x, true);
        }
      }
      // Update the Coherent and Normalized Accumulators
      nacc_coher++;

      for (int i = 0; i < ncyclo; i++) {
        if (0 != (debug & (DEBUG_COHER | DEBUG_NCOHER))) {
          M.info ("cycle\t"+i);
        }
        Cyc[i].blkphas = (Cyc[i].blkphas -
            Cyc[i].cycfreq*hin.getXDelta()*dcons_len) % 1.0;
        if (!log_gate) {
          float[] cxphase = new float[2];
          cxphase[0] = (float) Math.cos(Cyc[i].blkphas*TWOPI);
          cxphase[1] = (float) Math.sin(Cyc[i].blkphas*TWOPI);

          if (0 != (debug & DEBUG_COHER)) {
            M.info (" blkphas ("+cxphase[0]+","+cxphase[1]+")");
          }
          istat = accum_coher (wddat,
              Cyc[i].csxy,
              Cyc[i].syy,
              Cyc[i],
              cxphase,
              hCyc[i].getSubSize(),
              nacc_coher_max,
              ifrqbgn,
              ifrqend);
          if (i >= (ncyclo-1)) log_init_sxy = false;
        }

        if (nacc_coher >= nacc_coher_max) {
          if (!log_gate) {
            istat = accum_norm (Cyc[i].csxy,
                sxx,
                Cyc[i].syy,
                Cyc[i].cnmsxy,
                Cyc[i],
                hCyc[i].getSubSize(),
                nacc_coher_max,
                nacc_norm_max,
                rwhiten_avg,
                ifrqbgn,
                ifrqend);
            if (i >= (ncyclo-1)) log_init_nmsxy = false;
          }
          if (i >= (ncyclo-1)) {
            log_init_sxx   = true;
            log_init_sxy   = true;
            nacc_coher = 0;
            nacc_norm++;
          }
        }
      }

      // Compute the output if navg iterations
      if (nacc_norm >= nacc_norm_max) {
        for (int i = 0; i < ncyclo; i++) {
          // Ensure vector was initialized
          if (log_init_nmsxy) {
            DataOp.zero (Cyc[i].cnmsxy, 0, Cyc[i].cnmsxy.getSize());
          }
          else {
            // Smooth cyclic spectrum if requested
            if (log_smooth || log_magout)
              istat = clean_cycspec (Cyc[i].cnmsxy,
                  Cyc[i].mag,
                  hCyc[i].getSubSize(),
                  ifrqbgn,
                  ifrqend,
                  Cyc[i].nsmth,
                  log_smooth,
                  log_magout);
          }
          // Output cyclic spectra
          if (hCyc[i].getFormat().startsWith("C"))
            hCyc[i].write (Cyc[i].cnmsxy, 1);
          else
            hCyc[i].write (Cyc[i].mag, 1);
        }

        // Output power spectra
        if (log_init_psd) {
          double ftemp = dminpwr/(Math.max(1,navg*nfft));
          ftemp = 0.2021311785*ftemp*ftemp;
          DataOp.fill (psd, 0, ftemp, hpsd.getSubSize());
        }

        if (izerobins > 0) {
          DataOp.fill (psd, 0,
              psd.getF(izerobins), izerobins);
          DataOp.fill (psd, hpsd.getSubSize()-izerobins,
              psd.getF(hpsd.getSubSize()-izerobins-1), izerobins);
        }
        hpsd.write (psd, 1);

        // Build a report record
        if (hrpt.isOpen()) {
          double mean = Stat.mean/Math.max(1.0,(double) Stat.ntot);
          Rpt.put ("MEAN", mean);
          Rpt.put ("VARI", Stat.vari - mean*mean);
          Rpt.put ("MMAX", Stat.mmax/Math.max(1.0,(double)Stat.nsub));
          Rpt.put ("mmin", Stat.mmin/Math.max(1.0,(double)Stat.nsub));
          Rpt.put ("vmax", Stat.vmax/Math.max(1.0,(double)Stat.nsub));
          Rpt.put ("vmin", Stat.vmin/Math.max(1.0,(double)Stat.nsub));
          for (int i = 0; i < ncyclo; i++) {
            String f1 = new String ("FREQ(");
            String f2 = f1.concat(String.valueOf(i));
            String f3 = f2.concat(")");
            Rpt.put (f3, Cyc[i].cycfreq);
          }
          Rpt.put ("AVGS", navg);
          Rpt.put ("SAMP", hin.getOffset());
          switch (mode_realtime) {
            case RealTime.SPECIFIED:
              Rpt.put ("SIZE", -1.0);
              Rpt.put ("TTAG",  0.0);
              //ch = M$RFIND(name_realtime,aaval,RPT.ttag)
              //if (ch(1:1) .eq. ' ')
              //	         call M$WARNING('Time result not found: '//name_realtime)
              break;

            case RealTime.NOW:
              Rpt.put ("SIZE", -1.0);
              ttag += Time.current ();
              Rpt.put ("TTAG", ttag);
              break;

            case RealTime.DEFAULT:
              Rpt.put ("SIZE", hin.getSize());
              ttag = hin.getXStart() +
                  (hin.getOffset()-(1.0-doverlap)*(double)(nfft*navg))*hin.getXDelta();
              Rpt.put ("TTAG", ttag);
          }
          hrpt.setData (hrpt.getOffset(), Rpt);
          Stat.mean =   0.0;
          Stat.mmax = -1.0e20;
          Stat.mmin =  1.0e20;
          Stat.vari =   0.0;
          Stat.vmax = -1.0e20;
          Stat.vmin =  1.0e20;
          Stat.ntot = 0;
          Stat.nsub = nfft;
        }

        // Check for change in cycle frequency
        // Doesn't look like we can check Mc.break anymore, even with Mc.getBreak()
        for (int i = 0; i < ncyclo; i++) {

          if (Cyc[i].isrslt) {
            double ddval = MR.getD (Cyc[i].nmrslt);

            if (Math.abs(Cyc[i].cycfreq-ddval) >= 1.0e-10) {
              Cyc[i].cycfreq = ddval;
              Cyc[i].nshift = gen_xlate_taps (Cyc[i].cxlo,
                  Cyc[i].ntaps,
                  Cyc[i].cycfreq,
                  hin.getXDelta());
              Cyc[i].blkphas = 0.0;
              Cyc[i].nsmth = (int) ((dsmooth*Cyc[i].cycfreq)/
                  (dfsamp/(double)nfft));
              Cyc[i].nsmth   = Math.min(Math.max(1,Cyc[i].nsmth), Cyc[i].nshift);
              DataOp.fill (Cyc[i].cnmsxy, 0, 0.0, Cyc[i].cnmsxy.getSize());
              DataOp.fill (Cyc[i].csxy,   0, 0.0, Cyc[i].csxy.getSize());
              DataOp.fill (Cyc[i].syy,    0, 0.0, Cyc[i].syy.getSize());
            }
          }
        }

        // Check for change in coherent/normalized averages
        if (log_navg) {
          /*
          ddval = 0.0d0
          aaval = ' '
          ch    = ' '
          ch    = M$RFIND(name_navg(1:len_navg),aaval,ddval)
           */
          String aaval = MR.getS (name_navg);
          if (0 != aaval.compareTo (c_navg)) {
            c_navg = aaval;
            String noBars = c_navg.replace ('|', ' ');
            String noQuotes = noBars.replace ('"', ' ');
            nacc_coher_max = 0;
            nacc_norm_max = 0;
            try {
              StringTokenizer tokens = new StringTokenizer (noQuotes);
              nacc_coher_max = Integer.parseInt(tokens.nextToken());
              nacc_norm_max  = Integer.parseInt(tokens.nextToken());
            }
            catch (Exception e) {
              M.warning("Invalid integration parameter");
            }
            nacc_coher_max = Math.min(Math.max(1, nacc_coher_max),2147483647);
            nacc_norm_max  = Math.min(Math.max(1, nacc_norm_max), 2147483647);
            navg = (int) Math.min((double)nacc_coher_max
                * (double)nacc_norm_max, 2147483647);
            rwhiten_avg = rwhiten * (float)nacc_coher_max * (float)nacc_coher_max;

            win = Window.get (c_wind, nfft, 1.0/((float)(nfft * navg)));
          }
        }
        // Reset counters/logic states
        log_init_psd   = true;
        log_init_sxx   = true;
        log_init_sxy   = true;
        log_init_nmsxy = true;
        nacc_coher = 0;
        nacc_norm  = 0;
      }
    }
    else if (n==0)
      return NOOP;
    else
      return FINISH;
    return NORMAL;
  }

  public int close() {

    hin.close();
    for (int i = 0; i < ncyclo; i++)
      hCyc[i].close ();
    hpsd.close();
    hrpt.close();

    return NORMAL;
  }

  /**
   ************************************************************************
   * FUNCTION WIN_STAT_ROT
   ************************************************************************
   *	WIN_STAT_ROT --- Applies a time-domain window to data set;
   *                     computes sum, sumsquared;
   *	                 does complex data rotation for FFT processing
   */
  int win_stat_rot (float[] awin,    // in  - Window data
                    Data dadat,      // in  - data
                    Data daout,      // out - windowed data
                    int n) {         // number of total elements

    int retVal = n;

    if (0 != (debug & DEBUG_STATROT))
      M.info ("win_stat_rot nelem "+n+", n "+daout.getSize());

    float[] asfout = daout.castF(true);
    // Initialize statistical stuff
    dst_rsum  = 0.0;
    dst_isum  = 0.0;
    dst_sqsum = 0.0;

    String cfmt = dadat.getFormat();

    if (cfmt.startsWith("SF")) {
      float[] sfbuf = dadat.castF(true);
      for (int i = 0; i < n; i++) {
        dst_rsum  += sfbuf[i];
        dst_sqsum += sfbuf[i]*sfbuf[i];
        asfout[i] = awin[i]*sfbuf[i];
      }
      dadat.uncast(sfbuf,true);
    }
    else if (cfmt.startsWith("SB")) {
      byte[] sbbuf = dadat.castB(true);
      for (int i = 0; i < n; i++) {
        asfout[i] = sbbuf[i];
        dst_rsum  += asfout[i];
        dst_sqsum += asfout[i]*asfout[i];
        asfout[i] = awin[i]*asfout[i];
      }
      dadat.uncast(sbbuf, true);
    }
    else if (cfmt.startsWith("SI")) {
      short[] sibuf = dadat.castI(true);
      for (int i = 0; i < n; i++) {
        asfout[i] = sibuf[i];
        dst_rsum  += asfout[i];
        dst_sqsum += asfout[i]*asfout[i];
        asfout[i] = awin[i]*asfout[i];
      }
      dadat.uncast (sibuf, true);
    }
    else if (cfmt.startsWith("SL")) {
      int[] slbuf = dadat.castL(true);
      for (int i = 0; i < n; i++) {
        asfout[i] = slbuf[i];
        dst_rsum  += asfout[i];
        dst_sqsum += asfout[i]*asfout[i];
        asfout[i] = awin[i]*asfout[i];
      }
      dadat.uncast(slbuf, true);
    }
    else if (cfmt.startsWith("SD")) {
      double[] sdbuf = dadat.castD(true);
      for (int i = 0; i < n; i++) {
        asfout[i] = (float) sdbuf[i];
        dst_rsum  += asfout[i];
        dst_sqsum += asfout[i]*asfout[i];
        asfout[i] = awin[i]*asfout[i];
      }
      dadat.uncast(sdbuf, true);
    }
    else if (cfmt.startsWith("CB")) {
      byte[] sbbuf = dadat.castB(true);
      for (int i = 0; i < n; i+=2) {
        int ci = i*2;
        asfout[ci]   = sbbuf[ci];
        asfout[ci+1] = sbbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = awin[i]*asfout[ci];
        asfout[ci+1] = awin[i]*asfout[ci+1];

        ci += 2;
        asfout[ci]   = sbbuf[ci];
        asfout[ci+1] = sbbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = -awin[i+1]*asfout[ci];
        asfout[ci+1] = -awin[i+1]*asfout[ci+1];
      }
      dadat.uncast(sbbuf, true);
    }
    else if (cfmt.startsWith("CI")) {
      short[] sibuf = dadat.castI(true);
      for (int i = 0; i < n; i+=2) {
        int ci = i*2;
        asfout[ci]   = sibuf[ci];
        asfout[ci+1] = sibuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = awin[i]*asfout[ci];
        asfout[ci+1] = awin[i]*asfout[ci+1];

        ci += 2;
        asfout[ci]   = sibuf[ci];
        asfout[ci+1] = sibuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = -awin[i+1]*asfout[ci];
        asfout[ci+1] = -awin[i+1]*asfout[ci+1];
      }
      dadat.uncast(sibuf, true);
    }
    else if (cfmt.startsWith("CL")) {
      int[] slbuf = dadat.castL(true);
      for (int i = 0; i < n; i+=2) {
        int ci = i*2;
        asfout[ci]   = slbuf[ci];
        asfout[ci+1] = slbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = awin[i]*asfout[ci];
        asfout[ci+1] = awin[i]*asfout[ci+1];

        ci += 2;
        asfout[ci]   = slbuf[ci];
        asfout[ci+1] = slbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = -awin[i+1]*asfout[ci];
        asfout[ci+1] = -awin[i+1]*asfout[ci+1];
      }
      dadat.uncast(slbuf, true);
    }
    else if (cfmt.startsWith("CF")) {
      float[] sfbuf = dadat.castF(true);
      for (int i = 0; i < n; i+=2) {
        int ci = i*2;
        asfout[ci]   = sfbuf[ci];
        asfout[ci+1] = sfbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = awin[i]*asfout[ci];
        asfout[ci+1] = awin[i]*asfout[ci+1];

        ci += 2;
        asfout[ci]   = sfbuf[ci];
        asfout[ci+1] = sfbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = -awin[i+1]*asfout[ci];
        asfout[ci+1] = -awin[i+1]*asfout[ci+1];
      }
      dadat.uncast(sfbuf, true);
    }
    else if (cfmt.startsWith("CD")) {
      double[] sdbuf = dadat.castD(true);
      for (int i = 0; i < n; i+=2) {
        int ci = i*2;
        asfout[ci]   = (float) sdbuf[ci];
        asfout[ci+1] = (float) sdbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = awin[i]*asfout[ci];
        asfout[ci+1] = awin[i]*asfout[ci+1];

        ci += 2;
        asfout[ci]   = (float) sdbuf[ci];
        asfout[ci+1] = (float) sdbuf[ci+1];
        dst_rsum  += asfout[ci];
        dst_isum  += asfout[ci+1];
        dst_sqsum += asfout[ci]*asfout[ci]
                   + asfout[ci+1]*asfout[ci+1];
        asfout[ci]   = -awin[i+1]*asfout[ci];
        asfout[ci+1] = -awin[i+1]*asfout[ci+1];
      }
      dadat.uncast(sdbuf, true);
    }
    else {
      if (cfmt.startsWith("C"))
        for (int i = 0; i < 2*n; i++)
          asfout[i] = 0.0f;
      else
        for (int i = 0; i < n; i++)
          asfout[i] = 0.0f;
      retVal = -1;
    }

    int size = daout.getSize();
    if (n < size) {
      if (cfmt.startsWith("C")) {
        for (int i = n*2; i < 2*size; i++)
          asfout[i] = 0.0f;
      }
      else {
        for (int i = n; i < size; i++)
          asfout[i] = 0.0f;
      }
    }

    daout.uncast(asfout, true);
    return retVal;
  }

  /**
   ***********************************************************************
   * FUNCTION GEN_XLATE_TAPS
   ***********************************************************************
     GEN_XLATE_TAPS --- generates the frequency-domain translation taps
   */
  int gen_xlate_taps (float[] cxresp,   // Cyc[i].cxlo
                      int ntaps,        // Cyc[i].ntaps
                      double dfreq,     // Cyc[i].cycfreq
                      double dtincr) {

    Data cxbuf = new Data ("CF", nfft);
    float[] cxbuf1 = cxbuf.castF(true);

    // Create the complex LO response
    double dstrt = -0.5*(nfft-1);
    double darg  = TWOPI * Math.abs(dfreq) * dtincr;

    for (int i = 0; i < nfft; i++) {
      double ang = darg*(dstrt+(double)i);
      cxbuf1[2*i]     = (float) Math.cos(ang);
      cxbuf1[2*i + 1] = (float) Math.sin(ang);
    }

    // Perform complex fft on LO
    cxbuf.uncast (cxbuf1, true);
    Fft fft = new Fft (nfft, Fft.COMPLEX | Fft.FORWARD);
    fft.work (cxbuf);
    cxbuf1 = cxbuf.castF(true);

    //	Find maximum bin in spectrum
    float fmax = 0.0f;
    float ftest;
    int nmax = 0;
    for (int i = 0; i < nfft/2; i++) {
      //M.info (i+"\t("+cxbuf1[2*i]+","+cxbuf1[2*i+1]+")");
      ftest = (cxbuf1[2*i]*cxbuf1[2*i] + cxbuf1[2*i+1]*cxbuf1[2*i+1]);
      if (ftest >= fmax) {
        fmax = ftest;
        nmax = i;
      }
    }

    if (0 != (debug & DEBUG_EXTRA))
      M.info ("nmax,fmax: " + nmax + "\t" + fmax);

    // Extract the coefficients from unwindowed data using peak descent
    // (invert the order for convolution)
    int nctap = ntaps/2; // ntaps is always odd
    for (int i = 0; i < nctap; i++) {
      int ilow  = nmax - i;
      int ihigh = nmax + i;

      // wrap around
      if (ilow  < 0)     ilow  += nfft;
      if (ihigh >= nfft) ihigh -= nfft;

      cxresp[2*(nctap+i)  ] = cxbuf1[2*ilow  ];
      cxresp[2*(nctap+i)+1] = cxbuf1[2*ilow+1];

      cxresp[2*(nctap-i)  ] = cxbuf1[2*ihigh  ];
      cxresp[2*(nctap-i)+1] = cxbuf1[2*ihigh+1];
    }

    if (0 != (debug & DEBUG_EXTRA)) {
      M.info("key tap: 0 ("            +cxbuf1[0]          +","+ cxbuf1[1]             +")");
      M.info("key tap: 1 ("            +cxbuf1[2]          +","+ cxbuf1[3]             +")");
      M.info("key tap: "+(nfft-2) +" ("+cxbuf1[2*(nfft-2)] +","+ cxbuf1[2*(nfft-2) + 1]+")");
      M.info("key tap: "+(nfft-1) +" ("+cxbuf1[2*(nfft-1)] +","+ cxbuf1[2*(nfft-1) + 1]+")");
      for (int i = 0; i < ntaps; i++) {
        M.info ("i,ntaps,ctap: "+i+"\t"+ntaps+"\t("+cxresp[2*i]+","+cxresp[2*i+1]+")");
      }
    }

    // Scale and conjugate the taps
    fmax = 1.0f/(float)nfft;
    for (int i = 0; i < ntaps; i++) {
      cxresp[2*i  ] =   fmax * cxresp[2*i  ];
      cxresp[2*i+1] = - fmax * cxresp[2*i+1];
    }

    cxbuf.uncast (cxbuf1, true);

    return nmax;
  }

  /**
   ************************************************************************
   * FUNCTION ACCUM_SXX_PSD
   ************************************************************************
   *	ACCUM_SXX_PSD --- Accumulator for Sxx and PSD spectra
   */
  int accum_sxx_psd (Data dcfft,       // Input FFT'd data
      Data dfsxx,       // out - Sxx spectrum
      Data dfpsd,       // out - Power spectrum
      float fscale,     // FFT normalization/scale factor
      int n,            // number of input bins
      int iedge) {       // edge bins to zero out

    float[] fpsd = dfpsd.castF(true);
    float[] fsxx = dfsxx.castF(true);
    float[] cfft = dcfft.castF(true);

    if (0 != (debug & DEBUG_SXXPSD)) {
      M.info ("cfft "+dcfft.getSize()+", "+dcfft.getFormat());
      M.info ("fsxx "+dfsxx.getSize()+", "+dfsxx.getFormat());
      M.info ("fpsd "+dfpsd.getSize()+", "+dfpsd.getFormat());
      M.info ("n "+n);
      M.info ("fscale "+fscale);
    }

    if (fscale != 1.0f) {
      float fscalesq = fscale*fscale;
      if ((!log_init_sxx) && (!log_init_psd)) {
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = fpsd[i] + ftemp;
          fsxx[i] = fsxx[i] + fscalesq*ftemp;
          cfft[2*i  ] = fscale*cfft[2*i  ];
          cfft[2*i+1] = fscale*cfft[2*i+1];
        }
      }
      else if ((log_init_sxx) && (!log_init_psd)) {
        log_init_sxx = false;
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = fpsd[i] + ftemp;
          fsxx[i] = fscalesq * ftemp;
          cfft[2*i  ] = fscale*cfft[2*i  ];
          cfft[2*i+1] = fscale*cfft[2*i+1];
        }
      }
      else if ((!log_init_sxx) && (log_init_psd)) {
        log_init_psd = false;
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = ftemp;
          fsxx[i] = fsxx[i] + fscalesq*ftemp;
          cfft[2*i  ] = fscale*cfft[2*i  ];
          cfft[2*i+1] = fscale*cfft[2*i+1];
        }
      }
      else if ((log_init_sxx) && (log_init_psd)) {
        log_init_sxx = false;
        log_init_psd = false;
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = ftemp;
          fsxx[i] = fscalesq * ftemp;
          cfft[2*i  ] = fscale*cfft[2*i  ];
          cfft[2*i+1] = fscale*cfft[2*i+1];
        }
      }
    }
    else {
      if ((!log_init_sxx) && (!log_init_psd)) {
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = fpsd[i] + ftemp;
          fsxx[i] = fsxx[i] + ftemp;
        }
      }
      else if ((log_init_sxx) && (!log_init_psd)) {
        log_init_sxx = false;
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = fpsd[i] + ftemp;
          fsxx[i] = ftemp;
        }
      }
      else if ((!log_init_sxx) && (log_init_psd)) {
        log_init_psd = false;
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = ftemp;
          fsxx[i] = fsxx[i] + ftemp;
        }
      }
      else if ((log_init_sxx) && (log_init_psd)) {
        log_init_sxx = false;
        log_init_psd = false;
        for (int i = 0; i < n; i++) {
          float ftemp   = (cfft[2*i]*cfft[2*i] + cfft[2*i+1]*cfft[2*i+1]);
          fpsd[i] = ftemp;
          fsxx[i] = ftemp;
        }
      }
    }
    dfpsd.uncast (fpsd, true);
    dfsxx.uncast (fsxx, true);
    dcfft.uncast (cfft, true);

    if (iedge > 0) {
      for (int i = 0; i < iedge; i++) {
        fsxx[i    ] = 0.0f;
        cfft[2*i  ] = 0.0f;
        cfft[2*i+1] = 0.0f;

        fsxx[n-1-i      ] = 0.0f;
        cfft[2*(n-1-i)  ] = 0.0f;
        cfft[2*(n-1-i)+1] = 0.0f;
      }
    }

    return n;
  }


  /**
   ************************************************************************
   * FUNCTION ACCUM_COHER
   ************************************************************************
   *	ACCUM_COHER --- Computes cyclic spectrum and performs spectral
   *			averaging
   *                    Convol
   */
  int accum_coher (Data dcsxx,       // in  - Complex Sxx spectrum
                   Data dcsxy,       // out - Complex Sxy spectrum
                   Data dfsyy,       // out - Auto spectrum of shifted signal
                   CycRecord Cinf,   // CYC record information
                   float[] crotat,   // Block rotation factor
                   int n,            // Number of points in input/output
                   int nacc_max,     // Total Coherent accumulations
                   int ibinbgn,      // Min. Starting bin for computations
                   int ibinend) {    // Max. Stop bin for computations

    float[] csyy = new float[2];   // Complex Syy spectrum

    float[] csxy = dcsxy.castF(true);
    float[] fsyy = dfsyy.castF(true);
    float[] csxx = dcsxx.castF(true);

    //	Compute start/stop bin locations

    int nctap = Cinf.ntaps/2;
    int nbck  = Cinf.nshift - (Cinf.nshift/2);
    int noff  = -Cinf.nshift-nctap;
    int istrt = Math.max(Cinf.nshift,ibinbgn+nbck);
    int istop = Math.min(ibinend+nbck,n);
    istrt     = Math.min(istrt,istop);

    if (0 != (debug & DEBUG_COHER)) {
      M.info ("accum_coher:");
      M.info (" ibinbgn "+ibinbgn);
      M.info (" ibinend "+ibinend);
      M.info (" n       "+n);
      M.info (" ntaps   "+Cinf.ntaps);
      M.info (" nshift  "+Cinf.nshift);
      M.info (" nctap   "+nctap);
      M.info (" nbck    "+nbck);
      M.info (" noff    "+noff);
      M.info (" istrt   "+istrt);
      M.info (" istop   "+istop);
      M.info (" min     "+(nctap-1+istrt+noff));
      M.info (" max     "+(nctap+1+istop+noff));
      M.info ("nacc_max "+nacc_max);
    }

    // Single Coherent Accumulation (no Syy required)
    if (nacc_max <= 1) {

      // Ensure we stay within the bounds of the arrays

      if (istrt-nbck < 0) {
        istrt = nbck;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istrt "+istrt);
      }
      if (istop-nbck >= n) {
        istop = n-1+nbck;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istop "+istop);
      }
      if ((nctap-1+istrt+noff) < 0) {
        istrt = 1-nctap-noff;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istrt, min "+istrt+", "+(nctap-1+istrt+noff));
      }
      if ((nctap+1+istop+noff) >= n) {
        istop = n-nctap-2-noff;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istop, max "+istop+", "+(nctap+1+istop+noff));
      }
      if (istop >= n) {
        istop = n-1;
      }

      if (0 != (debug & DEBUG_COHER)) {
        M.info ("centered "+(istrt+nctap+noff)+", "+(istop+nctap+noff));
      }

      if (!log_init_sxy) {
        for (int i = istrt; i <= istop; i++) {
          csyy[0]  = Cinf.cxlo[2*(nctap)  ]*csxx[2*(nctap+i+noff)  ]
                   - Cinf.cxlo[2*(nctap)+1]*csxx[2*(nctap+i+noff)+1];
          csyy[1]  = Cinf.cxlo[2*(nctap)  ]*csxx[2*(nctap+i+noff)+1]
                   + Cinf.cxlo[2*(nctap)+1]*csxx[2*(nctap+i+noff)  ];

          csyy[0] += Cinf.cxlo[2*(nctap-1)  ]*csxx[2*(nctap-1+i+noff)  ]
                   - Cinf.cxlo[2*(nctap-1)+1]*csxx[2*(nctap-1+i+noff)+1];
          csyy[1] += Cinf.cxlo[2*(nctap-1)  ]*csxx[2*(nctap-1+i+noff)+1]
                   + Cinf.cxlo[2*(nctap-1)+1]*csxx[2*(nctap-1+i+noff)  ];

          csyy[0] += Cinf.cxlo[2*(nctap+1)  ]*csxx[2*(nctap+1+i+noff)  ]
                   - Cinf.cxlo[2*(nctap+1)+1]*csxx[2*(nctap+1+i+noff)+1];
          csyy[1] += Cinf.cxlo[2*(nctap+1)  ]*csxx[2*(nctap+1+i+noff)+1]
                   + Cinf.cxlo[2*(nctap+1)+1]*csxx[2*(nctap+1+i+noff)  ];

          float csyyre =  csxx[2*i]*csyy[0] + csxx[2*i+1]*csyy[1];
          float csyyim = -csxx[2*i]*csyy[1] + csxx[2*i+1]*csyy[0];
          csyy[0] = csyyre*crotat[0] - csyyim*crotat[1];
          csyy[1] = csyyre*crotat[1] + csyyim*crotat[0];

          float abscsyymax = Math.max((float) Math.sqrt(csyy[0]*csyy[0] + csyy[1]*csyy[1]), floorval);
          csxy[2*(i-nbck)  ] += csyy[0]/abscsyymax;
          csxy[2*(i-nbck)+1] += csyy[1]/abscsyymax;
        }
      }
      else {
        for (int i = istrt; i <= istop; i++) {
          csyy[0]  = Cinf.cxlo[2*nctap  ]*csxx[2*(nctap+i+noff)  ]
                   - Cinf.cxlo[2*nctap+1]*csxx[2*(nctap+i+noff)+1];
          csyy[1]  = Cinf.cxlo[2*nctap  ]*csxx[2*(nctap+i+noff)+1]
                   + Cinf.cxlo[2*nctap+1]*csxx[2*(nctap+i+noff)  ];

          csyy[0] += Cinf.cxlo[2*(nctap-1)  ]*csxx[2*(nctap-1+i+noff)  ]
                   - Cinf.cxlo[2*(nctap-1)+1]*csxx[2*(nctap-1+i+noff)+1];
          csyy[1] += Cinf.cxlo[2*(nctap-1)  ]*csxx[2*(nctap-1+i+noff)+1]
                   + Cinf.cxlo[2*(nctap-1)+1]*csxx[2*(nctap-1+i+noff)  ];

          csyy[0] += Cinf.cxlo[2*(nctap+1)  ]*csxx[2*(nctap+1+i+noff)  ]
                   - Cinf.cxlo[2*(nctap+1)+1]*csxx[2*(nctap+1+i+noff)+1];
          csyy[1] += Cinf.cxlo[2*(nctap+1)  ]*csxx[2*(nctap+1+i+noff)+1]
                   + Cinf.cxlo[2*(nctap+1)+1]*csxx[2*(nctap+1+i+noff)  ];

          float csyyre  = csxx[2*i]*csyy[0] + csxx[2*i+1]*csyy[1];
          float csyyim = -csxx[2*i]*csyy[1] + csxx[2*i+1]*csyy[0];
          csyy[0] = csyyre*crotat[0] - csyyim*crotat[1];
          csyy[1] = csyyre*crotat[1] + csyyim*crotat[0];

          float abscsyymax = Math.max((float) Math.sqrt(csyy[0]*csyy[0] + csyy[1]*csyy[1]), floorval);
          csxy[2*(i-nbck)  ] = csyy[0]/abscsyymax;
          csxy[2*(i-nbck)+1] = csyy[1]/abscsyymax;
        }
      }
    }
    // Multiple Coherent accumulations (Syy and Sxx required)
    else {

      // Ensure we stay within the array bounds
      if ((istrt-nbck) < 0) {
        istrt = nbck;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istrt "+istrt);
      }
      if ((istop-nbck) >= n) {
        istop = n-1+nbck;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istop "+istop);
      }
      if ((istrt+noff) < 0) {
        istrt = -noff;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istrt, min "+istrt+", "+(istrt+noff));
      }
      if ((Cinf.ntaps+istop+noff) >= n) {
        istop = n - Cinf.ntaps - noff;
        if (0 != (debug & DEBUG_COHER))
          M.info ("coher: new istop, max "+istop+","+(Cinf.ntaps+istop+noff));
      }
      if (istop >= n) {
        istop = n-1;
      }

      if (0 != (debug & DEBUG_COHER)) {
        M.info ("centered "+(istrt+nctap+noff)+", "+(istop+nctap+noff));
      }

      if (!log_init_sxy) {
        for (int i = istrt; i <= istop; i++) {
          csyy[0] = 0.0f;
          csyy[1] = 0.0f;
          for (int ii = 0; ii < Cinf.ntaps; ii++) {
            csyy[0] += Cinf.cxlo[2*ii  ]*csxx[2*(ii+i+noff)  ]
                     - Cinf.cxlo[2*ii+1]*csxx[2*(ii+i+noff)+1];
            csyy[1] += Cinf.cxlo[2*ii  ]*csxx[2*(ii+i+noff)+1]
                     + Cinf.cxlo[2*ii+1]*csxx[2*(ii+i+noff)  ];
          }
          float csyyre =  csxx[2*i]*csyy[0] + csxx[2*i+1]*csyy[1];
          float csyyim = -csxx[2*i]*csyy[1] + csxx[2*i+1]*csyy[0];

          csxy[2*(i-nbck)  ] += csyyre*crotat[0] - csyyim*crotat[1];
          csxy[2*(i-nbck)+1] += csyyre*crotat[1] + csyyim*crotat[0];

          fsyy[i] += (csyy[0]*csyy[0] + csyy[1]*csyy[1]);
        }
      }
      else {
        for (int i = istrt; i <= istop; i++) {
          csyy[0] = 0.0f;
          csyy[1] = 0.0f;
          for (int ii = 0; ii < Cinf.ntaps; ii++) {
            csyy[0] += Cinf.cxlo[2*ii  ]*csxx[2*(ii+i+noff)  ]
                     - Cinf.cxlo[2*ii+1]*csxx[2*(ii+i+noff)+1];
            csyy[1] += Cinf.cxlo[2*ii  ]*csxx[2*(ii+i+noff)+1]
                     + Cinf.cxlo[2*ii+1]*csxx[2*(ii+i+noff)  ];
          }
          float csyyre =  csxx[2*i]*csyy[0] + csxx[2*i+1]*csyy[1];
          float csyyim = -csxx[2*i]*csyy[1] + csxx[2*i+1]*csyy[0];

          csxy[2*(i-nbck)  ] = csyyre*crotat[0] - csyyim*crotat[1];
          csxy[2*(i-nbck)+1] = csyyre*crotat[1] + csyyim*crotat[0];

          fsyy[i] = (csyy[0]*csyy[0] + csyy[1]*csyy[1]);
        }
      }
    }

    if (0 != (debug & DEBUG_COHER)) {
      M.info ("accum_coher:");
      M.info ("  sxy");
      for (int j = 0; j < istrt+10; j++) {
        M.info ("\t"+j+"\t("+csxy[2*j]+","+csxy[2*j+1]+")");
      }
      for (int j = istop-10; j<istop; j++) {
        M.info ("\t"+j+"\t("+csxy[2*j]+","+csxy[2*j+1]+")");
      }
      if (nacc_coher_max > 1) {
        M.info ("  syy");
        for (int j = 0; j < istrt+10; j++) {
          M.info ("\t"+j+"\t"+fsyy[j]);
        }
        for (int j = istop-10; j<istop; j++) {
          M.info ("\t"+j+"\t"+fsyy[j]);
        }
      }
    }

    dcsxy.uncast (csxy, true);
    dfsyy.uncast (fsyy, true);
    dcsxx.uncast(csxx, true);

    return n;
  }

  /**
   ************************************************************************
   * FUNCTION ACCUM_NORM
   ************************************************************************
   *	ACCUM_NORM --- Accumulates coherence spectra
   */
  int accum_norm (Data dcsxy,           // in  - Complex Sxy spectrum
                  Data dfsxx,           // in  - Auto spectrum
                  Data dfsyy,           // in  - Auto spectrum of shifted signal
                  Data dcnmsxy,         // out - Normalized Complex Sxy spectrum
                  CycRecord Cinf,       // CYC record information
                  int n,                // Number of points in input/output
                  int ncoher_max,       // Total Normalized accumulations
                  int nnorm_max,        // Total Normalized accumulations
                  float rwhite,         // Whitening coefficient
                  int ibinbgn,          // Min. Starting bin for computations
                  int ibinend)          // Max. Stop bin for computations
  {

    float[] csxy = dcsxy.castF(true);
    float[] fsxx = dfsxx.castF(true);
    float[] fsyy = dfsyy.castF(true);
    float[] cnmsxy = dcnmsxy.castF(true);

    // Single Coherent Accumulation (no Syy needed)
    if (ncoher_max <= 1) {
      float rscale = 1.0f/((float) nnorm_max*Math.max(1,ncoher_max));
      M.info("rscale="+rscale);
      if (!log_init_nmsxy) {
        for (int i = ibinbgn; i <= ibinend; i++) {
          cnmsxy[2*i  ] += rscale*csxy[2*i];
          cnmsxy[2*i+1] += rscale*csxy[2*i+1];
        }
      }
      else {
        for (int i = ibinbgn; i <= ibinend; i++) {
          cnmsxy[2*i  ] = rscale*csxy[2*i];
          cnmsxy[2*i+1] = rscale*csxy[2*i+1];
        }
      }
    }

    // Multiple Coherent Accumulations (need Sxx and Syy)
    else {
      int nbck = Math.max(1,Cinf.nshift-Cinf.nshift/2);
      int istrt = ibinbgn;
      int istop = Math.min(ibinend, (n-nbck+1));
      istrt = Math.min (istrt, istop);
      float rscale = 1.0f/((float) Math.max (1,nnorm_max));
      // Stay within array bounds
      if (istrt+nbck < 0) {
        istrt = -nbck;
        if (0 != (debug & DEBUG_NCOHER))
          M.info ("ncoher: new istrt "+istrt);
      }
      if (istop+nbck >= n) {
        istop = n-nbck-1;
        if (0 != (debug & DEBUG_NCOHER))
          M.info ("ncoher: new istop "+istop);
      }
      if (!log_init_nmsxy) {
        for (int i = istrt; i <= istop; i++) {
          float fsxymax = (float) Math.sqrt(fsxx[i+nbck]*fsyy[i+nbck]+rwhite);
          fsxymax = Math.max (fsxymax, floorval);
          cnmsxy[2*i  ] += rscale*csxy[2*i]/fsxymax;
          cnmsxy[2*i+1] += rscale*csxy[2*i+1]/fsxymax;
        }
      }
      else {
        for (int i = istrt; i <= istop; i++) {
          float fsxymax = (float) Math.sqrt(fsxx[i+nbck]*fsyy[i+nbck]+rwhite);
          fsxymax = Math.max (fsxymax, floorval);
          cnmsxy[2*i  ] = rscale*csxy[2*i]/fsxymax;
          cnmsxy[2*i+1] = rscale*csxy[2*i+1]/fsxymax;
        }
      }
    }

    if (0 != (debug & DEBUG_NCOHER)) {
      M.info ("accum_norm cnmsxy");
      M.info ("rwhiten_avg "+rwhiten_avg);
      for (int j = 0; j < 10; j++)
        M.info ("\t"+j+"\t("+cnmsxy[2*j]+","+cnmsxy[2*j+1]+")");
      for (int j = n-10; j < n; j++)
        M.info ("\t"+j+"\t("+cnmsxy[2*j]+","+cnmsxy[2*j+1]+")");
    }

    dcsxy.uncast(csxy, true);
    dfsxx.uncast(fsxx, true);
    dfsyy.uncast(fsyy, true);
    dcnmsxy.uncast(cnmsxy, true);

    return n;
  }

  /**
   ***********************************************************************
   * FUNCTION CLEAN_CYCSPEC
   ************************************************************************
   * CLEAN_CYCSPEC --- smooths the cyclic spectrum and computes the cyclic
   *		    power spectrum (real) from the complex valued spectrum
   */
  int clean_cycspec (Data     dcspec,
                     Data     dfspec,
                     int      n,
                     int      ibgn,
                     int      iend,
                     int      nwin,
                     boolean  log_dosmooth,
                     boolean  log_dopower)
  {

    float[] cspec = dcspec.castF(true);
    float[] fspec = dfspec.castF(true);

    // Do smoothing and possibly power conversion
    if (log_dosmooth) {
      // Calculate smoothing factors
      float fbeta  = (float) (1.0/Math.max(1.0,Math.sqrt(sqrt2o2*((double)nwin))));
      float fbetai = 1.0f-fbeta;
      // Filter reverse for zero phase
      float[] ctemp = new float[2];
      ctemp[0] = 0.0f;
      ctemp[1] = 0.0f;
      for (int i = n-1; i >= 0; i--) {
        ctemp[0] = fbetai*ctemp[0] + fbeta*cspec[2*i  ];
        ctemp[1] = fbetai*ctemp[1] + fbeta*cspec[2*i+1];

        cspec[2*i] = ctemp[0];
        cspec[2*i+1] = ctemp[1];
      }
      // Filter forward
      if (log_dopower) {
        for (int i = 0; i < n; i++) {
          ctemp[0] = fbetai*ctemp[0] + fbeta*cspec[2*i  ];
          ctemp[1] = fbetai*ctemp[1] + fbeta*cspec[2*i+1];
          fspec[i] = (float) Math.sqrt(ctemp[0]*ctemp[0] * ctemp[1]*ctemp[1]);
        }
        for (int i = iend+1; i < n; i++)
          fspec[i] = 0.0f;
        for (int i = 0; i < ibgn; i++)
          fspec[i] = 0.0f;
      }
      else {
        for (int i = 0; i < n; i++) {
          ctemp[0] = fbetai*ctemp[0] + fbeta*cspec[2*i  ];
          ctemp[1] = fbetai*ctemp[1] + fbeta*cspec[2*i+1];
          cspec[2*i  ] = ctemp[0];
          cspec[2*i+1] = ctemp[1];
        }
        for (int i = iend+1; i < n; i++) {
          cspec[2*i  ] = 0.0f;
          cspec[2*i+1] = 0.0f;
        }
        for (int i = 0; i < ibgn; i++) {
          cspec[2*i  ] = 0.0f;
          cspec[2*i+1] = 0.0f;
        }
      }
    }

    // Just do power conversion
    if ((log_dopower) && (!log_dosmooth)) {
      for (int i = 0; i < ibgn; i++)
        fspec[i] = 0.0f;
      for (int i = ibgn; i <= iend; i++) {
        fspec[i] = (float) Math.sqrt(cspec[2*i]*cspec[2*i] + cspec[2*i+1]*cspec[2*i+1]);
      }
      for (int i = iend+1; i < n; i++)
        fspec[i] = 0.0f;
    }

    dcspec.uncast (cspec, true);
    dfspec.uncast (fspec, true);

    return n;
  }

}
