package cquest.utils;
//------------------------------------------------------------------------
//                           Fft.java
//------------------------------------------------------------------------
//
// HISTORY :
//------------------------------------------------------------------------
//|Date    |Author     |Method/Description                               |
//------------------------------------------------------------------------
//|10/05/98   |C.Couturier|creation                                         |
//------------------------------------------------------------------------
//
import cquest.*;

/** this class allows to carry out some fourier transforms or reverse fourier transforms
@author 	C. Couturier
@version 	%I%, %U%*/
public class Fft {

	  
  final int REAL = 0;
  final int IMAG = 1;

	/** constructor
	@author 	C. Couturier
	@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 	C. Couturier
	@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 	C. Couturier
	@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 	C. Couturier
	@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 	C. Couturier
	@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 	C. Couturier
	@version 	%I%, %U%*/
	public synchronized double[][] direct_1D_Estim_fft(double[][] input)
		{return n_Direct_1D_Estim_fft(input);}

/** 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 	C. Couturier
	@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 	C. Couturier
	@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 	C. Couturier
	@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 	C. Couturier
	@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[0][j]=input1[0][j];
			input[1][j]=input1[1][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 	C. Couturier
	@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[0][i] = Math.cos((double)i*theta);
       omega[1][i] = Math.sin((double)i*theta);
    }
    double y[][] = calcFFTD(x,omega);
    if (!forward)
      for(int i=0; i<n; i++)
      {
        y[0][i]/=n;
        y[1][i]/=n;
      }
    return y;
   }

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

      shift= (omega[0].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[0][k] = x[0][k<<1];
       z[1][k] = x[1][k<<1];
       w[0][k] = x[0][(k<<1)+1];
       w[1][k] = x[1][(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[0][k] = z[0][k];
				y[1][k] = z[1][k];
				tr = omega[0][shifted]*w[0][k] - omega[1][shifted]*w[1][k];
				w[1][k] = omega[1][shifted]*w[0][k] + omega[0][shifted]*w[1][k];
				w[0][k] = tr;
				y[0][k] += w[0][k];
				y[1][k] += w[1][k];
				y[0][k+(n>>1)] = z[0][k];
				y[1][k+(n>>1)] = z[1][k];
				y[0][k+(n>>1)] -= w[0][k];
				y[1][k+(n>>1)] -= w[1][k];
				shifted += shift;
			}
		return y;
	}
//------------------------------------------------------------------------
// shiftFft()
//------------------------------------------------------------------------
	/** to
		@author 	D.Stefan
		@version 	%I%, %U%*/

	public void shiftFft(double [] signal) {

	}
//------------------------------------------------------------------------
// shiftFft2D()
//------------------------------------------------------------------------
	/** to
		@author 	D.Stefan
		@version 	%I%, %U%*/

	public double[][] shiftFft2D(double[][] matrix) {
		int width = matrix.length;
		int height = matrix[0].length;
		double [][] temp = new double [width][height];
		for (int i=0; i<width; i++) {
			for (int j=0; j<height; j++) {
				temp[i][j]=matrix[((i+(width/2))%width)][((j+(height/2))%height)];
			}
		}
		return temp;
	}

}