package nxm.dsp.prim;

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

/**
  Perform a FIR filter using FFT convolution algorithm.

  @author Jeff Schoen
*/
public class fdfilter extends Primitive {

  private DataFile hi,hf,ho;
  private Data dbw,dbf,dbm;
  private int nfft,nfir,nbins,flags;
  private int bias,keep,over,pad1,pad2;
  private boolean mux,pad,cxm;
  private Fft fft,ifft;
  private double scale;

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

  public int open() {

    hi = MA.getDataFile ("IN","1000","S#,C#",0);
    hi.open();

    hf = MA.getDataFile ("FIR","1000","S#,C#",0);
    hf.open();

    cxm  = (hi.spa==2) || (hf.spa==2);  // complex mode
    mux  = MA.getState("/MUX");    // mux in&out for hilbert
    pad  = MA.getState("/PAD") || mux;  // pad transient with zeros
    if (mux && cxm) M.error("Cannot /MUX with complex input file");

    nfft = 8192;
    if (MA.getState("/FFT")) {
      nfft = MA.getL("/FFT");
      M.deprecate("FDFILTER: Use of /FFT is deprecated. Use /NFFT instead.");
    }
    nfft = MA.getL("/NFFT",nfft);

    nfir = (int)hf.size;
    over = nfir-1;
    keep = nfft-over;
    bias = (int)Math.round(hf.getXStart()/hf.getXDelta());

    if (pad) {        // pad transient
      pad1 = (over+bias);
      pad2 = -bias;
    } else {
      pad1 = pad2 = 0;
    }

    if (cxm) {
      flags = Fft.COMPLEX;
      scale = 1.0/nfft;
      nbins = nfft;
    } else {
      flags = Fft.REAL|Fft.UNPACKED;
      scale = 1.0/(4*nfft);
      nbins = nfft/2+1;
    }
    fft = new Fft(nfft,flags|Fft.FORWARD);
    ifft = new Fft(nfft,flags|Fft.INVERSE);

    ho = MA.getDataFile ("OUT",hi,0);
    ho.setSize( hi.getSize()-over+pad1+pad2 );
    ho.setXStart( hi.getXStart()+(bias+over-pad1)*hi.getXDelta() );
    if (cxm||mux) ho.setFormatMode(COMPLEX);
    ho.open();

    hi.setXferLength(nfft);
    hi.setConsLength(keep);

    dbw = new Data("CF",nbins);       // data buffer for work
    dbf = new Data("CF",nbins);       // data buffer for filter
    dbm = new Data("CF",nfft);        // data buffer for mux

    return NORMAL;
  }

  public int process() {

    if (hf.avail()>0) loadFilter();    // new filter available

    if (pad1>0) {                      // handle front pad I/O
      DataOp.zero (dbw,0,pad1);        // insert pad1 samples of zero
      dbw.boff = (pad1*dbw.bpe);       // byte offset for read buffer
      hi.setXferLength(nfft-pad1);     // short the transfer
      hi.setConsLength(keep-pad1);     // short the consume
    }

    int n = hi.read(dbw);              // read input data
    if (n==0) return NOOP;
    if (n<over) return FINISH;


    if (pad1>0) {                      // restore std I/O parameters
      dbw.boff = 0;
      hi.setXferLength(nfft);
      hi.setConsLength(keep);
      n += pad1;
      pad1 = Math.max(0,pad1-keep);
    }

    if (n<nfft) {                      // handle back pad I/O
      DataOp.zero (dbw,n,nfft-n);
      n = Math.min(nfft,n+pad2);
      //more pad2 stuff...???
    }
    if (cxm&&hi.spa==1) DataOp.cmplx (dbw,dbw,nfft,DataOp.SC);
    if (mux) DataOp.copy (dbw,-bias*4,4, dbm,0,8, 4,n-over);

    fft.work (dbw);                    // perform forward FFT
    dbw.setSPA(2);
    DataOp.mul (dbw,dbf,dbw,nbins,DataOp.CCC);    // complex multiply
    ifft.work (dbw);                   // perform inverse FFT

    if (mux) {
      dbw.setSPA(1);                   // set for IO libs
      DataOp.copy (dbw,0,4, dbm,4,8, 4,n-over);
      ho.write (dbm,n-over);           // just the good points
    } else {
      dbw.setSPA(ho.spa);              // set for IO libs
      ho.write (dbw,n-over);           // just the good points
    }

    return NORMAL;
  }

  public int close() {
    fft.free();
    ifft.free();
    hi.close();
    hf.close();
    ho.close();
    return NORMAL;
  }

  private void loadFilter() {
    float fscale = (float)scale;
    dbw.setSPA(hf.spa);
    hf.read(dbw,nfir);          // get the FIR taps
    if (MA.find("/SCALE")) {
      float cscale = MA.getF("/SCALE");
      fscale = fscale*cscale;
    }
    Data dbs = new Data(fscale);      // scale factor
    DataOp.mul (dbw,dbs,dbw, nfir);
    dbf.setSPA(hf.spa);
    DataOp.zero (dbf,0,nfft);        // zero
    DataOp.copy (dbw,0, dbf,nfft-over, over);    // put at tail of FFT
    DataOp.copy (dbw,over, dbf,0, 1);      // one point at front
    if (cxm&&hf.spa==1) DataOp.cmplx (dbf,dbf,nfft,DataOp.SC);
    fft.work (dbf);           // perform forward FFT
    dbf.setSPA(2);
  }

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

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

}
