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 divide by N final Complex answer[]=new Complex[data.length]; final double denom=data.length; for(i=0;i0) return i; } } /** * Reverse bits. * @author Don Cross */ private static int reverseBits(int index,final int numBits) { int i,rev; for(i=rev=0;i>=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