package fouriertransform;

/**
 * Fast Fourier Transform
 * @author Winston Prakash
 */
public class FastFourierTransform {
    public static double[] transform(int nData, double[] data, boolean inverse){
        double omega, tempr, tempi, scale;
        double xtemp, cos, sin, xr, xi;
        int j, k, n, M;
        
        //  bit-shift the input data. Remember that N
        //  is the number of samples. There is a real and
        //  a complex component to each data point.
        
        j=0;
        for(int i = 0; i < nData - 1; i++) {
            if (i < j) {
                tempr = data[2 * i];
                tempi = data[2 * i + 1];
                data[2 * i] = data[2 * j];
                data[2 * i + 1] = data[2 * j + 1];
                data[2 * j] = tempr;
                data[2 * j + 1] = tempi;
            }
            k = nData / 2;
            while (k <= j)  {
                j -= k;
                k >>= 1;
            }
            j += k;
        }
        
        //  Perform the FFT. Recursively build up the
        //  full transform from the sub-transforms
        
        if (inverse) {
            scale = -1.0;
        } else {
            scale = 1.0;
        }
        
        M = 2;
        while( M < 2 * nData ) {
            omega = scale * 2.0 * Math.PI / M;
            sin = Math.sin(omega);
            cos = Math.cos(omega) - 1.0;
            xr = 1.0;
            xi = 0.0;
            for (int m = 0; m < M-1; m += 2) {
                for (int i = m; i <2*nData; i += M*2) {
                    j = i + M;
                    tempr = xr * data[j] - xi * data[j + 1];
                    tempi = xr * data[j + 1] + xi * data[j];
                    data[j] = data[i] - tempr;
                    data[j + 1] = data[i + 1] - tempi;
                    data[i] += tempr;
                    data[i + 1] += tempi;
                }
                xtemp = xr;
                xr = xr + xr * cos - xi * sin;
                xi = xi + xtemp * sin + xi * cos;
            }
            M *= 2;
        }
        
        //  If this is a forward transform, multiply
        //  in the 1/N terms
        
        if (!inverse) {
            for(int i = 0; i < nData; i++) {
                data[2 * i] /= nData;
                data[2 * i + 1] /= nData;
            }
        }
        return data;
    }
}