package nxm.dsp.lib;

import nxm.sys.inc.InternalUseOnly;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataOp;
import nxm.sys.lib.Midas;
import nxm.sys.lib.MidasException;
import nxm.sys.lib.Shell;
import nxm.sys.libm.Fft;

/** A Frequency Impulse Response (FIR) filter. This replaces code that was internal
    to TUNER (and FILTER), and FIRWIND.

    @author  J. Czechowski IV (based on code from Jason Broussard, Jeff Schoen, and Ken Smith)
    @since   NeXtMidas 3.3.1
 */
@InternalUseOnly("since initial version")
public class FilterFIR {
  /** Filter type. */
  public static enum FilterType {
    /** Low pass filter   */ LowPass,
    /** High pass filter  */ HighPass,
    /** Bandpass filter  */ BandPass,
    /** Stopband filter */ StopBand
    };

  private final Midas      midas;      // Midas context to use.
  private       FilterType type;       // Filter type.
  private       int        ntaps;      // Number of taps.
  private       double     freq1;      // 1st transition frequency (Hz).
  private       double     freq2;      // 2nd transition frequency (Hz).
  private       double     freqF;      // Folding frequency or nyquist (Hz).
  private       boolean    complex;    // Complex filter?
  private       boolean    symmetric;  // Symmetric filter?
  private       String     windowName; // Window name (or "NONE").
  private       double     scaleFact;  // Scale factor (usually 1.0).

  /** Creates a new FIR filter.
      @param midas      The midas context to use.
      @param type       Filter type.
      @param ntaps      Number of taps.
      @param freq1      1st transition frequency (Hz).
      @param freq2      2nd transition frequency (Hz).
      @param freqF      Folding frequency or nyquist (Hz).
      @param complex    Complex filter?
      @param symmetric  Symmetric filter?
      @param windowName Window name (or "NONE").
      @param scaleFact  Scale factor (usually 1.0).
  */
  public FilterFIR (Midas midas, FilterType type, int ntaps, double freq1, double freq2,
                     double freqF, boolean complex, boolean symmetric, String windowName,
                     double scaleFact) {
    this.midas      = (midas != null)? midas : Shell.getSharedMidasContext();
    this.type       = type;
    this.ntaps      = ntaps;
    this.freq1      = freq1;
    this.freq2      = freq2;
    this.freqF      = freqF;
    this.complex    = complex;
    this.symmetric  = symmetric;
    this.windowName = (windowName == null)? "NONE" : windowName;
    this.scaleFact  = scaleFact;
    validate();
  }

  /** Validate the various settings.
     <pre>
    ==== Frequency Checks For Complex Data ====================================
       -                                  +
       n           f            f         n
       q           r            r         q
       u           e            e         u
       i           q            q         i
       s           1            2         s
       t           .            .         t
       |           .            .         |
       |           .            .         |
    -------------------------------------------------------------------------
       |           .            .         |
       |           .            .         |           HighPass
       |           ******************************
       |           *            .         |
    ****************            .         |
       |           .            .         |
    -------------------------------------------------------------------------
       |           .            .         |
       |           .            .         |           LowPass
    ****************            .         |
       |           *            .         |
       |           ******************************
       |           .            .         |
    -----------------------------------------------------------------------
       |           .            .         |
       |           .            .         |           BandPass
       |           **************         |
       |           *            *         |
    ****************            *****************
       |           .            .         |
    -------------------------------------------------------------------------
       |           .            .         |
       |           .            .         |           StopBand
    ****************            *****************
       |           *            *         |
       |           **************         |
    -------------------------------------------------------------------------
     </pre>
   */
  private void validate () {
    // ==== FILTER TYPE ========================================================
    if (type == null) throw new MidasException("Unknown filter type: null");

    // ==== FREQUENCIES ========================================================
    //
    // Validate that freq1 is in range (for LOWPASS & HIGHPASS there will not be any freq2)
    //
    if (complex) {
      if (freq1 < -freqF)
        throw new MidasException("Transition criteria (freq1 >= -freqF) not met, "
                                +"given freq1="+freq1+", freqF="+freqF);
      if (freq1 > +freqF)
        throw new MidasException("Transition criteria (freq1 <= +freqF) not met, "
                                +"given freq1="+freq1+", freqF="+freqF);
    }
    else {
      if (freq1 < 0.0  )
        throw new MidasException("Transition criteria (freq1 >= 0.0) not met, "
                                +"given freq1="+freq1);
      if (freq1 > freqF)
        throw new MidasException("Transition criteria (freq1 <= freqF) not met, "
                                +"given freq1="+freq1+", freqF="+freqF);
    }

    //
    // For filters that have a second transition "freq2" (BANDPASS & STOPBAND filters),
    // make sure freq 2 is in range
    //
    if ((type == FilterType.BandPass) || (type == FilterType.StopBand)) {
      if (freq2 < freq1)
        throw new MidasException("Transition criteria (freq2 >= freq1) not met, "
                                +"given freq1="+freq1+", freq2="+freq2);
      if (freq2 > freqF)
        throw new MidasException("Transition criteria (freq1 <= freqF) not met, "
                                +"given freq2="+freq2+", freqF="+freqF);
    }
  }

  /** Creates a new BandPass FIR filter.
      @param midas      The midas context to use.
      @param ntaps      Number of taps.
      @param cfreq      The requested IF (in HZ).
      @param bw         The requested bandwidth (in Hz).
      @param nyquist    The nyquist of the input signal.
      @param dec        Decimation (1 = no decimation).
      @param complexOut output complex (by adjusting fmin to -nyquist, and hence max bandwidth)
      @param complex    Complex filter?
      @param symmetric  Symmetric filter?
      @param windowName Window name (or "NONE").
      @param scaleFact  Scale factor (usually 1.0).
      @return The new filter.
   */
  @InternalUseOnly // used only by TUNER
  public static FilterFIR initBandPass (Midas midas, int ntaps, double cfreq, double bw,
                                        double nyquist, int dec, boolean complexOut,
                                        boolean complex, boolean symmetric,
                                        String windowName, double scaleFact) {
    if (midas == null) midas = Shell.getSharedMidasContext();

    double fmax  = nyquist;
    double fmin  = (complexOut)? -nyquist : 0;
    double maxBW = fmax - fmin;

    if (bw > maxBW/dec) {
      midas.warning("Bandwidth limited to output bandwidth");
      bw = maxBW/dec;
    }

    double freq1 = cfreq - bw/2;
    double freq2 = cfreq + bw/2;
    double freqF = nyquist;

    // Make sure that the filter lies within the input bandwidth

    if (freq1 < fmin) {
      freq1 = fmin;
      freq2 = fmin+bw;
      midas.warning("Filter adjusted to fit in available input bandwidth "+freq1+" to "+freq2);
    }
    if (freq2 > fmax) {
      freq1 = fmax-bw;
      freq2 = fmax;
      midas.warning("Filter adjusted to fit in available input bandwidth "+freq1+" to "+freq2);
    }

    return new FilterFIR(midas, FilterType.BandPass, ntaps, freq1, freq2, freqF,
                         complex, symmetric, windowName, scaleFact);
  }

  @Override
  public String toString () {
    return "FilterFIR: FilterType="+type
                   + " NTaps="+ntaps
                   + " Freq1="+freq1+" Hz"
                   + " Freq2="+freq2+" Hz"
                   + " FoldingFreq="+freqF+" Hz"
                   + " Complex="+complex
                   + " Symmetric="+symmetric
                   + " WindowName="+windowName
                   + " ScaleFactor="+scaleFact;
  }

  //////////////////////////////////////////////////////////////////////////////
  // FILTER CREATION METHODS
  //////////////////////////////////////////////////////////////////////////////
  /** Creates the filter by calling {@link #makeFilter(int)} with min length of 0.
      @return The filter
   */
  public float[] makeFilter () {
    return makeFilter(0);
  }

  /** Creates the filter with option to override length of output array.
      @param minLen The minimum length of the output array (0=Default). This
                    can be used in cases where the output filter will be
                    immediately applied to other computations which require
                    a longer array (e.g. FFT overlap, up-sampling, etc)
                    without necessitating an extra array-copy to a larger array.
      @return The filter.
   */
  public float[] makeFilter (int minLen) {
    Data    dbw  = DataOp.getWindow(windowName, ntaps, scaleFact/ntaps, midas);
    float[] wind = dbw.castF(true);
    int     off  = (symmetric)? -2      : 0;
    int     len  = (symmetric)? ntaps*2 : (ntaps+1)*2;     // *2 for (real,imag)
    float[] filt = new float[Math.max(minLen, len)];
    double  f1   = freq1 / freqF;
    double  f2   = freq2 / freqF;
    int     flags; // inverse FFT flags (set below)

    // convert freqs to bins of interest
    int bin1 = ((int)(0.5*ntaps*f1)+1)*2; // = ((int)floor(0.5*ntaps*f1)+1)*2
    int bin2;
    if ((type==FilterType.BandPass) || (type==FilterType.StopBand)) {
      if (f2==0) f2=1.00; // f2=100%
      bin2 = ((int)(0.5*ntaps*f2)+1)*2;
    }
    else {
      bin2 = ntaps+2;
    }
    bin1 += off;
    bin2 += off;
//System.out.println("    bin1 = "+bin1+" bin2 = "+bin2);

    while (complex && (bin1 < 0)) {
      int _bin1 = bin1;
      int _bin2 = bin2;

      switch (type) {
        case LowPass:  bin1 = len-_bin1;                     type = FilterType.HighPass; break;
        case HighPass: bin1 = len-_bin1;                     type = FilterType.LowPass;  break;
        case BandPass: bin1 = _bin2; bin2 = (ntaps*2)+_bin1; type = FilterType.StopBand; break;
        case StopBand: bin1 = _bin2; bin2 = (ntaps*2)+_bin1; type = FilterType.BandPass; break;
      }
//System.out.println("==> bin1 = "+bin1+" bin2 = "+bin2+" type = "+type);
    }

    if ((bin1 < 0) || (bin1 > bin2) || (bin2 > len)) {
      throw new MidasException("Can not create "+type+" filter for freq1="+freq1+" Hz "
                              +"freq2="+freq2+" Hz freqF="+freqF+" Hz.");
    }

    // make ones in pass bands
    if (complex) {
      flags = Fft.INVERSE|Fft.COMPLEX;

      int n = ntaps*2;
      switch (type) {
        case LowPass:  for (int i=0;    i<bin1; i+=2) { filt[i]=1; filt[n-i]=1; } break;
        case HighPass: for (int i=bin1; i<bin2; i+=2) { filt[i]=1; filt[n-i]=1; } break;
        case BandPass: for (int i=bin1; i<bin2; i+=2) { filt[i]=1; } break;
        case StopBand: for (int i=0;    i<bin1; i+=2) { filt[i]=1; }
                       for (int i=bin2; i<n;    i+=2) { filt[i]=1; } break;
      }
    }
    else {
      flags = Fft.INVERSE|Fft.REAL|Fft.UNPACKED;

      int n = ntaps+1;
      switch (type) {
        case LowPass:  for (int i=0;    i<bin1; i+=2) { filt[i]=1; } break;
        case HighPass: for (int i=bin1; i<bin2; i+=2) { filt[i]=1; } break;
        case BandPass: for (int i=bin1; i<bin2; i+=2) { filt[i]=1; } break;
        case StopBand: for (int i=0;    i<bin1; i+=2) { filt[i]=1; }
                       for (int i=bin2; i<n;    i+=2) { filt[i]=1; } break;
      }
    }

    // shift alternate points by 180 degrees to make inverse FFT non-causal
    for (int i=2; i<len; i+=4) filt[i] = -filt[i];

//Integer _bin1a = null;
//Integer _bin1b = null;
//Integer _bin2a = null;
//Integer _bin2b = null;
//int mode = 0;
//for (int i = 0; i < filt.length; i+=2) {
//  //System.out.println(i+":\t"+buf[i]+"\t"+buf[i+1]+"\t=>\t"+filtBuf[i]+"\t"+filtBuf[i+1]+"\t=\t"+bufx[i]+"\t"+bufx[i+1]);
//  switch (mode) {
//    case 0: if (filt[i] != 0) { _bin1a = i; mode++; } break;
//    case 1: if (filt[i] == 0) { _bin1b = i; mode++; } break;
//    case 2: if (filt[i] != 0) { _bin2a = i; mode++; } break;
//    case 3: if (filt[i] == 0) { _bin2b = i; mode++; } break;
//  }
//}
//System.out.println("_bin1a = "+_bin1a+" _bin1b = "+_bin1b+" _bin2a = "+_bin2a+" _bin2b = "+_bin2b);

    // inverse FFT
    Fft.process(filt,ntaps,flags);

    // multiply by window
    if (complex) for (int i=0,j=0; i<ntaps; i++) { filt[j++] *= wind[i]; filt[j++] *= wind[i]; }
    else         for (int i=0,j=0; i<ntaps; i++) { filt[j++] *= wind[i]; }

    dbw.uncast(wind,false);
    return filt;
  }

  //////////////////////////////////////////////////////////////////////////////
  // GET/SET METHODS
  //////////////////////////////////////////////////////////////////////////////
  /** Get filter type
   * @return The filter type. */
  public FilterType getFilterType  () { return type; }
  /** Get Number of taps
   *  @return Number of taps. */
  public int        getNTaps       () { return ntaps; }
  /** Get 1st transition freq (Hz
   *  @return The 1st transition freq (Hz). */
  public double     getFreq1       () { return freq1; }
  /** Get 2nd transition freq (Hz)
   *  @return The 2nd transition freq (Hz). */
  public double     getFreq2       () { return freq2; }
  /** Get Folding freq or nyquist (Hz)
   *  @return Folding freq or nyquist (Hz). */
  public double     getFoldingFreq () { return freqF; }
  /** Get Complex filter?
   *  @return Complex filter? */
  public boolean    getComplex     () { return complex; }
  /** Get Symmetric filter?
   *  @return Symmetric filter? */
  public boolean    getSymmetric   () { return symmetric; }
  /** Get window name (or "NONE")
   * @return The window name (or "NONE"). */
  public String     getWindowName  () { return windowName; }
  /** Get scale factor (usually 1)
   *  @return The scale factor (usually 1). */
  public double     getScaleFactor () { return scaleFact; }

  /** Set filter type
   *  @param t  The filter type. */
  public void setFilterType  (FilterType t ) { type       = t; }
  /** Set Number of taps
   *  @param n  Number of taps. */
  public void setNTaps       (int        n ) { ntaps      = n; }
  /** Set 1st transition freq (Hz)
   *  @param f  The 1st transition freq (Hz). */
  public void setFreq1       (double     f ) { freq1      = f; }
  /** Set 2nd transition freq (Hz)
   *  @param f  The 2nd transition freq (Hz). */
  public void setFreq2       (double     f ) { freq2      = f; }
  /** Set Folding freq or nyquist (Hz)
   *  @param f  Folding freq or nyquist (Hz). */
  public void setFoldingFreq (double     f ) { freqF      = f; }
  /** Set Complex filter?
   *  @param cx Complex filter? */
  public void setComplex     (boolean    cx) { complex    = cx; }
  /** Set Symmetric filter?
   *  @param s  Symmetric filter? */
  public void setSymmetric   (boolean    s ) { symmetric  = s; }
  /** Set window name (or "NONE")
   *  @param wn The window name (or "NONE"). */
  public void setWindowName  (String     wn) { windowName = (wn == null)? "NONE" : wn; }
  /** Set scale factor (usually 1)
   *  @param sf The scale factor (usually 1). */
  public void setScaleFactor (double     sf) { scaleFact  = sf; }


  /** Gets the bandwidth of the filter. Only applies to a {@link FilterType#BandPass}
      (or {@link FilterType#StopBand}) filter, all others return NaN.
      @return the bandwidth of this filter
   */
  public double getBandwidth () {
    switch (type) {
      case LowPass:  return Double.NaN;
      case HighPass: return Double.NaN;
      case BandPass: return (freq2 - freq1) / 2;
      case StopBand: return (freq2 - freq1) / 2;
      default:       throw new AssertionError("Unknown type "+type);
    }
  }

  /** Gets the center frequency of the filter. This is equal to:
      <code>(freq1 + freq2) / 2</code>
      @return the center frequency of this filter
   */
  public double getCenterFreq () {
    return (freq1 + freq2) / 2;
  }

  /** Filter delay. This uses the following equation for computing delay in a
      linear-phase filter:
     <pre>
               (ntaps - 2)     (ntaps - 2)
       delay = -----------  =  ----------- * xdelta
                2.0 * Fs           2.0
     </pre>
     @param xdelta The sampling interval (xdelta from the filter file).
     @return the filter delay
   */
  public double getDelay (double xdelta) {
    return (ntaps - 2) / 2.0 * xdelta;
  }

  /** Sets the freq1, freq2, freqF, and bandwidth.
      @param freq1      1st transition frequency (Hz).
      @param freq2      2nd transition frequency (Hz).
      @param freqF      Folding frequency or nyquist (Hz).
      @throws MidasException if the combination of freq1, freq2, and freqF are invalid.
   */
  public void setFreq (double freq1, double freq2, double freqF) {
    this.freq1 = freq1;
    this.freq2 = freq2;
    this.freqF = freqF;
    validate();
  }
}
