package nxm.dsp.prim;

import nxm.sys.lib.BaseFile;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.Parser;
import nxm.sys.lib.Primitive;

/**
  Creates finite impulse response (FIR) filter using the FIR Kaiser method.

  @author Rich Holley / Jeff Schoen
  @since  NeXtMidas 2.7.1
*/
public class firkais extends Primitive {

  /** List of filter types */
  public static final String typeList = "LowPass,HighPass,BandPass";
  /** filter type constants (see {@link #typeList} for details). */
  public static final int LOW=1, HIGH=2, BAND=3;

  private DataFile ho;
  private Data dbo;
  private int type,ntap;
  private double freq1,freq2,freq3,freqf,attn,delta;
  private float scale;
  private boolean cx,dbg;
  private int maxTaps = 32768;

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

  public int open () {

    // Get parameters
    type  = MA.getSelectionIndex("TYPE",typeList,0);
    freqf = MA.getD("FOLDING");
    freq1 = MA.getD("FREQUENCY")/freqf;
    freq2 = MA.getD("BANDWIDTH")/freqf;
    freq3 = MA.getD("TRANSITION",freqf*.01)/freqf;
    attn  = MA.getD("ATTENUATION");
    ntap  = MA.getL("NTAP");
    ntap  = MA.getL("/NTAPS",ntap); // old switch syntax
    cx    = MA.getState("/CMPLX");
    dbg   = MA.getState("/DEBUG");
    scale = MA.getF("/SCALE",1.0F);
    delta = 0.5/freqf;

    // Open output file
    ho = MA.getDataFile("OUT","1000",cx?"CF":"SF",0);
    ho.setSize(ntap);
    ho.setXDelta(delta);
    ho.open(BaseFile.OUTPUT);

    // get data buffers
    dbo = ho.getDataBuffer(maxTaps*2,FLOAT);
    float[] fbo = dbo.castF(false);

    // generate FIR taps
    ntap = generate(type,freq1,freq2,freq3,attn,fbo,ntap, dbg);

    // scale and convert to real
    if (!cx && type==BAND) scale *= 2;
    if (cx) for (int i=0; i<ntap*2; i++)  fbo[i] *= scale;
    else    for (int i=0; i<ntap; i++) fbo[i] = fbo[i+i]*scale;

    // write data buffers
    dbo.uncast(fbo,true);
    ho.write(dbo,ntap);
    ho.setXStart(-.5*(ntap-1)*delta);
    ho.close();

    MR.put(MA.getU("/TAPRSLT"),ntap);

    return FINISH;
  }

  /** generate FIR taps
   *  @param type filter type
   *  @param f1 Lowpass cutoff / Highpass Cutin / Bandpass center freq.(Hz)
   *  @param f2 Filter bandwidth (Hz)
   *  @param f3 Filter cutoff transition bandwidth (Hz)
   *  @param attn Out-of-band attenuation (dB)
   *  @param filt Data buffer for filter
   *  @param ntap Number of taps (-1 for minimum)
   *  @param dbg  debug flag
   *  @return number of taps
   */
  public static int generate (int type, double f1, double f2, double f3,
                double attn, float[] filt, int ntap, boolean dbg) {

    double fcntr = f1;
    double fpass = (type==BAND)? f2*.5 : f1;
    double ftnbw = f3;
    double fsamp = 2.0; // all values previously normalized by nyquist freq

    // compute number of taps required for filter response
    double mtap = (attn-8)/(2.285*2*PI*ftnbw/fsamp);
    if (ntap <= 0) ntap = (int)(mtap+1.5);
    if (ntap <= 0) return 0;

    // compute "BETA" factor for KAISER window
    double beta;
    if (attn > 50)
      beta = 0.1102*(attn-8.7);
    else if (attn >= 21)
      beta = 0.58420*Math.pow(attn-21.0,0.4) + 0.07886*(attn-21.0);
    else
      beta = 0.0;

    // compute alpha
    double alpha = 0.5 * (ntap-1);

    // compute radian cutoff frequency
    double fcut = (type==HIGH)? fpass - .5*ftnbw : fpass + .5*ftnbw;
           fcut /= fsamp;
    double wcut = 2*PI*fcut;

    // compute constant denominator term
    double den = bessel_io(beta);

    // setup BPF local oscillator stuff
    if (type!=BAND) fcntr = 0;
    double strt = -.5*(ntap-1);
    double lo   = 2*PI*fcntr/fsamp;

    // print debug info
    if (dbg) {
      StringBuilder sb = new StringBuilder(300); // usually around 265-300 chars
      sb.append("Type    : "+type+ " ("+Parser.get(typeList, type)+ ")\n");
      sb.append("Ntap    : "+ntap+"\n");
      sb.append("Freq1   : "+f1+"\n");
      sb.append("Freq2   : "+f2+"\n");
      sb.append("Freq3   : "+f3+"\n");
      sb.append("Ftnbw   : "+ftnbw+"\n");
      sb.append("Attn    : "+attn+"\n");
      sb.append("Fpass   : "+fpass+"\n");
      sb.append("Beta    : "+beta+"\n");
      sb.append("Alpha   : "+alpha+"\n");
      sb.append("Fcut    : "+fcut+"\n");
      sb.append("Wcut    : "+wcut+"\n");
      sb.append("Den     : "+den+"\n");
      sb.append("Strt    : "+strt+"\n");
      sb.append("Lo      : "+lo+"\n");
      System.out.print(sb);
    }

    // compute lowpass filter response
    for (int i=0; i<ntap; i++) {
      double pcos = Math.cos((strt+i)*lo);
      double psin = Math.sin((strt+i)*lo);
      double resp = 2*fcut*sinc(wcut*(i-alpha));
      if (type==HIGH) resp = sinc(PI*(i-alpha)) - resp;
      double factor = Math.sqrt(1.0-Math.pow((i-alpha)/alpha,2.0));
      double wind = bessel_io(beta*factor) / den;
      filt[i+i+0] = (float)(wind*resp*pcos);
      filt[i+i+1] = (float)(wind*resp*psin);
    }

    return ntap;
  }

  private static double sinc (double arg) {
    double num = Math.sin(arg);
    if (Math.abs(arg)<1e-20 && Math.abs(num)<1e-20) return 1.0;
    return num/arg;
  }

  private static double bessel_io (double x) {
    double ds=1, d=0, s=1;
    for (;;) {
      d += 2;
      ds = (ds*x*x)/(d*d);
      s += ds;
      if (ds-(s*0.2e-8) < 0) break;
    }
    return s;
  }

}
