package nxm.dsp.test;

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

/** Raised-Cosine filters are commonly used in digital data communication systems
    to limit Intersymbol Interference (ISI). The Impulse Response of a Raised-Cosine
    filter is zero at each adjacent symbol period.<br>
    Formula is on Wikipedia:<br>
    raisedCos = ((1/T) * sinc((n-center)/T) * cos(PI*rollOff*(n-center)/T) / (1-(2*rollOff*(n-center)/T)<br>
    except where (n-center) = +/- T/(2*rollOff) where<br>
    raisedCos = (PI/(4T)*sinc(1/(2*rollOff))<br>
    Note: T is symbol length<br>
    sinc is normalized sinc (i.e. sin(PI*x)/(PI*x))

    @author  Jeanette McCall
    @since   NeXtMidas 3.9.1
*/
@ProvisionalUseOnly("This has been minimally tested and verified against reliable data. " +
                    "It is currently primarily intended for educational purposes")
public class raisedcos extends Primitive {

  /** Maximum number of taps */ public static final int MAXTAPS=131071;
  /** Default number of taps */ public static final int DEFTAPS = 1025;

  private DataFile ho;
  private Data data;
  private int ntap;          // number of taps (i.e. number of output points)
  private boolean debug;
  private int symbolLength;  // number of samples per symbol (also notated as T)
  private int span;          // number of symbols to span
  private double rollOff;    // roll off (also notated as Beta)
  private boolean normalize; // normalize this so all teh filter taps sum to 1
  /** PI / (symbol length), where T is a common notation for symbol length */
  transient double pi_invT;

  /** deafult constructor */
  public raisedcos () {}

  public int open() {

    debug = MA.getState("/DEBUG");

    // Open output file
    ho = MA.getDataFile("OUT","1000","SF",0);
    ho.open(BaseFile.OUTPUT);

    rollOff = MA.getD("ROLLOFF",MA.getD("BETA"));

    symbolLength = MA.getL("SYMBOLLENGTH",MA.getL("T",1));

    span = MA.getL("SPAN",1);

    ntap = (symbolLength * span) + 1;

    normalize = !MA.isPresent("/NONORM");

    data = new Data("SF",ntap);

    if (debug)  {
      System.out.println("ntap: "+ntap);
      System.out.println("rollOff: "+rollOff);
      System.out.println("symbol length: "+symbolLength);
    }

    pi_invT = Math.PI / symbolLength;
    raisedCos();
    ho.setSize(ntap);
    ho.write (data,ntap);
    ho.close();


    return FINISH;
  }

  /** Raised Cosine - Formula is on Wikipedia:<br>
      raisedCos = ((1/T) * sinc((n-center)/T) * cos(PI*rollOff*(n-center)/T) / (1-(2*rollOff*(n-center)/T)<br>
      except where (n-center) = +/- T/(2*rollOff) where<br>
      raisedCos = (PI/(4T)*sinc(1/(2*rollOff))<br>
      Note: T is symbol length<br>
      sinc is normalised sinc (i.e. sin(PI*x)/(PI*x))
   */
  private void raisedCos () {
    float outbuf[] = data.castF(false);

    // (PI*rollOff)/symbolLength - where T is the common notation for symbolLength
    double rollOffPIinvT = rollOff * pi_invT;
    // 1/symbolLength - where T is the common notation for symbolLength
    double invT = 1.0/symbolLength;
    // ((2*rollOff)/symbolLength)**2 - i.e. (2B/T)**2 where B is one of the chars used for rolloff and T is symbol length
    double twoRollOffInvT_squared = Math.pow(2*rollOff*invT,2);
    // special value needed for +/- T/(2*rollOff)
    double specialCaseValue = specialRaisedCosCase(rollOff, symbolLength);

    int center = (ntap / 2);
    double sum = 0;
    for (int n=0; n < ntap; n++) {
      int shiftedN = n - center;
      // divisor = 1 - (2*rollOff*shiftedN/T)**2
      double divisor = 1-(shiftedN*shiftedN*twoRollOffInvT_squared);
      if (Math.abs(divisor) > 0.000001) {
        // cos(PI*rollOff*shiftedN/T)
        double cosComponent = Math.cos(shiftedN*rollOffPIinvT);
        double sincComponent = sincT(shiftedN);
        outbuf[n]=(float)(invT * sincComponent * cosComponent /divisor);
        if (debug) showAllPartsOfCalc(n, shiftedN, invT, cosComponent, divisor, sincComponent, outbuf[n]);
      }
      else {
        outbuf[n] = (float)(specialCaseValue);
      }
      sum += outbuf[n];
    }
    if (normalize) {
      for (int i=0; i < ntap; i++) {
        outbuf[i] /= sum;
      }
    }

    data.uncast(outbuf, true);
  }

  /**  Normalized sinc(n/T) where T is the symbol length<br>
   *   In out use case it is:<br>
   *   sin((n-center)*PI/T))/((n-center)*PI/T)<br>
   *   This account for the case where n equals 0
   *   @param n tap being calculated
   *   @return normalized sinc for n/T
   */
  private double sincT (int n) {
    if (Math.abs(n/(double)symbolLength) < 0.000001) {return 1.0;}
    double sinc = Math.sin(n*pi_invT) / (n*pi_invT);
    return sinc;
  }

  /**  Normalized sinc(n) <br>
   *   In out use case it is:<br>
   *   This account for the case where n equals 0
   *   @param t value to take sinc of
   *   @return normalized sinc for t
   */
  private static double sinc (double t) {
    if (t<0.0000001) {return 1.0;}
    double sinc = Math.sin(t*Math.PI) / (t*Math.PI);
    return sinc;
  }

  /** (PI/(4T)*sinc(1/(2*rollOff)) for +/- T/(2*rollOff)
   *  @param rollOff roll off
   *  @param T       symbol length
   *  @return        raised cosine special case value
   */
  private static double specialRaisedCosCase (double rollOff, int T) {
    return (Math.PI / (4 * T)) * sinc(1/(2*rollOff));
  }

  private void showAllPartsOfCalc (int n, int shiftedN, double invT, double cosComponent, double divisor, double sincComponent, float tap) {
    System.out.println("n:"+n+" shiftedN:"+shiftedN+" raisedCos = ((1/T) * sinc((n-center)/T) * cos(PI*rollOff*(n-center)/T) / (1-(2*rollOff*(n-center)/T)");
    System.out.println("sinc component:"+sincComponent);
    System.out.println("cos component:"+cosComponent);
    System.out.println("numerator = "+invT+" * "+sincComponent+" * "+cosComponent+" = "+(invT * sincComponent * cosComponent ));
    System.out.println("divisor:"+divisor);
    System.out.println("tap value:"+tap);
  }

}
