/***************************************************************************
 // @(#)vsFFT1D.java
 //
 // Perform a non recursive in place FFT using doubles.
 //
 // @version 1.0, 4 JULY 1997
 // @author  D. Lyon, V. Silva
 ****************************************************************************/
package transforms;

/**
 * A 1D FFT for integral power of two FFTs
 */
public class FFT1d {
    private final int FORWARD_FFT = -1;
    private final int REVERSE_FFT = 1;
    private int direction = FORWARD_FFT;

    private static final float twoPI = (float) (2 * Math.PI);

    // These arrays hold the complex 1D FFT data.
    private double realPart[] = null;
    private double imaginaryPart[] = null;

    public static final double gauss(
            double x,
            double xc, double sigma) {
        double oneOnSigmaSquaredOn2 = 1 / (sigma * sigma) / 2;
        return
                Math.exp(-((x - xc) * (x - xc)) *
                oneOnSigmaSquaredOn2) / Math.PI * oneOnSigmaSquaredOn2;
    }
    public static final double[]getGauss(int n){
        double g[] = new double[n];
        for (int x=0; x < n; x++)
            g[x]=gauss(x,n/2,2);
        return g;
    }
    public static final void testGauss() {
        print(getGauss(10));
    }
    public static final void print(double d[]){
       for (int i=0; i < d.length; i++)
           System.out.println(d[i]);
    }

    /**
     * A way to visually test the 1D FFT on a small amount of data.
     */
    public void visuallyTest() {
        System.out.println("Testing FFT engine in 1D class...");
        int N = 8;


        double[] in_r = new double[N];
        double[] in_i = new double[N];

        double[] holder_r = new double[N];
        double[] holder_i = new double[N];

        // Create a test ramp signal.
        for (int i = 0; i < N; i++) {
            in_r[i] = (double) i / N;
            in_i[i] = (double) 0;
        }

        forwardFFT(in_r, in_i);

        // Copy to new array because IFFT will destroy the FFT results.
        for (int i = 0; i < N; i++) {
            holder_r[i] = in_r[i];
            holder_i[i] = in_i[i];
        }

        for (int i = 0; i < N; i++) {
            in_r[i] = in_r[i];
            in_i[i] = in_i[i];
        }

        reverseFFT(in_r, in_i);
        double[] mag = magnitudeSpectrum(in_r, in_i);

        System.out.println("i    x1[i]    re[i]             im[i]           tv[i]");
        for (int i = 0; i < N; i++) {
            System.out.println(i +
                    "\t" + i +
                    "\t" + holder_r[i] +
                    "\t\t" +
                    holder_i[i] +
                    "\t\t" + mag[i]);
        }


    }

    public void timeFFT(int n) {
        // Time a 256K FFT
        double realPart[] = new double[n];
        double imaginaryPart[] = new double[n];

        // Create a test ramp signal.

        synthesizeRamp(realPart, imaginaryPart);

        ip.gui.Timer t1 = new ip.gui.Timer();
        runExperiment(t1, realPart, imaginaryPart);


    }

    /**
     * Destroy the input data with a linear ramp.
     * @param realPart input datas real component
     * @param imaginaryPart  input datas' imaginary component.
     */
    public static final void synthesizeRamp(double[] realPart, double[] imaginaryPart) {
        int n = realPart.length;
        for (int i = 0; i < n; i++) {
            realPart[i] = (double) i / n;
            imaginaryPart[i] = (double) 0;
        }
    }

    private void runExperiment(ip.gui.Timer t1,
                               double[] realPart,
                               double[] imaginaryPart) {
        t1.start();
        double w[] = getGauss(realPart.length);

        for (int i = 0; i < 1000; i++) {
            window(realPart,w);
            forwardFFT(realPart, imaginaryPart);
            getMaxPSDLocation(realPart, imaginaryPart);
        }
        // Stop the timer and report.
        t1.print(realPart.length +
                " done 1000 times point 1D FFT (using double):");
    }
    private static final void window(double r[], double w[]){
          for (int i=0; i < r.length; i++)
              r[i]=r[i]*w[i];
    }

    /**
     * 1D FFT utility functions.
     */
    public void swap(int i, int numBits) {
        double tempr;

        int j = bitr(i, numBits);

        tempr = realPart[j];
        realPart[j] = realPart[i];
        realPart[i] = tempr;
        tempr = imaginaryPart[j];
        imaginaryPart[j] = imaginaryPart[i];
        imaginaryPart[i] = tempr;
    }

    // swap Zi with Zj
    private void swapInt(int i, int j) {
        double tempr;
        int ti;
        int tj;
        ti = i - 1;
        tj = j - 1;
        // System.out.print("swapInt " + ti+","+tj+"\t");
        // System.out.println(Integer.toBinaryString(ti)+","+Integer.toBinaryString(tj));
        // System.out.println(Integer.toBinaryString(ti)+"bitr="+Integer.toBinaryString(bitr(ti)));
        tempr = realPart[tj];
        realPart[tj] = realPart[ti];
        realPart[ti] = tempr;
        tempr = imaginaryPart[tj];
        imaginaryPart[tj] = imaginaryPart[ti];
        imaginaryPart[ti] = tempr;
    }

    /**
     * get the location of the maximum partial
     * @param realPart
     * @param imaginaryPart
     * @return  location of max(realPart[i],imaginaryPart[i])
     */
    public static final int getMaxPSDLocation(
            double realPart[],
            double imaginaryPart[]) {
        double max;
        max = -0.99e30;
        int location = 0;
        for (int i = 0; i < realPart.length; i++)
            if (magnitude(realPart[i], imaginaryPart[i]) > max) {
                max = realPart[i];
                location = i;
            }
        return location;
    }

    /**
     * Compute the sum of the squares of a complex number
     * @param r real part
     * @param imag imaginary part
     * @return  r*r + imag * imag.
     */
    public static final double magnitude(double r, double imag) {
        return r * r + imag * imag;
    }

    double getMaxValue(double in[]) {
        double max;
        max = -0.99e30;

        for (int i = 0; i < in.length; i++)
            if (in[i] > max)
                max = in[i];
        return max;

    }

    void bitReverse2() {
        /* bit reversal */
        int n = realPart.length;
        int j = 1;
        int k;

        for (int i = 1; i < n; i++) {
            if (i < j) {
                swapInt(i, j);
            } //if

            k = n / 2;
            while (k >= 1 && k < j) {
                j = j - k;
                k = k / 2;
            }
            j = j + k;

        } // for
    }


    int bitr(int j, int numBits) {
        int ans = 0;
        for (int i = 0; i < numBits; i++) {
            ans = (ans << 1) + (j & 1);
            j = j >> 1;
        }
        return ans;
    }

    public void forwardFFT(double[] in_r, double[] in_i) {
        direction = FORWARD_FFT;
        fft(in_r, in_i);
    }

    public void reverseFFT(double[] in_r, double[] in_i) {
        direction = REVERSE_FFT;
        fft(in_r, in_i);
    }

    /**
     *
     * @param d input argument to the log2 function
     * @return  base 2 log.
     */
    public int log2(double d) {
        return (int) (Math.log(d) / Math.log(2));
    }

    /**
     * FFT engine.
     */
    public void fft(double _realPart[], double _imaginaryPart[]) {
        int id;
        // radix 2 number if sample, 1D of course.
        int localN;
        double wtemp, Wjk_r, Wjk_i, Wj_r, Wj_i;
        double theta, tempr, tempi;

        int numBits = log2(_realPart.length);
        // Truncate input data to a power of two
        int length = 1 << numBits; // length = 2**nu
        int n = length;

        // Copy passed references to variables to be used within
        // transforms.fft routines & utilities
        realPart = _realPart;
        imaginaryPart = _imaginaryPart;

        bitReverse2();
        for (int m = 1; m <= numBits; m++) {
            // localN = 2^m;
            localN = 1 << m;
            Wjk_r = 1;
            Wjk_i = 0;

            theta = twoPI / localN;


            Wj_r = Math.cos(theta);
            Wj_i = direction * Math.sin(theta);

            int nby2 = localN / 2;
            for (int j = 0; j < nby2; j++) {
                // This is the FFT innermost loop
                // Any optimizations that can be made here will yield
                // great rewards.
                for (int k = j; k < n; k += localN) {
                    id = k + nby2;
                    tempr = Wjk_r * realPart[id] - Wjk_i * imaginaryPart[id];
                    tempi = Wjk_r * imaginaryPart[id] + Wjk_i * realPart[id];

                    // Zid = Zi -C
                    realPart[id] = realPart[k] - tempr;
                    imaginaryPart[id] = imaginaryPart[k] - tempi;
                    realPart[k] += tempr;
                    imaginaryPart[k] += tempi;
                }

                // (eq 6.23) and (eq 6.24)

                wtemp = Wjk_r;

                Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
                Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
            }
        }
    }

    /**
     *
     * @return  gets the real component from FFT.
     */
    public double[] getRealData() {
        return realPart;
    }

    /**
     *
     * @return gets the imaginary component from FFT.
     */
    public double[] getImaginaryData() {
        return imaginaryPart;
    }

    /**
     * Compute the power spectral density of the input
     * arrays
     * @param in_r  real part of an fft
     * @param in_i  imaginary part of an fft
     * @return  the psd.
     */
    public double[] magnitudeSpectrum(double[] in_r, double[] in_i) {
        int N = in_r.length;
        double[] mag = new double[N];

        for (int i = 0; i < in_r.length; i++)
            mag[i] = Math.sqrt(in_r[i] * in_r[i] + in_i[i] * in_i[i]);

        return (mag);
    }

    public static void testFFT() {
        FFT1d f = new FFT1d();
        //f.visuallyTest();
        f.timeFFT(1024);
    }
    public static void main(String args[]){
        testFFT();
    }
}