fastfouriertransformer.java

Apache的common math数学软件包

     * <p>
     * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
     * </p>
     * 
     * @param f the complex data array to be inversely transformed
     * @return the complex inversely transformed array
     * @throws MathException if any math-related errors occur
     * @throws IllegalArgumentException if any parameters are invalid
     */
    public Complex[] inversetransform2(Complex f[]) throws MathException,
        IllegalArgumentException {

        computeOmega(-f.length);    // pass negative argument
        double scaling_coefficient = 1.0 / Math.sqrt(f.length);
        return scaleArray(fft(f), scaling_coefficient);
    }

    /**
     * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
     *
     * @param f the real data array to be transformed
     * @param isInverse the indicator of forward or inverse transform
     * @return the complex transformed array
     * @throws MathException if any math-related errors occur
     * @throws IllegalArgumentException if any parameters are invalid
     */
    protected Complex[] fft(double f[], boolean isInverse) throws
        MathException, IllegalArgumentException {

        verifyDataSet(f);
        Complex F[] = new Complex[f.length];
        if (f.length == 1) {
            F[0] = new Complex(f[0], 0.0);
            return F;
        }

        // Rather than the naive real to complex conversion, pack 2N
        // real numbers into N complex numbers for better performance.
        int N = f.length >> 1;
        Complex c[] = new Complex[N];
        for (int i = 0; i < N; i++) {
            c[i] = new Complex(f[2*i], f[2*i+1]);
        }
        computeOmega(isInverse ? -N : N);
        Complex z[] = fft(c);

        // reconstruct the FFT result for the original array
        computeOmega(isInverse ? -2*N : 2*N);
        F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);
        F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0);
        for (int i = 1; i < N; i++) {
            Complex A = z[N-i].conjugate();
            Complex B = z[i].add(A);
            Complex C = z[i].subtract(A);
            Complex D = omega[i].multiply(Complex.I);
            F[i] = B.subtract(C.multiply(D));
            F[2*N-i] = F[i].conjugate();
        }

        return scaleArray(F, 0.5);
    }

    /**
     * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
     *
     * @param data the complex data array to be transformed
     * @return the complex transformed array
     * @throws MathException if any math-related errors occur
     * @throws IllegalArgumentException if any parameters are invalid
     */
    protected Complex[] fft(Complex data[]) throws MathException,
        IllegalArgumentException {

        int i, j, k, m, N = data.length;
        Complex A, B, C, D, E, F, z, f[] = new Complex[N];

        // initial simple cases
        verifyDataSet(data);
        if (N == 1) {
            f[0] = data[0];
            return f;
        }
        if (N == 2) {
            f[0] = data[0].add(data[1]);
            f[1] = data[0].subtract(data[1]);
            return f;
        }

        // permute original data array in bit-reversal order
        j = 0;
        for (i = 0; i < N; i++) {
            f[i] = data[j];
            k = N >> 1;
            while (j >= k && k > 0) {
                j -= k; k >>= 1;
            }
            j += k;
        }

        // the bottom base-4 round
        for (i = 0; i < N; i += 4) {
            A = f[i].add(f[i+1]);
            B = f[i+2].add(f[i+3]);
            C = f[i].subtract(f[i+1]);
            D = f[i+2].subtract(f[i+3]);
            E = C.add(D.multiply(Complex.I));
            F = C.subtract(D.multiply(Complex.I));
            f[i] = A.add(B);
            f[i+2] = A.subtract(B);
            // omegaCount indicates forward or inverse transform
            f[i+1] = omegaCount < 0 ? E : F;
            f[i+3] = omegaCount > 0 ? E : F;
        }

        // iterations from bottom to top take O(N*logN) time
        for (i = 4; i < N; i <<= 1) {
            m = N / (i<<1);
            for (j = 0; j < N; j += i<<1) {
                for (k = 0; k < i; k++) {
                    z = f[i+j+k].multiply(omega[k*m]);
                    f[i+j+k] = f[j+k].subtract(z);
                    f[j+k] = f[j+k].add(z);
                }
            }
        }
        return f;
    }

    /**
     * Calculate the n-th roots of unity.
     * <p>
     * The computed omega[] = { 1, w, w^2, ... w^(n-1) } where
     * w = exp(-2 \pi i / n), i = sqrt(-1). Note n is positive for
     * forward transform and negative for inverse transform. </p>
     * 
     * @param n the integer passed in
     * @throws IllegalArgumentException if n = 0
     */
    protected void computeOmega(int n) throws IllegalArgumentException {
        if (n == 0) {
            throw new IllegalArgumentException
                ("Cannot compute 0-th root of unity, indefinite result.");
        }
        // avoid repetitive calculations
        if (n == omegaCount) { return; }
        if (n + omegaCount == 0) {
            for (int i = 0; i < Math.abs(omegaCount); i++) {
                omega[i] = omega[i].conjugate();
            }
            omegaCount = n;
            return;
        }
        // calculate everything from scratch
        omega = new Complex[Math.abs(n)];
        double t = 2.0 * Math.PI / n;
        double cost = Math.cos(t);
        double sint = Math.sin(t);
        omega[0] = new Complex(1.0, 0.0);
        for (int i = 1; i < Math.abs(n); i++) {
            omega[i] = new Complex(
                omega[i-1].getReal() * cost + omega[i-1].getImaginary() * sint,
                omega[i-1].getImaginary() * cost - omega[i-1].getReal() * sint);
        }
        omegaCount = n;
    }

    /**
     * Sample the given univariate real function on the given interval.
     * <p>
     * The interval is divided equally into N sections and sample points
     * are taken from min to max-(max-min)/N. Usually f(x) is periodic
     * such that f(min) = f(max) (note max is not sampled), but we don't
     * require that.</p>
     *
     * @param f the function to be sampled
     * @param min the lower bound for the interval
     * @param max the upper bound for the interval
     * @param n the number of sample points
     * @return the samples array
     * @throws MathException if any math-related errors occur
     * @throws IllegalArgumentException if any parameters are invalid
     */
    public static double[] sample(
        UnivariateRealFunction f, double min, double max, int n)
        throws MathException, IllegalArgumentException {

        if (n <= 0) {
            throw new IllegalArgumentException("Number of samples not positive.");
        }
        verifyInterval(min, max);

        double s[] = new double[n];
        double h = (max - min) / n;
        for (int i = 0; i < n; i++) {
            s[i] = f.value(min + i * h);
        }
        return s;
    }

    /**
     * Multiply every component in the given real array by the
     * given real number. The change is made in place.
     *
     * @param f the real array to be scaled
     * @param d the real scaling coefficient
     * @return a reference to the scaled array
     */
    public static double[] scaleArray(double f[], double d) {
        for (int i = 0; i < f.length; i++) {
            f[i] *= d;
        }
        return f;
    }

    /**
     * Multiply every component in the given complex array by the
     * given real number. The change is made in place.
     *
     * @param f the complex array to be scaled
     * @param d the real scaling coefficient
     * @return a reference to the scaled array
     */
    public static Complex[] scaleArray(Complex f[], double d) {
        for (int i = 0; i < f.length; i++) {
            f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary());
        }
        return f;
    }

    /**
     * Returns true if the argument is power of 2.
     * 
     * @param n the number to test
     * @return true if the argument is power of 2
     */
    public static boolean isPowerOf2(long n) {
        return (n > 0) && ((n & (n - 1)) == 0);
    }

    /**
     * Verifies that the data set has length of power of 2.
     * 
     * @param d the data array
     * @throws IllegalArgumentException if array length is not power of 2
     */
    public static void verifyDataSet(double d[]) throws IllegalArgumentException {
        if (!isPowerOf2(d.length)) {
            throw new IllegalArgumentException
                ("Number of samples not power of 2, consider padding for fix.");
        }       
    }

    /**
     * Verifies that the data set has length of power of 2.
     * 
     * @param o the data array
     * @throws IllegalArgumentException if array length is not power of 2
     */
    public static void verifyDataSet(Object o[]) throws IllegalArgumentException {
        if (!isPowerOf2(o.length)) {
            throw new IllegalArgumentException
                ("Number of samples not power of 2, consider padding for fix.");
        }       
    }

    /**
     * Verifies that the endpoints specify an interval.
     * 
     * @param lower lower endpoint
     * @param upper upper endpoint
     * @throws IllegalArgumentException if not interval
     */
    public static void verifyInterval(double lower, double upper) throws
        IllegalArgumentException {

        if (lower >= upper) {
            throw new IllegalArgumentException
                ("Endpoints do not specify an interval: [" + lower +
                ", " + upper + "]");
        }       
    }
}