package fr.ensea.math;
// was jsci.maths ?
import com.imsl.math.*;

/**
 * The Fourier math library.
 * This class cannot be subclassed or instantiated because all methods are static.
 * Adapted from JSci (www.jsci.org)
 *
 * Normalization rules :
 *
 * X_q = sum_n x_n e^(-I q x_n) avec x_n = n a (a = pas du reseau, 1 ici) => pas de 1/sqrt(N)
 * 
 * x_n = 1/N * sum_q X_q e^(I q x_n) avec q = m q_0 et q_0 = (2 pi) / (N a), afin d'assurer les PBC.
 *
 * @version 0.7
 * @author Mark Hale, Sylvain Reynal
 */
public final class FFT  {

    public final static double TWO_PI=6.2831853071795864769252867665590057683943387987502;
 	
	/**
	* Fourier transform for Complex data
	* @return an array containing positive frequencies in ascending order
	* followed by negative frequencies in ascending order.
	* @author Don Cross
	*/
	public static Complex[] transform(final Complex data[]) {
		final double arrayRe[]=new double[data.length];
		final double arrayIm[]=new double[data.length];
		int i,j,k,n;
		int numBits,blockSize,blockEnd;
		double deltaAngle,angleRe,angleIm,tmpRe,tmpIm;
		double alpha,beta; // used in recurrence relation

		if(!isPowerOf2(data.length))
			throw new IllegalArgumentException("The number of samples must be a power of 2.");

		numBits=numberOfBitsNeeded(data.length);
		// Simultaneous data copy and bit-reversal ordering into output
		for(i=0;i<data.length;i++) {
			j=reverseBits(i,numBits);
			arrayRe[j]=data[i].real();
			arrayIm[j]=data[i].imag();
		}
		// FFT
		blockEnd=1;
		for(blockSize=2;blockSize<=data.length;blockSize<<=1) {
			deltaAngle=TWO_PI/blockSize;
			alpha=Math.sin(0.5*deltaAngle);
			alpha=2.0*alpha*alpha;
			beta=Math.sin(deltaAngle);
			for(i=0;i<data.length;i+=blockSize) {
				angleRe=1.0;angleIm=0.0;
				for(j=i,n=0;n<blockEnd;j++,n++) {
					k=j+blockEnd;
					tmpRe=angleRe*arrayRe[k]-angleIm*arrayIm[k];
					tmpIm=angleRe*arrayIm[k]+angleIm*arrayRe[k];
					arrayRe[k]=arrayRe[j]-tmpRe;
					arrayIm[k]=arrayIm[j]-tmpIm;
					arrayRe[j]+=tmpRe;arrayIm[j]+=tmpIm;
					tmpRe=alpha*angleRe+beta*angleIm;
					tmpIm=alpha*angleIm-beta*angleRe;
					angleRe-=tmpRe;angleIm-=tmpIm;
				}
			}
			blockEnd=blockSize;
		}
		final Complex answer[]=new Complex[data.length];
		for(i=0;i<data.length;i++)
			answer[i]=new Complex(arrayRe[i],arrayIm[i]);
		return answer;
	}
	/**
	* Fourier transform for real data
	* @return an array containing positive frequencies in ascending order
	* followed by negative frequencies in ascending order.
	*/
	public static Complex[] transform(final double data[]) {
		final double arrayRe[]=new double[data.length];
		final double arrayIm[]=new double[data.length];
		int i,j,k,n;
		int numBits,blockSize,blockEnd;
		double deltaAngle,angleRe,angleIm,tmpRe,tmpIm;
		double alpha,beta; // used in recurrence relation

		if(!isPowerOf2(data.length))
			throw new IllegalArgumentException("The number of samples must be a power of 2.");

		numBits=numberOfBitsNeeded(data.length);
		// Simultaneous data copy and bit-reversal ordering into output
		for(i=0;i<data.length;i++) {
			j=reverseBits(i,numBits);
			arrayRe[j]=data[i];
			arrayIm[j]=0.0;
		}
		// FFT
		blockEnd=1;
		for(blockSize=2;blockSize<=data.length;blockSize<<=1) {
			deltaAngle=TWO_PI/blockSize;
			alpha=Math.sin(0.5*deltaAngle);
			alpha=2.0*alpha*alpha;
			beta=Math.sin(deltaAngle);
			for(i=0;i<data.length;i+=blockSize) {
				angleRe=1.0;angleIm=0.0;
				for(j=i,n=0;n<blockEnd;j++,n++) {
					k=j+blockEnd;
					tmpRe=angleRe*arrayRe[k]-angleIm*arrayIm[k];
					tmpIm=angleRe*arrayIm[k]+angleIm*arrayRe[k];
					arrayRe[k]=arrayRe[j]-tmpRe;
					arrayIm[k]=arrayIm[j]-tmpIm;
					arrayRe[j]+=tmpRe;arrayIm[j]+=tmpIm;
					tmpRe=alpha*angleRe+beta*angleIm;
					tmpIm=alpha*angleIm-beta*angleRe;
					angleRe-=tmpRe;angleIm-=tmpIm;
				}
			}
			blockEnd=blockSize;
		}
		final Complex answer[]=new Complex[data.length];
		for(i=0;i<data.length;i++)
			answer[i]=new Complex(arrayRe[i],arrayIm[i]);
		return answer;
	}
	/**
	* Inverse Fourier transform for complex data
	* @return an array containing the positive time part of the signal
	* followed by the negative time part.
	* @author Don Cross
	*/
	public static Complex[] inverseTransform(final Complex data[]) {
		final double arrayRe[]=new double[data.length];
		final double arrayIm[]=new double[data.length];
		int i,j,k,n;
		int numBits,blockSize,blockEnd;
		double deltaAngle,angleRe,angleIm,tmpRe,tmpIm;
		double alpha,beta; // used in recurrence relation

		if(!isPowerOf2(data.length))
			throw new IllegalArgumentException("Data length must be a power of 2.");

		numBits=numberOfBitsNeeded(data.length);
		// Simultaneous data copy and bit-reversal ordering into output
		for(i=0;i<data.length;i++) {
			j=reverseBits(i,numBits);
			arrayRe[j]=data[i].real();
			arrayIm[j]=data[i].imag();
		}
		// FFT
		blockEnd=1;
		for(blockSize=2;blockSize<=data.length;blockSize<<=1) {
			deltaAngle=-TWO_PI/blockSize;
			alpha=Math.sin(0.5*deltaAngle);
			alpha=2.0*alpha*alpha;
			beta=Math.sin(deltaAngle);
			for(i=0;i<data.length;i+=blockSize) {
				angleRe=1.0;angleIm=0.0;
				for(j=i,n=0;n<blockEnd;j++,n++) {
					k=j+blockEnd;
					tmpRe=angleRe*arrayRe[k]-angleIm*arrayIm[k];
					tmpIm=angleRe*arrayIm[k]+angleIm*arrayRe[k];
					arrayRe[k]=arrayRe[j]-tmpRe;
					arrayIm[k]=arrayIm[j]-tmpIm;
					arrayRe[j]+=tmpRe;arrayIm[j]+=tmpIm;
					tmpRe=alpha*angleRe+beta*angleIm;
					tmpIm=alpha*angleIm-beta*angleRe;
					angleRe-=tmpRe;angleIm-=tmpIm;
				}
			}
			blockEnd=blockSize;
		}
		// Normalize -> divide by N
		final Complex answer[]=new Complex[data.length];
		final double denom=data.length;
		for(i=0;i<data.length;i++)
			answer[i]=new Complex(arrayRe[i]/denom,arrayIm[i]/denom);
		return answer;
	}
	/**
	* Inverse Fourier transform for real datas.
	* @return an array containing the positive time part of the signal
	* followed by the negative time part.
	*/
	public static Complex[] inverseTransform(final double data[]) {
		final double arrayRe[]=new double[data.length];
		final double arrayIm[]=new double[data.length];
		int i,j,k,n;
		int numBits,blockSize,blockEnd;
		double deltaAngle,angleRe,angleIm,tmpRe,tmpIm;
		double alpha,beta; // used in recurrence relation

		if(!isPowerOf2(data.length))
			throw new IllegalArgumentException("Data length must be a power of 2.");

		numBits=numberOfBitsNeeded(data.length);
		// Simultaneous data copy and bit-reversal ordering into output
		for(i=0;i<data.length;i++) {
			j=reverseBits(i,numBits);
			arrayRe[j]=data[i];
			arrayIm[j]=0.0;
		}
		// FFT
		blockEnd=1;
		for(blockSize=2;blockSize<=data.length;blockSize<<=1) {
			deltaAngle=-TWO_PI/blockSize;
			alpha=Math.sin(0.5*deltaAngle);
			alpha=2.0*alpha*alpha;
			beta=Math.sin(deltaAngle);
			for(i=0;i<data.length;i+=blockSize) {
				angleRe=1.0;angleIm=0.0;
				for(j=i,n=0;n<blockEnd;j++,n++) {
					k=j+blockEnd;
					tmpRe=angleRe*arrayRe[k]-angleIm*arrayIm[k];
					tmpIm=angleRe*arrayIm[k]+angleIm*arrayRe[k];
					arrayRe[k]=arrayRe[j]-tmpRe;
					arrayIm[k]=arrayIm[j]-tmpIm;
					arrayRe[j]+=tmpRe;arrayIm[j]+=tmpIm;
					tmpRe=alpha*angleRe+beta*angleIm;
					tmpIm=alpha*angleIm-beta*angleRe;
					angleRe-=tmpRe;angleIm-=tmpIm;
				}
			}
			blockEnd=blockSize;
		}
		// Normalize
		final Complex answer[]=new Complex[data.length];
		final double denom=data.length;
		for(i=0;i<data.length;i++)
			answer[i]=new Complex(arrayRe[i]/denom,arrayIm[i]/denom);
		return answer;
	}


	/**
	* Returns true if x is a power of 2.
	* @author Don Cross
	*/
	private static boolean isPowerOf2(final int x) {
		final int BITS_PER_WORD=32;
		for(int i=1,y=2;i<BITS_PER_WORD;i++,y<<=1) {
			if(x==y)
				return true;
		}
		return false;
	}
	/**
	* Number of bits needed.
	* @author Don Cross
	*/
	private static int numberOfBitsNeeded(final int pwrOf2) {
		if(pwrOf2<2)
			throw new IllegalArgumentException();
		for(int i=0;;i++) {
			if((pwrOf2&(1<<i))>0)
				return i;
		}
	}
	/**
	* Reverse bits.
	* @author Don Cross
	*/
	private static int reverseBits(int index,final int numBits) {
		int i,rev;
		for(i=rev=0;i<numBits;i++) {
			rev=(rev<<1)|(index&1);
			index>>=1;
		}
		return rev;
	}

	/**
	* Sorts the output from the Fourier transfom methods into
	* ascending frequency/time order.
	*/
	public static Complex[] sort(final Complex output[]) {
		final Complex ret[]=new Complex[output.length];
		final int Nby2=output.length/2;
		for(int i=0;i<Nby2;i++) {
			ret[Nby2+i]=output[i];
			ret[i]=output[Nby2+i];
		}
		return ret;
	}


	private FFT() {}
}

