FFT

FFT.java

public class FFT {
    int n, m;

    // Lookup tables. recompute when size of FFT changes.
    double[] cos;
    double[] sin;

    double[] window;

    public FFT(int n) {
        this.n = n;
        this.m = (int)(Math.log(n) / Math.log(2));

        if (n != (1 << m)) {
            Log.e(TAG, "FFT length must be power of 2");
            return;
        }

        // pre-compute tables
        cos = new double[n/2];
        sin = new double[n/2];

        for (int i = 0; i < n/2; i++) {
            cos[i] = Math.cos(i * -2 * Math.PI / n);
            sin[i] = Math.sin(i * -2 * Math.PI / n);
        }

        makewindow();
    }

    private void makewindow() {
        // blackman window
        // w(n)=0.42-0.5cos{(2*PI*n)/(N-1)}+0.08cos{(4*PI*n)/(N-1)}
        window = new double[n];
        for (int i = 0; i < n; i++) {
            window[i] = 0.42 - 0.5 * Math.cos(i * 2 * Math.PI / (n - 1)) + 0.08 * Math.cos(i * 4 * Math.PI / (n - 1));
        }
    }

    public double[] getWindow() {
        return window;
    }

   /***************************************************************
    * fft.c
    * Douglas L. Jones
    * University of Illinois at Urbana-Champaign
    * January 19, 1992
    * http://cnx.rice.edu/content/m12016/latest/
    *
    *   fft: in-place radix-2 DIT DFT of a complex input
    *
    *   input:
    * n: length of FFT: must be a power of two
    * m: n = 2**m
    *   input/output
    * x: double array of length n with real part of data
    * y: double array of length n with imag part of data
    *
    *   Permission to copy and use this program is granted
    *   as long as this header is included.
    ****************************************************************/

    public void fft(double[] x, double[] y) {
        int i, j, k, n1, n2, a;
        double c, s, e, t1, t2;

        // Bit-reverse
        j = 0;
        n2 = n / 2;
        for (i = 1; i < n - 1; i++) {
            n1 = n2;
            while (j >= n1) {
                j = j - n1;
                n1 = n1 / 2;
            }
            j = j + n1;

            if (i < j) {
                t1 = x[i];
                x[i] = x[j];
                x[j] = t1;
                t1 = y[i];
                y[i] = y[j];
                y[j] = t1;
            }
        }

        // FFT
        n1 = 0;
        n2 = 1;

        for (i = 0; i < m; i++) {
            n1 = n2;
            n2 = n2 + n2;
            a = 0;

            for (j = 0; j < n1; j++) {
                c = cos[a];
                s = sin[a];
                a += 1 << (m - i - 1);
                for (k = j; k < n; k = k + n2) {
                    t1 = c * x[k + n1] - s * y[k + n1];
                    t2 = s * x[k + n1] + c * y[k + n1];
                    x[k + n1] = x[k] - t1;
                    y[k + n1] = y[k] - t2;
                    x[k] = x[k] + t1;
                    y[k] = y[k] + t2;
                }
            }
        }
    }
}