InplaceFFT.java

Below is the syntax highlighted version of InplaceFFT.java from §9.7 Optimization.

/******************************************************************************
 *  Compilation:  javac InplaceFFT.java
 *  Execution:    java InplaceFFT n
 *  Dependencies: Complex.java
 *
 *  Compute the FFT of a length n complex sequence in-place.
 *  Uses a non-recursive version of the Cooley-Tukey FFT.
 *  Runs in O(n log n) time.
 *
 *  Reference:  Algorithms 1.5.2 and 1.6.1 in Computational Frameworks
 *              for the Fast Fourier Transform by Charles Van Loan.
 *
 *
 *  Limitations
 *  -----------
 *   -  assumes n is a power of 2
 *
 *  
 ******************************************************************************/

public class InplaceFFT {

    // compute the FFT of x[]
    // precondition: the length of x[] is a power of 2
    public static void fft(Complex[] x) {

        // check that length is a power of 2
        int n = x.length;
        if (Integer.highestOneBit(n) != n) {
            throw new IllegalArgumentException("n is not a power of 2");
        }

        // bit reversal permutation
        int shift = 1 + Integer.numberOfLeadingZeros(n);
        for (int k = 0; k < n; k++) {
            int j = Integer.reverse(k) >>> shift;
            if (j > k) {
                Complex temp = x[j];
                x[j] = x[k];
                x[k] = temp;
            }
        }

        // butterfly updates
        for (int L = 2; L <= n; L = L+L) {
            for (int j = 0; j < L/2; j++) {
                double jth = 2 * Math.PI * j / L;
                Complex w = new Complex(Math.cos(jth), -Math.sin(jth));
                for (int k = 0; k < n/L; k++) {
                    Complex tao = w.times(x[k*L + j + L/2]);
                    x[k*L + j + L/2] = x[k*L + j].minus(tao); 
                    x[k*L + j]       = x[k*L + j].plus(tao); 
                }
            }
        }
    }


    // test client
    public static void main(String[] args) { 
        int n = Integer.parseInt(args[0]);
        Complex[] x = new Complex[n];

        // original data
        for (int i = 0; i < n; i++) {
            x[i] = new Complex(i, 0);
            // x[i] = new Complex(-2*Math.random() + 1, 0);
        }
        for (int i = 0; i < n; i++)
            StdOut.println(x[i]);
        StdOut.println();

        // FFT of original data
        fft(x);
        for (int i = 0; i < n; i++)
            StdOut.println(x[i]);
        StdOut.println();
    }

}