package jnt.scimark2;/** Computes FFT's of complex, double precision data where n is an integer power of 2.  * This appears to be slower than the Radix2 method,  * but the code is smaller and simpler, and it requires no extra storage.  * <P>  *  * @author Bruce R. Miller bruce.miller@nist.gov,  * @author Derived from GSL (Gnu Scientific Library),   * @author GSL's FFT Code by Brian Gough bjg@vvv.lanl.gov  */  /* See {@link ComplexDoubleFFT ComplexDoubleFFT} for details of data layout.   */public class FFT {  public static final double num_flops(int N)  {	 double Nd = (double) N;	 double logN = (double) log2(N);	 return (5.0*Nd-2)*logN + 2*(Nd+1);   }  /** Compute Fast Fourier Transform of (complex) data, in place.*/  public static void transform (double data[]) {    transform_internal(data, -1); }  /** Compute Inverse Fast Fourier Transform of (complex) data, in place.*/  public static void inverse (double data[]) {    transform_internal(data, +1);      // Normalize    int nd=data.length;    int n =nd/2;    double norm=1/((double) n);    for(int i=0; i<nd; i++)      data[i] *= norm;  }  /** Accuracy check on FFT of data. Make a copy of data, Compute the FFT, then    * the inverse and compare to the original.  Returns the rms difference.*/  public static double test(double data[]){    int nd = data.length;    // Make duplicate for comparison    double copy[] = new double[nd];    System.arraycopy(data,0,copy,0,nd);    // Transform & invert    transform(data);    inverse(data);    // Compute RMS difference.    double diff = 0.0;    for(int i=0; i<nd; i++) {      double d = data[i]-copy[i];      diff += d*d; }    return Math.sqrt(diff/nd); }  /** Make a random array of n (complex) elements. */  public static double[] makeRandom(int n){    int nd = 2*n;    double data[] = new double[nd];    for(int i=0; i<nd; i++)      data[i]= Math.random();    return data; }  /** Simple Test routine. */  public static void main(String args[]){    if (args.length == 0) {      int n = 1024;      System.out.println("n="+n+" => RMS Error="+test(makeRandom(n))); }    for(int i=0; i<args.length; i++) {      int n = Integer.parseInt(args[i]);      System.out.println("n="+n+" => RMS Error="+test(makeRandom(n))); }  }  /* ______________________________________________________________________ */  protected static int log2 (int n){    int log = 0;    for(int k=1; k < n; k *= 2, log++);    if (n != (1 << log))      throw new Error("FFT: Data length is not a power of 2!: "+n);    return log; }  protected static void transform_internal (double data[], int direction) {	if (data.length == 0) return;    	int n = data.length/2;    if (n == 1) return;         // Identity operation!    int logn = log2(n);    /* bit reverse the input data for decimation in time algorithm */    bitreverse(data) ;    /* apply fft recursion */	/* this loop executed log2(N) times */    for (int bit = 0, dual = 1; bit < logn; bit++, dual *= 2) {      double w_real = 1.0;      double w_imag = 0.0;      double theta = 2.0 * direction * Math.PI / (2.0 * (double) dual);      double s = Math.sin(theta);      double t = Math.sin(theta / 2.0);      double s2 = 2.0 * t * t;      /* a = 0 */      for (int b = 0; b < n; b += 2 * dual) {        int i = 2*b ;        int j = 2*(b + dual);        double wd_real = data[j] ;        double wd_imag = data[j+1] ;                  data[j]   = data[i]   - wd_real;        data[j+1] = data[i+1] - wd_imag;        data[i]  += wd_real;        data[i+1]+= wd_imag;      }            /* a = 1 .. (dual-1) */      for (int a = 1; a < dual; a++) {        /* trignometric recurrence for w-> exp(i theta) w */        {          double tmp_real = w_real - s * w_imag - s2 * w_real;          double tmp_imag = w_imag + s * w_real - s2 * w_imag;          w_real = tmp_real;          w_imag = tmp_imag;        }        for (int b = 0; b < n; b += 2 * dual) {          int i = 2*(b + a);          int j = 2*(b + a + dual);          double z1_real = data[j];          double z1_imag = data[j+1];                        double wd_real = w_real * z1_real - w_imag * z1_imag;          double wd_imag = w_real * z1_imag + w_imag * z1_real;          data[j]   = data[i]   - wd_real;          data[j+1] = data[i+1] - wd_imag;          data[i]  += wd_real;          data[i+1]+= wd_imag;        }      }    }  }  protected static void bitreverse(double data[]) {    /* This is the Goldrader bit-reversal algorithm */    int n=data.length/2;	int nm1 = n-1;	int i=0; 	int j=0;    for (; i < nm1; i++) {      //int ii = 2*i;      int ii = i << 1;      //int jj = 2*j;      int jj = j << 1;      //int k = n / 2 ;      int k = n >> 1;      if (i < j) {        double tmp_real    = data[ii];        double tmp_imag    = data[ii+1];        data[ii]   = data[jj];        data[ii+1] = data[jj+1];        data[jj]   = tmp_real;        data[jj+1] = tmp_imag; }      while (k <= j) 	  {        //j = j - k ;		j -= k;        //k = k / 2 ;         k >>= 1 ; 	  }      j += k ;    }  }}