package nxm.dsp.prim;

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

/** Port of the Fortran host primitive NEWRATE from X-Midas. Please see
    explain file for a detailed description of the "Principal of Operation"
    and "Program Design".<br>
    <br>
    <b>Implementation Notes:</b><br>
    <br>
    Note this port was done in Feb 2007 based on the 2.8 version of the Fortran
    code for NEWRATE. This activity predated the C++ implementation of NEWRATE
    in X-Midas by several years (hence why the C++ version was not used as the
    source).
    @author  John Marsh
 */
public class newrate extends Primitive {
  private static final int INBUFSIZE   = 128*1024;   // input  buffer size (floats)
  private static final int OUTBUFSIZE  = 128*1024;   // output buffer size (floats)
  private static final int FILTBUFSIZE = 256*1024;   // filter buffer size (floats)
  private static final int NTAPS       = 28;         // # taps in default filter
  private static final int NFILTERS    = 2048;       // # default filters in table

  //  Sync file additions
  private static final int PBUFSIZE_MAX = 16384;
  private  double t_now = 0.0;
  private  double tstep = 0.0;
  private  int nstep = 0;
  private  double bufstart = 0;
  private  double x_start = 0;
  private  double pstep_nominal = 0;
  private  int npts = 0;
  private  double p_start = 0;
  private  double pp_delta = 0;
  private  double pri = 0;
  private  boolean nowrap = false;
  private  int ncleft = 0;
  private  boolean matchXM = false;
  //  End Sync file additions

  private DataFile    inputFile;
  private DataFile    outputFile;
  private DataFile    syncFile;
  private DataFile    filterFile;
  private Data        inputData;
  private Data        outputData;
  private Data        filterData;
  private Data        syncData;

  private float[][]   filterBuffer;
  private float[]     inBuffer      = new float[INBUFSIZE];
  private float[]     outBuffer;


  /**Sample index at start of input buffer*/      private long    inBase;
  /**Current input sample read ptr (NB double!)*/ private double  inReadP;
  /**Start sample for current loop*/              private double  inReadPStart;
  /**Current input (sample) write pointer*/       private long    inWriteP;
  /**size of input buffer in samples*/            private int     inBufSamples;
  /**Sample index at start of outputbuffer*/      private long    outBase;
  /**Current output (sample) write pointer*/      private long    outWriteP;
  /**Size of output buffer in samples*/           private int     outBufSamples;
  /**Size of phase buffer in samples*/            private int     phaseBufSamples;

  /**TUNE switch present*/                        private boolean tuning;
  /**SYNC file supplied*/                         private boolean syncing;
  /**OFFSET switch present*/                      private boolean offsetting;
  /**Complex Filter?*/                            private boolean cxFilter;

  /**target sample rate*/                         private double  newRate;
  /**Phase offset to add to outout*/              private double  phaseOffset  = 0.0;
  /**Complex tuner step [ tune_cx ]*/             private double  rTune;
                                                  private double  qTune;
  /**Complex oscillator [ tune_twiddle ]*/        private double  rTuneTwiddle = 1.0;
                                                  private double  qTuneTwiddle = 0.0;
  /**Force round-to-nearest if filter complex*/   private double  round        = 0.0;
  /**step between input samples*/                 private double  inStep;

  /**number of taps per filter    [ ntaps ]*/     private int     nTaps;
  /**number interpolation filters [ ninterps ]*/  private int     nInterps;
  /**number output samples we have room for*/     private int     outRoomLeft;
  /**Input floats per sample (R=1,CX=2)*/         private int     inFPS;
  /**Output floats per sample (R=1,CX=2)*/        private int     outFPS;
  /**Filter floats per sample (R=1,CX=2)*/        private int     filterFPS;
  /**number outSteps taken (num iters) so far*/   private long    outCount     = 0;
  /**samples available for resampling*/           private int     inHave;

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

  @Override
  public int open () {
    // Input File - REAL, Type 1000
    inputFile = MA.getDataFile ("IN", "1000", "-#",0);
    inputFile.open();
    inputData = inputFile.getDataBuffer(INBUFSIZE, FLOAT);

    if (inputFile.getFormatMode() == COMPLEX) inFPS = 2;
    else if (inputFile.getFormatMode() == SCALAR) inFPS = 1;
    else throw new MidasException("Input File Format must be C(omplex) or S(calar)");

    // Output File - output type based on input type
    outputFile = MA.getDataFile ("OUT", inputFile, 0);
    outputFile.open();
    outputData = outputFile.getDataBuffer(OUTBUFSIZE, FLOAT);
    outBuffer  = outputData.castF(false);
    outBase    = 0;
    outWriteP  = 0;


    nowrap = MA.find("/NOWRAP");
    matchXM = MA.find("/XMTO");
    //System.out.print("nowrap " );
    //System.out.println( nowrap );
    if( nowrap ) {
      throw new MidasException( "NOWRAP not yet implemented.");
    }

    // Target Sample Rate
    setRate(MA.getD("RATE"));

    // Filter File (optional) - type 2000, real or complex
    if (MA.find("FILTER")) {
      filterFile = MA.getDataFile ("FILTER", "2000", "-#", 0);
      filterFile.open();

      if (filterFile.getFormatMode() == COMPLEX) {
        outputFile.setFormatMode(COMPLEX);
        cxFilter  = true;
        filterFPS = 2;
      } else if (filterFile.getFormatMode() != SCALAR) {
        throw new MidasException("Unsupported filter file format");
      }
      else {
        filterFPS = 1;
      }
      userFilter();
    } else {   // Use default built-in 28-tap filter
      cxFilter = false;
      nInterps = 2048;
      defaultFilter();
    }

    if (outputFile.getFormatMode() == COMPLEX) {
      outFPS = 2;
    } else {
      outFPS = 1;
    }

    // /TUNE=<tuneFreq> - double value
    //                    If present, use it to construct a "complex oscillator"
    //                    to apply to the output.  Each output sample is multiplied
    //                    by the complex [rTuneTwiddle, qTuneTwiddle], which itself
    //                    is multiplied by [rTune, qTune] ready for the next
    //                    output sample.
    tuning = MA.find("/TUNE");
    if (tuning) {
      setTune(MA.getD("/TUNE"));
    }

    // /OFFSET=<offsetPhase> - double value
    //                         If present, ...
    offsetting = MA.find("/OFFSET");
    if (offsetting) {
      phaseOffset = MA.getD("/OFFSET");
      phaseOffset -= (long)phaseOffset;
      if (phaseOffset < 0.0) phaseOffset++;
    }

    outputFile.setXStart(inputFile.getXStart() +
        (nTaps/2.0 - 1.0)*inputFile.getXDelta() +
        phaseOffset * outputFile.getXDelta() );

    // Sync File (optional) - type 1000
    syncing = MA.find("SYNC");
    if (syncing) {
      syncFile = MA.getDataFile ("SYNC", "1000", "S#", 0);
      syncFile.open();
      pri = 1.0 / newRate;

      phaseBufSamples = PBUFSIZE_MAX;
      firstClock( outputFile.getXStart(), 0.0 );
      ncleft = nstep;
      inReadP=( t_now + outCount * outputFile.getXDelta() + p_start - outputFile.getXStart()) / inputFile.getXDelta();
      inStep = (tstep + outputFile.getXDelta())/inputFile.getXDelta();
      inHave = 0;
    }
    else {
      // inStep is the number of input samples to step for each output sample
      // this is NOT an integer.
      inReadP = inStep * phaseOffset;
      inHave  = 0;
    }

    inBufSamples  = INBUFSIZE  / inFPS;
    outBufSamples = OUTBUFSIZE / outFPS;
    outRoomLeft   = outBufSamples;
    inWriteP      = 0;

    return NORMAL;
  }

  @Override
  public int process () {
    // process() produces samples at a constant output rate. Sync data is
    // used to calculate an adjusted sample rate between sync points so that
    // the clock phase matches the sync phase at each sync point. Output data
    // values are effectively interpolations of the input data samples. The
    // number of interpolated positions is determined by the filter size, and
    // the size of the filter table. The default filter uses 2048 interpolation
    // positions. Internal buffer management routines retain needed input data.
    // Output is accumulated to reduce the overhead write operations. The actual
    // filtering loop uses the maximin possible number of iterations, subject to
    // input data available, room in the output buffer, and the possible need
    // to set a new sample rate at the next phase point.
    //
    // Care has been taken to minimize cumulative errors in abscissa values when
    // processing large files, and to reduce the number of operations in the
    // inner processing loop. This has resulted in some increase in program
    // complexity. Comments in the data declaration section indicate the
    // intended function of some of the control variables.
    int     status  = NORMAL;
    int     iters   = 0;

    // Input buffer management:  In general we have samples
    // already in the buffer, some of which have been resampled, and some
    // of which are needed for future resampling.
    //
    // 1.  If there is no room in the input buffer, then we copy any samples that
    //     we might need for future resampling to the beginning of the buffer.
    //     This makes room in the buffer.
    if ((inWriteP-inBase == inBufSamples)) {
      int offset = (int)(inReadP-inBase) * inFPS; // note: when inStep > nTaps, size could be negative - no copying needed
      int size   = inHave * inFPS;
      for (int ii=0; ii<size; ii++) {
        inBuffer[ii] = inBuffer[ii + offset];
      }
      inBase = inWriteP - inHave;
    }

    // 2.  If there is room in the buffer for more samples, AND we need more
    //     then we try to fill the buffer.
    if (inWriteP-inBase < inBufSamples && inHave <= nTaps) {
      // When inStep > nTaps, we always want to read in as much as we can (not keeping any old input samples)
      long read  = (inHave < 0)? inBufSamples : inBufSamples - (inWriteP-inBase);
      int  inGot = inputFile.read(inputData, (int)read);
      if (inGot < 0 && inHave < nTaps) {
        status = FINISH;
      }
      else {
        int     offset = (int)(inWriteP-inBase) * inFPS;
        int     size   = inGot * inFPS;
        float[] input  = inputData.castF(true);

        int numSkip = 0; // amount of data we need to receive just to get up to inReadP
        if (inHave < 0) {
          // Handle case where inStep > nTaps and inReadP is ahead if inWriteP
          // Skip data up to the sample matching inReadP and start writing at beginning of inBufffer
          // - since 3.9.4 for NXM-3613
          numSkip = -inHave * inFPS;
          size -= numSkip;  // reduce size by input points we are skipping
          offset = 0;
          inBase = (long)inReadP;
        }
        // We could have the scenario where not enough data was read in to get to
        // the next sample we want to process, so check size before copying
        if (size > 0) System.arraycopy(input, numSkip, inBuffer, offset, size);
        inputData.uncast(input, false);
        inWriteP += inGot;
        inHave    = (int)(inWriteP - (long)inReadP);
      }
    }

    // 3.  If we have enough samples to resample (more than nTaps), then
    //     proceed to perform as many iterations as possible
    if (inHave < nTaps && status == NORMAL) { // not enough data to do anything with
      status = NOOP;
    }
    if (inHave >= nTaps) {
      // We can now iterate the interpolation process until the earliest of:
      //   (1) the number of points left in the output buffer,
      //   (2) the number of points for which we have input data, or
      //   (3) the number of steps left for the current syncfile interval.

      // FORTRAN code to calculate number of iterations:
      // nrep = min(blksize2-indexb,int((siglast-aindex+.999)/astep)+1)
      double steps = (inWriteP-inReadP-nTaps+.999999)/inStep;
      iters = Math.min(outRoomLeft, (int)steps+1);

      if (syncing) {
        iters = Math.min( iters , ncleft );

      }
      inReadPStart = inReadP;
      applyFilter(iters, inStep);
      outRoomLeft -= iters;
    }

    // Avoid writing small blocks by waiting until we are at least 90% full
    if (status == FINISH || (outWriteP - outBase) > outBufSamples * .9) {
      outputData.uncast(outBuffer, true);
      outputFile.write(outputData, (int)(outWriteP-outBase));
      outBuffer   = outputData.castF(false);
      outBase     = outWriteP;
      outRoomLeft = outBufSamples;
    }

    // Some preparation for next process() call
    //inReadP  += iters * inStep;
    inReadP = inReadPStart + ( iters * inStep );
    outCount += iters;
    inHave = (int)(inWriteP - (long)inReadP);

    if( syncing ) {
      ncleft = ncleft - iters;
      if( ncleft < 1 ) {
        t_now = t_now + nstep * tstep;
        nextClock();
        ncleft = nstep;
        inReadP = ( t_now + outCount*outputFile.getXDelta() + p_start - outputFile.getXStart())/inputFile.getXDelta();
        inStep = ( tstep + outputFile.getXDelta()) / inputFile.getXDelta();
      }
    }

    return status;
  }

  @Override
  public int close () {
    inputFile.close();
    outputData.uncast(outBuffer, false);
    outputFile.close();
    if (syncFile   != null) {
      //System.out.println("Closing sync file." );
      syncFile.close();
    }
    if (filterFile != null) filterFile.close();

    return NORMAL;
  }

  /** Built-in table of 17x28 coefficients for first half of filter used by defaultFilter(). */
  private static final double[][] HFILTER = {
    { -0.1723066E-12,  0.2602363E-12, -0.4423129E-12,  0.5346834E-12,
      -0.6608047E-12,  0.8437695E-12, -0.9716672E-12,  0.1005418E-11,
      -0.1132650E-11,  0.1079137E-11, -0.1073808E-11,  0.1087463E-11,
      -0.1034284E-11,  1.000000     , -0.1008971E-11,  0.1066258E-11,
      -0.1106115E-11,  0.1100453E-11, -0.1066702E-11,  0.1133982E-11,
      -0.9734435E-12,  0.8544276E-12, -0.6910028E-12,  0.5044853E-12,
      -0.4209966E-12,  0.3179679E-12, -0.1270095E-12,  0.1190159E-12 },
      { -0.4104828E-03,  0.5881413E-03, -0.8470448E-03,  0.1188470E-02,
        -0.1626680E-02,  0.2186433E-02, -0.2900885E-02,  0.3831123E-02,
        -0.5082040E-02,  0.6867848E-02, -0.9688202E-02,  0.1503587E-01,
        -0.3021841E-01,  0.9986488    ,  0.3166670E-01, -0.1503463E-01,
        0.9438746E-02, -0.6555507E-02,  0.4764554E-02, -0.3532743E-02,
        0.2634133E-02, -0.1954591E-02,  0.1433210E-02, -0.1028657E-02,
        0.7162614E-03, -0.4755949E-03,  0.2884227E-03, -0.4237333E-03 },
        { -0.7187069E-03,  0.9908061E-03, -0.1444174E-02,  0.2060965E-02,
          -0.2865905E-02,  0.3907888E-02, -0.5254533E-02,  0.7026941E-02,
          -0.9431731E-02,  0.1288739E-01, -0.1836151E-01,  0.2871460E-01,
          -0.5760054E-01,  0.9933048    ,  0.6595361E-01, -0.3120765E-01,
          0.1975870E-01, -0.1388523E-01,  0.1022675E-01, -0.7694323E-02,
          0.5829360E-02, -0.4402320E-02,  0.3291889E-02, -0.2416211E-02,
          0.1728239E-02, -0.1187160E-02,  0.7531443E-03, -0.9410921E-03 },
          { -0.9825213E-03,  0.1322678E-02, -0.1954004E-02,  0.2824262E-02,
            -0.3969249E-02,  0.5462254E-02, -0.7402023E-02,  0.9968592E-02,
            -0.1346557E-01,  0.1850434E-01, -0.2648273E-01,  0.4149505E-01,
            -0.8264041E-01,  0.9845923    ,  0.1021276    , -0.4779923E-01,
            0.3030100E-01, -0.2138262E-01,  0.1583498E-01, -0.1198874E-01,
            0.9144412E-02, -0.6956737E-02,  0.5244220E-02, -0.3886144E-02,
            0.2811235E-02, -0.1959618E-02,  0.1271215E-02, -0.1486347E-02 },
            { -0.1229729E-02,  0.1633989E-02, -0.2442849E-02,  0.3558575E-02,
              -0.5028745E-02,  0.6950616E-02, -0.9456534E-02,  0.1278278E-01,
              -0.1731891E-01,  0.2385583E-01, -0.3418219E-01,  0.5348989E-01,
              -0.1054820    ,  0.9727868    ,  0.1398093    , -0.6444895E-01,
              0.4075991E-01, -0.2877768E-01,  0.2134885E-01, -0.1620298E-01,
              0.1239382E-01, -0.9460293E-02,  0.7155922E-02, -0.5322010E-02,
              0.3867490E-02, -0.2714455E-02,  0.1784308E-02, -0.2023978E-02 },
              { -0.1473194E-02,  0.1938066E-02, -0.2921977E-02,  0.4275934E-02,
                -0.6060412E-02,  0.8396355E-02, -0.1144337E-01,  0.1548994E-01,
                -0.2100578E-01,  0.2894415E-01, -0.4145028E-01,  0.6467696E-01,
                -0.1261328    ,  0.9580158    ,  0.1787723    , -0.8094347E-01,
                0.5097277E-01, -0.3594548E-01,  0.2666636E-01, -0.2024795E-01,
                0.1550026E-01, -0.1184472E-01,  0.8971456E-02, -0.6684138E-02,
                0.4867971E-02, -0.3428578E-02,  0.2271249E-02, -0.2533638E-02 },
                { -0.1720918E-02,  0.2237771E-02, -0.3391272E-02,  0.4974463E-02,
                  -0.7060205E-02,  0.9790856E-02, -0.1334940E-01,  0.1807334E-01,
                  -0.2450471E-01,  0.3373951E-01, -0.4824132E-01,  0.7499048E-01,
                  -0.1445622    ,  0.9403764    ,  0.2188271    , -0.9711936E-01,
                  0.6083014E-01, -0.4280343E-01,  0.3172203E-01, -0.2407327E-01,
                  0.1842601E-01, -0.1408198E-01,  0.1066863E-01, -0.7953897E-02,
                  0.5798256E-02, -0.4090765E-02,  0.2721652E-02, -0.3000148E-02 },
                  { -0.1978698E-02,  0.2539513E-02, -0.3856772E-02,  0.5658871E-02,
                    -0.8029593E-02,  0.1113157E-01, -0.1516861E-01,  0.2052077E-01,
                    -0.2779895E-01,  0.3822265E-01, -0.5452402E-01,  0.8438267E-01,
                    -0.1607647    ,  0.9199924    ,  0.2597637    , -0.1128011    ,
                    0.7021537E-01, -0.4926194E-01,  0.3644509E-01, -0.2762774E-01,
                    0.2113172E-01, -0.1613911E-01,  0.1222079E-01, -0.9108335E-02,
                    0.6638190E-02, -0.4684233E-02,  0.3121979E-02, -0.3410907E-02 },
                    { -0.2246133E-02,  0.2837869E-02, -0.4311746E-02,  0.6319633E-02,
                      -0.8955512E-02,  0.1240187E-01, -0.1688161E-01,  0.2280971E-01,
                      -0.3085708E-01,  0.4234878E-01, -0.6024230E-01,  0.9279005E-01,
                      -0.1747232    ,  0.8969753    ,  0.3013911    , -0.1278384    ,
                      0.7903950E-01, -0.5526447E-01,  0.4079632E-01, -0.3087931E-01,
                      0.2359173E-01, -0.1799905E-01,  0.1361728E-01, -0.1014043E-01,
                      0.7382850E-02, -0.5205662E-02,  0.3470977E-02, -0.3759954E-02 },
                      { -0.2521820E-02,  0.3126621E-02, -0.4746639E-02,  0.6946272E-02,
                        -0.9827083E-02,  0.1358697E-01, -0.1846842E-01,  0.2491469E-01,
                        -0.3364485E-01,  0.4607555E-01, -0.6534874E-01,  0.1001599    ,
                        -0.1864318    ,  0.8714455    ,  0.3435165    , -0.1420850    ,
                        0.8721583E-01, -0.6075507E-01,  0.4474083E-01, -0.3380381E-01,
                        0.2578715E-01, -0.1964824E-01,  0.1484772E-01, -0.1104254E-01,
                        0.8029398E-02, -0.5655342E-02,  0.3768718E-02, -0.4042819E-02 },
                        { -0.2800683E-02,  0.3400850E-02, -0.5156998E-02,  0.7531353E-02,
                          -0.1063397E-01,  0.1467525E-01, -0.1991343E-01,  0.2681384E-01,
                          -0.3613512E-01,  0.4936980E-01, -0.6980349E-01,  0.1064519    ,
                          -0.1959051    ,  0.8435427    ,  0.3859360    , -0.1553897    ,
                          0.9465915E-01, -0.6568033E-01,  0.4824200E-01, -0.3637554E-01,
                          0.2769969E-01, -0.2107188E-01,  0.1590099E-01, -0.1180826E-01,
                          0.8573212E-02, -0.6029339E-02,  0.4014550E-02, -0.4259268E-02 },
                          { -0.3079207E-02,  0.3654434E-02, -0.5534245E-02,  0.8064164E-02,
                            -0.1136279E-01,  0.1564901E-01, -0.2119403E-01,  0.2848180E-01,
                            -0.3830048E-01,  0.5219878E-01, -0.7356650E-01,  0.1116262    ,
                            -0.2031609    ,  0.8134087    ,  0.4284509    , -0.1676149    ,
                            0.1012989    , -0.6999786E-01,  0.5127081E-01, -0.3857515E-01,
                            0.2931919E-01, -0.2226568E-01,  0.1677460E-01, -0.1243608E-01,
                            0.9014320E-02, -0.6328789E-02,  0.4209657E-02, -0.4408102E-02 },
                            { -0.3351431E-02,  0.3882369E-02, -0.5872210E-02,  0.8536462E-02,
                              -0.1200258E-01,  0.1649466E-01, -0.2229383E-01,  0.2989854E-01,
                              -0.4011697E-01,  0.5453512E-01, -0.7660818E-01,  0.1156579    ,
                              -0.2082350    ,  0.7812024    ,  0.4708541    , -0.1786220    ,
                              0.1070651    , -0.7366632E-01,  0.5380140E-01, -0.4038628E-01,
                              0.3063499E-01, -0.2322282E-01,  0.1746474E-01, -0.1292418E-01,
                              0.9351747E-02, -0.6553452E-02,  0.4354787E-02, -0.4491215E-02 },
                              { -0.3610339E-02,  0.4079527E-02, -0.6164601E-02,  0.8939894E-02,
                                -0.1254210E-01,  0.1719870E-01, -0.2319715E-01,  0.3104514E-01,
                                -0.4156235E-01,  0.5635464E-01, -0.7890430E-01,  0.1185305    ,
                                -0.2111745    ,  0.7470909    ,  0.5129383    , -0.1882789    ,
                                0.1118954    , -0.7665114E-01,  0.5581383E-01, -0.4179769E-01,
                                0.3164055E-01, -0.2393967E-01,  0.1797045E-01, -0.1327314E-01,
                                0.9586181E-02, -0.6704286E-02,  0.4451525E-02, -0.4512248E-02 },
                                { -0.3850176E-02,  0.4243306E-02, -0.6407434E-02,  0.9269213E-02,
                                  -0.1297470E-01,  0.1775285E-01, -0.2389400E-01,  0.3190988E-01,
                                  -0.4262333E-01,  0.5764329E-01, -0.8044270E-01,  0.1202436    ,
                                  -0.2120453    ,  0.7112565    ,  0.5544899    , -0.1964535    ,
                                  0.1157289    , -0.7891820E-01,  0.5728737E-01, -0.4279630E-01,
                                  0.3232756E-01, -0.2441106E-01,  0.1828864E-01, -0.1348123E-01,
                                  0.9716844E-02, -0.6781076E-02,  0.4499688E-02, -0.4474290E-02 },
                                  { -0.4063686E-02,  0.4370135E-02, -0.6595247E-02,  0.9517455E-02,
                                    -0.1329186E-01,  0.1814694E-01, -0.2437261E-01,  0.3247959E-01,
                                    -0.4328584E-01,  0.5838729E-01, -0.8121328E-01,  0.1208016    ,
                                    -0.2109199    ,  0.6738845    ,  0.5953013    , -0.2030252    ,
                                    0.1185167    , -0.8044429E-01,  0.5821074E-01, -0.4337718E-01,
                                    0.3269470E-01, -0.2463780E-01,  0.1842126E-01, -0.1355092E-01,
                                    0.9746339E-02, -0.6786256E-02,  0.4500832E-02, -0.4382680E-02 },
                                    { -0.4243565E-02,  0.4456881E-02, -0.6723019E-02,  0.9678334E-02,
                                      -0.1348607E-01,  0.1837216E-01, -0.2462313E-01,  0.3274378E-01,
                                      -0.4353938E-01,  0.5857775E-01, -0.8121268E-01,  0.1202183    ,
                                      -0.2078819    ,  0.6351669    ,  0.6351669    , -0.2078819    ,
                                      0.1202183    , -0.8121268E-01,  0.5857775E-01, -0.4353938E-01,
                                      0.3274378E-01, -0.2462313E-01,  0.1837216E-01, -0.1348607E-01,
                                      0.9678334E-02, -0.6723019E-02,  0.4456881E-02, -0.4243565E-02 } };

  /** Generate the default 28-tap filter table.
      Uses built-in table of 17x28 coefficients for first half of filter
      and expands this to 33x28 filter
   */
  private void defaultFilter () {
    // copy 28*17 HFILTER into 28*33, spaced 64:
    //
    // We start with 17 filters each of length 28.  The aim is to expand this
    // to 2048 filters of length 28.  First we scatter filters 1 thru 17 into
    // the first half of the buffer,  with a gap of 63 empty filter spaces
    // between them, and at the same time we scatter inverted filters into the
    // top half of the buffer, starting at the top and working towards the middle.
    // Note: The code ends up writing to the middle position twice, but that's
    // OK since it is symmetrical.
    final int HALFSIZE = 17;

    filterBuffer = new float[NFILTERS+1][NTAPS];

    for (int ii=0; ii<HALFSIZE; ii++) {
      for (int jj=0; jj<NTAPS; jj++) {
        filterBuffer[ii*64][jj]              = (float)HFILTER[ii][jj]; // work up from the bottom
        filterBuffer[(32-ii)*64][NTAPS-jj-1] = (float)HFILTER[ii][jj]; // work down from the top
      }
    }

    // Now use linear interpolation to create 28*2048 filters
    // Note: there will be one extra filter row at end which will not be used;
    // it is identical to the first filter row.
    for (int ii=0; ii<32; ii++) {
      for (int kk=1; kk<64; kk++) {
        for (int jj=0; jj<NTAPS; jj++) {
          filterBuffer[kk+ii*64][jj] = (float)(((64.0-kk)*(filterBuffer[ ii   *64][jj]) +
              kk*(filterBuffer[(ii+1)*64][jj]))/64.0);
        }
      }
    }
    nTaps = NTAPS;
  }

  /** Read in a user-supplied filter file and expand it into the filterBuffer
      ready for use by applyFilter(). <br>
      <br>
      The input filter file is a type 2000 file consisting of nFilters frames,
      each holding a single interpolating filter of length nTaps (the frame size)
      elements. Each filter exactly spans the delay range from -0.5 to 0.5. We
      interpolate to find filters for as many delays as will fit into the
      available memory.
   */
  private void userFilter () {
    int     nFilters;             // number of filters in filter file
    int     nSpread;              // Spacing factor if filter buffer is big
    int     nGot;                 // Number of frames read from filter file
    int     nElements;            // Number of elements in filter file
    float[] filters;              // filter file data buffer

    nFilters   = (int)filterFile.getSize();
    nTaps      = filterFile.getFrameSize();
    nElements  = nTaps * nFilters;
    nSpread    = (FILTBUFSIZE/filterFPS - nTaps)/(nTaps * (nFilters-1));
    filterData = filterFile.getDataBuffer(nFilters);

    nGot = filterFile.read(filterData, nFilters);
    if (nGot != nFilters) {
      throw new MidasException("Failed to read "+nFilters+" frames from filter file "
          +filterFile.getName());
    }
    filters      = filterData.castF(true);
    filterBuffer = new float[(nFilters-1)*nSpread + 1][nTaps*filterFPS];

    // Copy filters into filterBuffer, spreading if space is available
    for (int ii=0; ii<nFilters; ii++) {
      for (int jj=0; jj<nTaps*filterFPS; jj++) {
        filterBuffer[nSpread*ii][jj] = filters[nTaps*filterFPS*ii + jj];
      }
    }

    if (nSpread > 1) {                  // if spread-out, then interpolate
      float w1, w2;                     // interpolation weights

      for (int ii=0; ii<nFilters-1; ii++) {
        for (int jj=1; jj<nSpread; jj++) {
          w1 = (float)(nSpread-jj)/nSpread;
          w2 = (float)jj/nSpread;
          for (int kk=0; kk<nTaps*filterFPS; kk++) {
            filterBuffer[nSpread*ii+jj][kk] = w1 * filterBuffer[nSpread*ii][kk] +
                w2 * filterBuffer[nSpread*(ii+1)][kk];
          }
        }
      }

      filterData.uncast(filters, false);
    }

    nInterps = (nFilters-1) * nSpread;
    round    = 0.5;
  }

  ////////////////////////////////////////////////////////////////////////////////
  //
  //first_clock
  //
  //Sync file HCB and sync period (pri) passed by /newrate_syninfo/ common
  //block.  Pbufsize may be set to limit sync input for improved operation
  //in a piped macro.
  //Time_start is the earliest time that can be processed from the input
  //file.  Abscissa registration relative to the sync file is assumed.
  //Time_last is the time offset of the last symbol.  It is set to
  //0 by first_clock.
  //Tce is the offset between the sync time and the time for the
  //selected sample from the input file.  It is implemented by
  //adding tce to the start of the sync file.  P_start is, therefore,
  //set equal to the abscissa start of the sync file + tce.
  //Time_now is the time offset of the current sample.
  //Tstep is the offset increment per symbol.  The instantaneous
  //symbol period is tstep+pri.
  //Nstep is the number of periods to be processed with the current
  //value of tstep.
  //Count is both and input and a returned value.  It is initialized
  //by first_clock and used as input by next_clock.  The sync phase
  //file will be read at: time_last + pri*count + p_start.
  //The next sample to be taken frmo the input file is at
  //t_now + pri*count + p_start.
  //
  ////////////////////////////////////////////////////////////////////////////////

  void firstClock (double time_start, double tce ) {

    // Set pbufsize within legal limits
    phaseBufSamples = Math.min( Math.max ( phaseBufSamples, 2 ),PBUFSIZE_MAX );

    // read one buffer at start
    syncData = syncFile.getDataBuffer( phaseBufSamples, FLOAT );
    npts = syncFile.read(syncData, phaseBufSamples );

    if( syncData.getSize() < 1 ) {
      throw new MidasException("Empty sync file.");
    }
    p_start = syncFile.getXStart() + tce;
    bufstart = p_start;
    pp_delta = syncFile.getXDelta();
    //System.out.format( "pp_delta %2.12f pri %2.12f", pp_delta, pri );
    pstep_nominal=pp_delta/pri;
    x_start = time_start;


    // The initial value of count is an estimate for the purpose of
    // reading the sync file.  It may be corrected before exiting.
    //
    outCount=Math.round((x_start - p_start)/pri-0.5);


    if( !matchXM ) {
      // Test to account for file start problem.  Adjust in pri steps at start.
      int nfix = 0;
      if (t_now + outCount*pri+p_start < x_start) {
        nfix = (int)(((x_start-(t_now+(double)(outCount)*pri+(p_start)))/pri)+1);
      }
      else {
        nfix = -(int)(((t_now+outCount*pri+p_start-x_start)/pri));
      }
      //System.out.format( "nfix %d\n", nfix );
      nstep -= nfix;
      outCount += nfix;
    }

    nextClock();

    if( matchXM ) {
      // Test to account for file start problem.  Adjust in pri steps at start.
      int nfix = 0;
      if (t_now + outCount*pri+p_start < x_start) {
        nfix = (int)(((x_start-(t_now+(double)(outCount)*pri+(p_start)))/pri)+1);
      }
      else {
        nfix = -(int)(((t_now+outCount*pri+p_start-x_start)/pri));
      }
      //System.out.format( "nfix %d\n", nfix );
      nstep -= nfix;
      outCount += nfix;
    }

    //x_start += (double)(nfix-1) * inputFile.getXDelta();
    //outputFile.setXStart(x_start);
    /*
    NumberFormat formatter = new DecimalFormat();
    formatter = new DecimalFormat("0.######E0  ");
    System.out.print( formatter.format( t_now ));
    System.out.print( formatter.format( outCount ));
    System.out.print( formatter.format( pri ));
    System.out.print( formatter.format( p_start ));
    System.out.print( formatter.format( x_start ));
    System.out.println("");
    System.out.print( formatter.format(( x_start-(t_now+outCount*pri+p_start))/pri ));
    System.out.print( formatter.format(( t_now+outCount*pri+p_start-x_start)/pri ));
    System.out.println("");
    */
  }

  ////////////////////////////////////////////////////////////////////////////////
  //
  //next_clock
  //
  ////////////////////////////////////////////////////////////////////////////////
  //void nextClock ( double time_last ) {
  void nextClock() {

    double t_last = t_now;
    double time_now = t_last + pri*outCount+p_start;


    // test to see if time_now is within the phase buffer
    long bindex = ( long )((time_now-bufstart)/pp_delta+1)-1;

    // test if we are in leadin extrapolation region
    if (bindex < 0) {
      //adjust time_now so that it is aligned with symbol clock
      t_now = -( syncData.getF( 0 ) % 1.0 ) * pri;

      syncData.getFormat();

      double ioff = Math.round((t_last-t_now)/pri);
      tstep = 0;
      t_now -= ioff*pri;
      nstep = (int) (1.0-outCount-Math.round((t_now)/pri+0.5));

      return;
    }

    // test if we need more data
    if (bindex >= npts-1) {
      if (npts < phaseBufSamples ) {
        // we are past end of phase file so extrapolate
        tstep = 0;
        nstep = Integer.MAX_VALUE;
        t_now = -syncData.getF(npts-1)*pri+bufstart-p_start+pp_delta*(npts-1);
        t_now += Math.round((t_last-t_now)/pri)*pri;
        return;
      }
      double p_index = Math.round((time_now-p_start)/pp_delta+0.5);
      bufstart = p_start+(p_index-1)*pp_delta;
      //read one buffer at start
      npts = syncFile.read(syncData, phaseBufSamples );

      bindex = ( long )((time_now-bufstart)/pp_delta);
    }

    // we are after the first point and before last point so interpolate
    double phase1 = syncData.getF( ( int )bindex );
    double time1=bufstart+bindex*pp_delta;
    double phase2 = syncData.getF( (int)( bindex + 1 ));
    double pdif=phase2-phase1;


    if( nowrap ) {
      //System.out.format("\npdif %2.12f  psn %2.12f \n", pdif, pstep_nominal );
      pdif=pdif+pstep_nominal;
      //System.out.format("pdif %2.12f\n", pdif );
    }
    else {
      pdif=pdif+Math.round(pstep_nominal-pdif);
    }
    double phasint = phase1+pdif*(time_now-time1)/pp_delta;
    // from the interpolated phase adjust time_next to the nearest symbol
//    double temp = (phasint-Math.round(phasint));
//    temp *= (pri); // variable NOT used?


    t_now=t_last-(phasint-Math.round(phasint))*pri;

    double tleft=time1+pp_delta-time_now;
    if (tleft < 0) {
      tleft += pp_delta;
    }
    tstep = pp_delta/pdif-pri;

    nstep= (int)(tleft/pri+1);

  }

  /** Iterate through available input data applying interpolation filters.
      @param nReps  Number of output points (iterations)          [ nreps     ]
      @param inStep Input samples step increment                  [ astep     ]
   */
  private void applyFilter (int nReps, double inStep) {
    int     inOff;      // Integer Offset into input buffer     [ ioff      ]
    int     outOff;     // Offset into output buffer            [ indexb    ]
    double  inReadPos;  // Current read position (double)       [ aindex    ]
    float[] filter;
    long lInBase = inBase;

    outOff = (int)(outWriteP - outBase) * outFPS;

    for (int ii=0; ii<nReps; ii++) {
      // On each interation, calculate the offset into the input buffer by adding
      // the step size to the current input buffer pointer.  The result is not
      // (necessarily) an integer, and we use the integral part to identify
      // the offset of the samples to be interpolated in the input buffer, and
      // the fractional part to index into the table of interpolation filters.
      long ltemp;

      inReadPos = inReadP + inStep*ii;                // Current (double) read position
      ltemp     = (long)inReadPos;
      inOff     = (int)(ltemp - lInBase) * inFPS;  // Offset into input buffer
      //System.out.format( "ii %d  inReadPos %2.12f  inReadP %2.12f  ltemp  %d  inOff %d\n",
      //                    ii, inReadPos, inReadP, ltemp, inOff );

      filter    = filterBuffer[(int)(nInterps * (inReadPos - ltemp) + round)];
      // inBuffer[inOff]   - start of relevant input samples               [sigbuf   ]
      // filter[]          - selected interpolation filter                 [xxbuf5   ]
      // outBuffer[outOff] - output buffer for resulting sample            [xxbuf2   ]

      float rsum = 0.0F, qsum = 0.0F;

      if (cxFilter) {
        if (inFPS == 2) {       // complex input, complex filter, complex output
          for (int jj=0; jj<nTaps; jj++) {
            rsum += inBuffer[inOff+jj*2  ]*filter[jj*2  ] -
                inBuffer[inOff+jj*2+1]*filter[jj*2+1];
            qsum += inBuffer[inOff+jj*2  ]*filter[jj*2+1] +
                inBuffer[inOff+jj*2+1]*filter[jj*2  ];
          }
        } else {              // real input, complex filter, complex output
          for (int jj=0; jj<nTaps; jj++) {
            rsum += inBuffer[inOff+jj]*filter[jj*2  ];
            qsum += inBuffer[inOff+jj]*filter[jj*2+1];
          }
        }
      } else {
        if (inFPS == 2) {       // complex input, real filter, complex output
          for (int jj=0; jj<nTaps; jj++) {
            rsum += inBuffer[inOff+jj*2  ]*filter[jj];
            qsum += inBuffer[inOff+jj*2+1]*filter[jj];
          }
        } else {              // real input, real filter, real output
          for (int jj=0; jj<nTaps; jj++) {
            rsum += inBuffer[inOff + jj] * filter[jj];
          }
        }
      }

      if (tuning) {
        float temp1   = (float)(rsum * rTuneTwiddle - qsum * qTuneTwiddle);
        qsum          = (float)(rsum * qTuneTwiddle + qsum * rTuneTwiddle);
        rsum          = temp1;
        double temp2  = (rTuneTwiddle * rTune - qTuneTwiddle * qTune);
        qTuneTwiddle  = rTuneTwiddle * qTune + qTuneTwiddle * rTune;
        rTuneTwiddle  = temp2;
      }

      outBuffer[outOff++] = rsum;
      if (outFPS == 2) {
        outBuffer[outOff++] = qsum;
      }
    }

    outWriteP += nReps;
  }

  /** Sets the output sample rate
   *  @param rate Desired sample rate in Hz
   */
  public void setRate (double rate) {
    if (rate <= 0.0) {
      throw new MidasException("Invalid new rate [" + rate + "]");
    }
    newRate = rate;
    outputFile.setXDelta(1.0 / newRate);
    inStep  = outputFile.getXDelta() / inputFile.getXDelta();

  }

  /** Shifts output frequency by f if output format is complex
   *  @param tune amount of tune (shift of output frequency), if output format is complex
   */
  public void setTune (double tune) {
    if (outFPS != 2) {
      M.warning("TUNE switch ignored for real output");
      tuning = false;
    }
    else {
      tuning = true;
      rTune = Math.cos(tune * outputFile.getXDelta() * TWOPI);
      qTune = Math.sin(tune * outputFile.getXDelta() * TWOPI);
    }
  }

}
