package fft;


//                           Fft.java
//------------------------------------------------------------------------
//

/** this class allows to carry out some fourier transforms or reverse fourier transforms
@author 	H.Ratiney
@version 	%I%, %U%*/
public class Fft {
  final static int REAL = 0;
  final static int IMAG = 1;

  /** constructor
  @author 	H. Ratiney
  @version 	%I%, %U%*/
  public Fft() {
  }

  static {
    System.loadLibrary("fftw");
  }

  /** fourier transform (FORWARD and 1D) provided in the package 'fftw' (library : fftw).
  	'measure' means that FFTW actually runs and measures the execution time of several
  	FFTs in order to find the best way to compute the transform of size n. This may
  	take some time, depending on your installation and on the precision of the timer
  	in your machine
  	@param 		FID[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FFT [2 (REAL/IMAG)][length], length can be whatever >0
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  private synchronized native double[][] n_Direct_1D_Measure_fft(double[][] valuesJ);

  /** fourier transform (FORWARD and 1D) provided in the package 'fftw' (library : fftw).
  	FFTW_ESTIMATE does not run any computation, and just builds a reasonable plan,
  	which may be sub-optimal
  	@param 		FID[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FFT [2 (REAL/IMAG)][length], length can be whatever >0
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  private synchronized native double[][] n_Direct_1D_Estim_fft(double[][] valuesJ);

  /** REVERSE fourier transform (BACKWARD and 1D) provided in the package 'fftw' (library : fftw).
  	FFTW_ESTIMATE does not run any computation, and just builds a reasonable plan,
  	which may be sub-optimal
  	@param 		FFT[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FIF [2 (REAL/IMAG)][length], length can be whatever >0
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  private synchronized native double[][] n_Reverse_1D_Estim_fft(double[][] valuesJ);

  /** REVERSE fourier transform (BACKWARD and 1D) provided in the package 'fftw' (library : fftw).
  	'measure' means that FFTW actually runs and measures the execution time of several
  	FFTs in order to find the best way to compute the transform of size n. This may
  	take some time, depending on your installation and on the precision of the timer
  	in your machine
  	@param 		FFT[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FID
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  private synchronized native double[][] n_Reverse_1D_Measure_fft(double[][] valuesJ);

  /** fourier transform (FORWARD and 1D) provided in the package 'fftw' (library : fftw).
  	FFTW_ESTIMATE does not run any computation, and just builds a reasonable plan,
  	which may be sub-optimal
  	@param 		FID[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FFT [2 (REAL/IMAG)][length], length can be whatever >0
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  public synchronized double[][] direct_1D_Estim_fft(double[][] input) {

    // Correction to correct lifting of spectrum due to DFT:
    // 1st datapoint should be average between 1st and last datapoint.
    // See e.g. page 132 in E. Oran Brigham, 1974, The Fast Fourier Transform.
    // Changes made by Yoeri on 12-4-2002

    double[][] daReturn = null;

    if (input != null) {
      if (input[REAL].length > 0) {
        double[] daFirst = new double[2];
        daFirst[REAL] = input[REAL][0];
        daFirst[IMAG] = input[IMAG][0];

        int iLast = input[REAL].length - 1;
        //input[REAL][0] = (input[REAL][0] + input[REAL][iLast]) / 2;
        //input[IMAG][0] = (input[IMAG][0] + input[IMAG][iLast]) / 2;

        try {
          daReturn = n_Direct_1D_Estim_fft(input);
        }
        catch (Exception e) {
          daReturn = calcFFTD(input, true); // This backup was only in FIDResult, put here by Yoeri.
        }
        finally {
          input[REAL][0] = daFirst[REAL];
          input[IMAG][0] = daFirst[IMAG];
        }
      }
    }

    if (daReturn == null)
      System.out.println("Error in FFT");

    return daReturn;

    //	   return n_Direct_1D_Estim_fft(input);      // Old method, without averaging.
  }

  /** fourier transform (FORWARD and 1D) provided in the package 'fftw' (library : fftw).
  	'measure' means that FFTW actually runs and measures the execution time of several
  	FFTs in order to find the best way to compute the transform of size n. This may
  	take some time, depending on your installation and on the precision of the timer
  	in your machine
  	@param 		FID[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FFT [2 (REAL/IMAG)][length], length can be whatever >0
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  public synchronized double[][] direct_1D_Measure_fft(double[][] input) {
    return n_Direct_1D_Measure_fft(input);
  }

  /** REVERSE fourier transform (BACKWARD and 1D) provided in the package 'fftw' (library : fftw).
  	FFTW_ESTIMATE does not run any computation, and just builds a reasonable plan,
  	which may be sub-optimal
  	@param 		FFT[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FIF [2 (REAL/IMAG)][length], length can be whatever >0
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  public synchronized double[][] reverse_1D_Estim_fft(double[][] input) {
    return n_Reverse_1D_Estim_fft(range(input));
  }

  /** REVERSE fourier transform (BACKWARD and 1D) provided in the package 'fftw'
  	'measure' means that FFTW actually runs and measures the execution time of several
  	FFTs in order to find the best way to compute the transform of size n. This may
  	take some time, depending on your installation and on the precision of the timer
  	in your machine
  	@param 		FFT[2 (REAL/IMAG)][length], length can be whatever >0
  	@return 	FID
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  public synchronized double[][] reverse_1D_Measure_fft(double[][] input) {
    return n_Reverse_1D_Measure_fft(range(input));
  }

  /** FFTW computes an unnormalized transform, that is, the equation IFFT(FFT(X)) = n X holds.
  	In other words, applying the forward and then the backward transform will multiply the
  	input by n.
  	An FFTW_FORWARD transform corresponds to a sign of -1 in the exponent of the DFT. Note
  	also that we use the standard "in-order" output ordering--the k-th output corresponds to
  	the frequency k/n (or k/T, where T is your total sampling period). For those who like to
  	think in terms of positive and negative frequencies, this means that the positive
  	frequencies are stored in the first half of the output and the negative frequencies are
  	stored in backwards order in the second half of the output. (The frequency -k/n is the
  	same as the frequency (n-k)/n.)
  	@param 		fft
  	@return 	signal (in the right direction)
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  private synchronized double[][] range(double[][] input1) {
    int len = input1[0].length;
    int m = (int) Math.floor(len / 2 + 0.5);
    double[][] input = new double[2][len];
    double rpart, ipart;
    int j;
    for (j = 0; j < len; j++) {
      input[REAL][j] = input1[REAL][j];
      input[IMAG][j] = input1[IMAG][j];
    }

    for (int i = 0; i < m; i++) {
      rpart = input[REAL][i];
      ipart = input[IMAG][i];
      j = i + m;
      if (j < len) {
        input[REAL][i] = input[REAL][j] / len;
        input[IMAG][i] = input[IMAG][j] / len;
        input[REAL][j] = rpart / len;
        input[IMAG][j] = ipart / len;
      }
      else {
        input[REAL][i] = input[REAL][i] / len;
        input[IMAG][i] = input[IMAG][i] / len;
      }
    }
    return input;
  }

  /** fourier transform (reverse or direct) for signals of a power 2 length
  (less efficient that the other methods (cf. fftw)
  	@param 		FFT or FID [2 (REAL/IMAG)][length], length is a power of 2
  	@param 		true means forward, otherwise backward
  	@return 	FFT or FID [2 (REAL/IMAG)][length], length is a power of 2
  	@author 	H. Ratiney
  	@version 	%I%, %U%*/
  public synchronized double[][] calcFFTD(double[][] x, boolean forward) {
    int n = x[0].length;
    double omega[][] = new double[2][n >> 1];
    double theta = -2 * Math.PI / (double) n;
    if (!forward)
      theta = -theta;
    for (int i = 0; i < (n >> 1); i++) {
      omega[REAL][i] = Math.cos((double) i * theta);
      omega[IMAG][i] = Math.sin((double) i * theta);
    }
    double y[][] = calcFFTD(x, omega);
    if (!forward)
      for (int i = 0; i < n; i++) {
        y[REAL][i] /= n;
        y[IMAG][i] /= n;
      }
    return y;
  }

  /** method used by calcFFTD */
  private synchronized static double[][] calcFFTD(double[][] x, double[][] omega) {
    int n = x[REAL].length, k = 0, shift, shifted = 0;
    if (n == 1) {
      double z1[][] = new double[2][1];
      z1[REAL][0] = x[REAL][0];
      z1[IMAG][0] = x[IMAG][0];
      return z1;
    }

    shift = (omega[REAL].length / (n >> 1));

    double z[][] = new double[2][n >> 1];
    double w[][] = new double[2][n >> 1];

    for (k = 0; k < n >> 1; k++) {
      z[REAL][k] = x[REAL][k << 1];
      z[IMAG][k] = x[IMAG][k << 1];
      w[REAL][k] = x[REAL][(k << 1) + 1];
      w[IMAG][k] = x[IMAG][(k << 1) + 1];
    }

    z = calcFFTD(z, omega);
    w = calcFFTD(w, omega);

    double y[][] = new double[2][n];
    double tr;
    for (k = 0; k < n >> 1; k++) {
      y[REAL][k] = z[REAL][k];
      y[IMAG][k] = z[IMAG][k];
      tr = omega[REAL][shifted] * w[REAL][k] - omega[IMAG][shifted] * w[IMAG][k];
      w[IMAG][k] = omega[IMAG][shifted] * w[REAL][k] + omega[REAL][shifted] * w[IMAG][k];
      w[REAL][k] = tr;
      y[REAL][k] += w[REAL][k];
      y[IMAG][k] += w[IMAG][k];
      y[REAL][k + (n >> 1)] = z[REAL][k];
      y[IMAG][k + (n >> 1)] = z[IMAG][k];
      y[REAL][k + (n >> 1)] -= w[REAL][k];
      y[IMAG][k + (n >> 1)] -= w[IMAG][k];
      shifted += shift;
    }
    return y;
  }
}