package nxm.dsp.prim;

import nxm.sys.inc.DataTypes;
import nxm.sys.lib.BaseFile;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.Primitive;

/**
  Creates FIR filter using Kaiser method<br>
  Note: This class moved from the UCL option tree into the DSP option tree
  in NeXtMidas 3.1.0.

  @author  Jason Broussard
  @since   NeXtMidas 2.1.1
*/
public class kaiswind extends Primitive {

  /** List of filter types */
  public static final String typeList = "LowPass,HighPass,BandPass";
  /** Lowpass filter         */ public static final int LOW     =      1;
  /** Highpass filter        */ public static final int HIGH    =      2;
  /** Bandpass filter        */ public static final int BAND    =      3;
  /** Maximum number of taps */ public static final int MAXTAPS = 131072;

  private DataFile ho,hi;
  private Data dbw;
  private int type,ntap,intap,maxtaps=131072;
  boolean debug, inref, cmplx;
  private double freq1,freq2,freq3,fold,attndb,tnbw,fcntr,fpass;
  private static final int REAL = 0;
  private static final int IMAG = 1;
  private String taprslt;

  /** deault constructor */
  public kaiswind () {}

  public int open() {

    debug = MA.getState("/DEBUG");
    if (debug) System.out.println("Running kaiswin with db="+attndb+" and ntap="+ntap);
    cmplx = MA.getState("/CMPLX",false);
    intap  = MA.getL("/NTAPS",maxtaps);
    if (intap > maxtaps) {
      M.warning("Maximum number of taps exceeded. Using default maximum of 131072");
      intap = maxtaps;
    }

    taprslt = MA.getU("/TAPRSLT");
    // Get parameters
    type  = MA.getSelectionIndex("TYPE",typeList,1);
    if (type<=0) M.error("Unknown filter type: "+MA.getU("TYPE"));
    freq1 = MA.getD("FREQ");
    freq2 = MA.getD("BW");
    freq3 = MA.getD("TRANS",-1.0);//transition bandwidth
    if (MA.getState("INPUTREF")) {
      hi = MA.getDataFile (MA.getU("INPUTREF"),"1000,2000","S#,C#",0);
      hi.open();
      fold = 0.5*(1.0/hi.getXDelta());
      hi.close();
    } else {
      fold = MA.getD("FOLD");
      if (fold <= 0.0) M.error("Must specify folding frequency");
    }
    attndb = MA.getD("ATTEN",60.0);

    // error check
//    if (freq1<0 || freq1>1) M.error("Transition freq1 criteria (0<=freq<=nyquist) not met");
//    if (freq2<0 || freq2>1) M.error("Transition freq2 criteria (0<=freq<=nyquist) not met");
    if (type==BAND) {
      fcntr = freq1;
      fpass = 0.5*freq2;
    }
    else {
      fcntr = 0.0;
      fpass = freq1;
    }
    if (freq3 == -1.0) tnbw = 0.01*fold;
    else	       tnbw = freq3;

    // Open output file
    ho = MA.getDataFile("OUT","1000","SF",0);
    if (cmplx) {
      ho.setFormatMode(COMPLEX);
      dbw = new Data("CF",maxtaps);
    }
    else{
      dbw = new Data("SF",maxtaps);
    }
    ho.setXStart ((double)(intap-1)*(-0.25/fold));
    ho.setXDelta (0.5/fold);
    ho.setSize(intap);
    //ho.switch=optional ??
    ho.open(BaseFile.OUTPUT);

    if (intap != maxtaps)
      ntap = -1*intap;
    else
      ntap = intap;
    if (debug)  {
      System.out.println("intap: "+intap);
      System.out.println("ntap: "+ntap);
      System.out.println("type: "+type);
      System.out.println("freq1: "+freq1);
      System.out.println("freq2: "+freq2);
      System.out.println("freq3: "+freq3);
      System.out.println("attndb: "+attndb);
      System.out.println("fcntr: "+fcntr);
      System.out.println("fpass: "+fpass);
      System.out.println("tnbw: "+tnbw);
    }

    kaiswin(dbw, ntap);
    ho.setXStart ((double)(ntap-1)*(-0.25/fold));
    ho.setXDelta (0.5/fold);
    ho.setSize(ntap);
    ho.write (dbw,ntap);
    ho.close();
    //write number of taps to result
    MR.put(taprslt, ntap);

    return FINISH;
  }

  private void kaiswin (Data odata, int ntaps) {
    float[][] buf = new float[maxtaps][2];
    float outbuf[] = odata.castF(false);
    double beta,bessio_of_beta,bess_arg_top,fsamp=2.0*fold,fcut,wcut,strt,lo;
    int i,mtaps;

    mtaps = ntaps;
    if (fsamp == 0.0 || tnbw == 0.0) {
      ntaps = 0;
    }
    else {
      float tmpr8 = (float)((attndb-8)/(2.285*2.0*Math.PI*tnbw/fsamp)+0.5) + 1.0f;
      //check tmpr8 for largeness, warn if maxint is exceeded
      ntaps = (int)tmpr8;
      if (ntaps > mtaps) {
        ntaps = 0;
      }
      else if (tmpr8 > DataTypes.LMAX) {
        M.error("Unreasonable input parameters given, please reconsider.");
      }
    }

    if (mtaps < 0) ntaps = -1*mtaps;
    if (type==HIGH) {
      if (ntaps%2 == 0) ntaps = ntaps + 1;
      if (mtaps < 0) ntaps = -1*mtaps;
    }
    if (attndb > 50.0)
      beta = 0.1102*(attndb-8.7);
    else if (attndb >= 21 && attndb <= 50)
      beta = 0.58420*Math.pow((attndb-21.0),0.4) + 0.07886*(attndb-21.0);
    else
      beta = 0.0;

    double dalpha = 0.5 * (ntaps - 1.0);
    //Compute radian cutoff freq
    if (ntaps > 0) {
      if (type==HIGH)
        fcut = (fpass-0.5*tnbw)/fsamp;
      else
        fcut = (fpass+0.5*tnbw)/fsamp;
      wcut = 2.0*Math.PI*fcut;
    }
    else {
      fcut = 0.0;
      wcut = 0.0;
    }

    //Compute constant denominator term
    bessio_of_beta = window_bessio_d(beta);
    if (bessio_of_beta <= 0.0) bessio_of_beta = 1e-20; //????
    if (type!=BAND) fcntr = 0.0;
    strt = -0.5*(ntaps-1) - 1;
    lo = 2.0*Math.PI*(fcntr/fsamp);
    if (debug) {
      System.out.println("fsamp: "+fsamp);
      System.out.println("beta: "+beta);
      System.out.println("dalpha: "+dalpha);
      System.out.println("fcut: "+fcut);
      System.out.println("wcut: "+wcut);
      System.out.println("bessio_of_beta: "+bessio_of_beta);
      System.out.println("strt: "+strt);
      System.out.println("lo: "+lo);
      System.out.println("ultimately ntaps: "+ntaps);
    }
    int jj=0;
    for (i=0;i<ntaps;i++) {
      double dresp, dwind;

      double dpcos = Math.cos(lo*(strt+(i+1)));
      double dpsin = Math.sin(lo*(strt+(i+1)));
      if (type==HIGH) {
        dresp = firwindow_sinc(Math.PI*(i-dalpha))-2.0*fcut*firwindow_sinc(wcut*(i-dalpha));
      }
      else {
        dresp = 2.0*fcut*firwindow_sinc(wcut*(i-dalpha));
      }
      bess_arg_top = beta * Math.sqrt(1.0 - Math.pow(((i-dalpha)/dalpha),2.0));
      dwind = (float)(window_bessio_d(bess_arg_top)/bessio_of_beta);
      buf[i][REAL] = (float)(dwind*dresp*dpcos);
      buf[i][IMAG] = (float)(dwind*dresp*dpsin);
      if (cmplx) {
        outbuf[jj++] = buf[i][REAL];
        outbuf[jj++] = buf[i][IMAG];
      }
      else {
        float scale = 1;
        if (type==BAND) scale = 2;
        outbuf[i] = buf[i][REAL]*scale;
      }
    }
    odata.uncast(outbuf, true);
    ntap = ntaps;
  }

  /** Double precision Modified Zeroth order bessel function routine for calculating Kaiser window */
  private double window_bessio_d (double x) {
    double t = 1.e-12; // Once new iteration contributes less than this fraction of the sum to the sum, exit
    double y = 0.5*x;
    double e = 1.0;
    double de = 1.0;
    double xi,sde;
    for (int i=1; i<26; i++) {
      xi = i;
      de *= y/xi;
      sde = de*de;
      e += sde;
      if ((e*t-sde) > 0) break;
    }
    return e;
  }

  private double firwindow_sinc (double darg) {
    double dnum,dsinc;

    dnum = Math.sin(darg);
    if (darg < 1e-20 && darg > -1e-20 &&
        dnum < 1e-20 && dnum > -1e-20)
      dsinc = 1.0;
    else
      dsinc = dnum / darg;
    return dsinc;
  }

}
