commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From celes...@apache.org
Subject svn commit: r1241807 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/transform/ test/java/org/apache/commons/math/transform/
Date Wed, 08 Feb 2012 08:10:56 GMT
Author: celestin
Date: Wed Feb  8 08:10:56 2012
New Revision: 1241807

URL: http://svn.apache.org/viewvc?rev=1241807&view=rev
Log:
Speed improvements to o.a.c.m.transform.FastFourierTransformer. Patch contributed by Kurt
Ostfeld (MATH-732).

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/TransformUtils.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/transform/FastFourierTransformerTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java?rev=1241807&r1=1241806&r2=1241807&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java
Wed Feb  8 08:10:56 2012
@@ -18,11 +18,11 @@ package org.apache.commons.math.transfor
 
 import java.io.Serializable;
 import java.lang.reflect.Array;
+import java.util.Arrays;
 
 import org.apache.commons.math.analysis.FunctionUtils;
 import org.apache.commons.math.analysis.UnivariateFunction;
 import org.apache.commons.math.complex.Complex;
-import org.apache.commons.math.complex.RootsOfUnity;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
@@ -85,7 +85,53 @@ import org.apache.commons.math.util.Fast
 public class FastFourierTransformer implements Serializable {
 
     /** Serializable version identifier. */
-    static final long serialVersionUID = 20120501L;
+    static final long serialVersionUID = 20120802L;
+
+    /**
+     * {@code W_SUB_N_R[i]} is the real part of
+     * {@code exp(- 2 * i * pi / n)}:
+     * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}.
+     */
+    private static final double[] W_SUB_N_R =
+            {  0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1
+            , 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1
+            , 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1
+            , 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1
+            , 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1
+            , 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1
+            , 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1
+            , 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
+            , 0x1.0p0, 0x1.0p0, 0x1.0p0 };
+
+    /**
+     * {@code W_SUB_N_I[i]} is the imaginary part of
+     * {@code exp(- 2 * i * pi / n)}:
+     * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}.
+     */
+    private static final double[] W_SUB_N_I =
+            {  0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1
+            , -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5
+            , -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9
+            , -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13
+            , -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17
+            , -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21
+            , -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25
+            , -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29
+            , -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33
+            , -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37
+            , -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41
+            , -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45
+            , -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49
+            , -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53
+            , -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57
+            , -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 };
 
     /**
      * {@code true} if the unitary version of the DFT should be used.
@@ -95,9 +141,6 @@ public class FastFourierTransformer impl
      */
     private final boolean unitary;
 
-    /** The roots of unity. */
-    private RootsOfUnity roots = new RootsOfUnity();
-
     /**
      * Creates a new instance of this class, with various normalization
      * conventions.
@@ -136,6 +179,210 @@ public class FastFourierTransformer impl
         return new FastFourierTransformer(true);
     }
 
+    public static void bitReversalShuffle2(double[] a, double[] b) {
+        final int n = a.length;
+        assert(b.length == n);
+        final int halfOfN = n >> 1;
+
+        int j = 0;
+        for (int i = 0; i < n; i++) {
+            if (i < j) {
+                // swap indices i & j
+                double temp = a[i];
+                a[i] = a[j];
+                a[j] = temp;
+
+                temp = b[i];
+                b[i] = b[j];
+                b[j] = temp;
+            }
+
+            int k = halfOfN;
+            while (k <= j && k > 0) {
+                j -= k;
+                k >>= 1;
+            }
+            j += k;
+        }
+    }
+
+    /**
+     * Computes the standard transform of the specified complex data. The
+     * computation is done in place. The input data is laid out as follows
+     * <ul>
+     * <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point,
+     * </li>
+     * <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data
+     * point.</li>
+     * </ul>
+     *
+     * @param dataRI the two dimensional array of real and imaginary parts of
+     * the data
+     * @param inverse {@code true} if the inverse standard transform must be
+     * performed
+     * @throws DimensionMismatchException if the number of rows of the specified
+     * array is not two, or the array is not rectangular
+     * @throws MathIllegalArgumentException if the number of data points is not
+     * a power of two
+     */
+    public static void transformInPlace(final double[][] dataRI,
+        boolean inverse) throws
+        DimensionMismatchException, MathIllegalArgumentException {
+
+        if (dataRI.length != 2) {
+            throw new DimensionMismatchException(dataRI.length, 2);
+        }
+        final double[] dataR = dataRI[0];
+        final double[] dataI = dataRI[1];
+        if (dataR.length != dataI.length) {
+            throw new DimensionMismatchException(dataI.length, dataR.length);
+        }
+
+        final int n = dataR.length;
+        if (!ArithmeticUtils.isPowerOfTwo(n)) {
+            throw new MathIllegalArgumentException(
+                LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
+                Integer.valueOf(n));
+        }
+
+        if (n == 1) {
+            return;
+        } else if (n == 2) {
+            final double srcR0 = dataR[0];
+            final double srcI0 = dataI[0];
+            final double srcR1 = dataR[1];
+            final double srcI1 = dataI[1];
+
+            // X_0 = x_0 + x_1
+            dataR[0] = srcR0 + srcR1;
+            dataI[0] = srcI0 + srcI1;
+            // X_1 = x_0 - x_1
+            dataR[1] = srcR0 - srcR1;
+            dataI[1] = srcI0 - srcI1;
+
+            if (inverse) {
+                dataR[0] /= 2;
+                dataI[0] /= 2;
+                dataR[1] /= 2;
+                dataI[1] /= 2;
+            }
+            return;
+        }
+
+        bitReversalShuffle2(dataR, dataI);
+
+        // Do 4-term DFT.
+        if (inverse) {
+            for (int i0 = 0; i0 < n; i0 += 4) {
+                final int i1 = i0 + 1;
+                final int i2 = i0 + 2;
+                final int i3 = i0 + 3;
+
+                final double srcR0 = dataR[i0];
+                final double srcI0 = dataI[i0];
+                final double srcR1 = dataR[i2];
+                final double srcI1 = dataI[i2];
+                final double srcR2 = dataR[i1];
+                final double srcI2 = dataI[i1];
+                final double srcR3 = dataR[i3];
+                final double srcI3 = dataI[i3];
+
+                // 4-term DFT
+                // X_0 = x_0 + x_1 + x_2 + x_3
+                dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
+                dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
+                // X_1 = x_0 - x_2 + j * (x_3 - x_1)
+                dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1);
+                dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3);
+                // X_2 = x_0 - x_1 + x_2 - x_3
+                dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
+                dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
+                // X_3 = x_0 - x_2 + j * (x_1 - x_3)
+                dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3);
+                dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1);
+            }
+        } else {
+            for (int i0 = 0; i0 < n; i0 += 4) {
+                final int i1 = i0 + 1;
+                final int i2 = i0 + 2;
+                final int i3 = i0 + 3;
+
+                final double srcR0 = dataR[i0];
+                final double srcI0 = dataI[i0];
+                final double srcR1 = dataR[i2];
+                final double srcI1 = dataI[i2];
+                final double srcR2 = dataR[i1];
+                final double srcI2 = dataI[i1];
+                final double srcR3 = dataR[i3];
+                final double srcI3 = dataI[i3];
+
+                // 4-term DFT
+                // X_0 = x_0 + x_1 + x_2 + x_3
+                dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
+                dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
+                // X_1 = x_0 - x_2 + j * (x_3 - x_1)
+                dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3);
+                dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1);
+                // X_2 = x_0 - x_1 + x_2 - x_3
+                dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
+                dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
+                // X_3 = x_0 - x_2 + j * (x_1 - x_3)
+                dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1);
+                dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3);
+            }
+        }
+
+        int lastN0 = 4;
+        int lastLogN0 = 2;
+        while (lastN0 < n) {
+            int n0 = lastN0 << 1;
+            int logN0 = lastLogN0 + 1;
+            double wSubN0R = W_SUB_N_R[logN0];
+            double wSubN0I = W_SUB_N_I[logN0];
+            if (inverse) {
+                wSubN0I = -wSubN0I;
+            }
+
+            // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0
* 2).
+            for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex
+= n0) {
+                int destOddStartIndex = destEvenStartIndex + lastN0;
+
+                double wSubN0ToRR = 1;
+                double wSubN0ToRI = 0;
+
+                for (int r = 0; r < lastN0; r++) {
+                    double grR = dataR[destEvenStartIndex + r];
+                    double grI = dataI[destEvenStartIndex + r];
+                    double hrR = dataR[destOddStartIndex + r];
+                    double hrI = dataI[destOddStartIndex + r];
+
+                    // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr
+                    dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI *
hrI;
+                    dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI *
hrR;
+                    // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr
+                    dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI *
hrI);
+                    dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI *
hrR);
+
+                    // WsubN0ToR *= WsubN0R
+                    double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I;
+                    double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R;
+                    wSubN0ToRR = nextWsubN0ToRR;
+                    wSubN0ToRI = nextWsubN0ToRI;
+                }
+            }
+
+            lastN0 = n0;
+            lastLogN0 = logN0;
+        }
+
+        if (inverse) {
+            final double scaleFactor = 1.0 / ((double) n);
+            for (int i = 0; i < n; i++) {
+                dataR[i] *= scaleFactor;
+                dataI[i] *= scaleFactor;
+            }
+        }
+    }
 
     /**
      * Returns the forward transform of the specified real data set.
@@ -146,11 +393,19 @@ public class FastFourierTransformer impl
      * not a power of two
      */
     public Complex[] transform(double[] f) {
+        final double[][] dataRI = new double[][] {
+            Arrays.copyOf(f, f.length), new double[f.length]
+        };
+
+        transformInPlace(dataRI, false);
+
         if (unitary) {
             final double s = 1.0 / FastMath.sqrt(f.length);
-            return TransformUtils.scaleArray(fft(f, false), s);
+            TransformUtils.scaleArray(dataRI[0], s);
+            TransformUtils.scaleArray(dataRI[1], s);
         }
-        return fft(f, false);
+
+        return TransformUtils.createComplexArray(dataRI);
     }
 
     /**
@@ -173,11 +428,7 @@ public class FastFourierTransformer impl
             double min, double max, int n) {
 
         final double[] data = FunctionUtils.sample(f, min, max, n);
-        if (unitary) {
-            final double s = 1.0 / FastMath.sqrt(n);
-            return TransformUtils.scaleArray(fft(data, false), s);
-        }
-        return fft(data, false);
+        return transform(data);
     }
 
     /**
@@ -189,12 +440,17 @@ public class FastFourierTransformer impl
      * not a power of two
      */
     public Complex[] transform(Complex[] f) {
-        roots.computeRoots(-f.length);
+        final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);
+
+        transformInPlace(dataRI, false);
+
         if (unitary) {
             final double s = 1.0 / FastMath.sqrt(f.length);
-            return TransformUtils.scaleArray(fft(f), s);
+            TransformUtils.scaleArray(dataRI[0], s);
+            TransformUtils.scaleArray(dataRI[1], s);
         }
-        return fft(f);
+
+        return TransformUtils.createComplexArray(dataRI);
     }
 
     /**
@@ -206,8 +462,19 @@ public class FastFourierTransformer impl
      * not a power of two
      */
     public Complex[] inverseTransform(double[] f) {
-        final double s = 1.0 / (unitary ? FastMath.sqrt(f.length) : f.length);
-        return TransformUtils.scaleArray(fft(f, true), s);
+        final double[][] dataRI = new double[][] {
+            Arrays.copyOf(f, f.length), new double[f.length]
+        };
+
+        transformInPlace(dataRI, true);
+
+        if (unitary) {
+            final double s = FastMath.sqrt(f.length);
+            TransformUtils.scaleArray(dataRI[0], s);
+            TransformUtils.scaleArray(dataRI[1], s);
+        }
+
+        return TransformUtils.createComplexArray(dataRI);
     }
 
     /**
@@ -229,8 +496,7 @@ public class FastFourierTransformer impl
     public Complex[] inverseTransform(UnivariateFunction f,
             double min, double max, int n) {
         final double[] data = FunctionUtils.sample(f, min, max, n);
-        final double s = 1.0 / (unitary ? FastMath.sqrt(n) : n);
-        return TransformUtils.scaleArray(fft(data, true), s);
+        return inverseTransform(data);
     }
 
     /**
@@ -242,141 +508,19 @@ public class FastFourierTransformer impl
      * not a power of two
      */
     public Complex[] inverseTransform(Complex[] f) {
-        roots.computeRoots(f.length);
-        final double s = 1.0 / (unitary ? FastMath.sqrt(f.length) : f.length);
-        return TransformUtils.scaleArray(fft(f), s);
-    }
-
-    /**
-     * Returns the FFT of the specified real data set. Performs the base-4
-     * Cooley-Tukey FFT algorithm.
-     *
-     * @param f the real data array to be transformed
-     * @param isInverse {@code true} if inverse transform is to be carried out
-     * @return the complex transformed array
-     * @throws MathIllegalArgumentException if the length of the data array is
-     * not a power of two
-     */
-    protected Complex[] fft(double[] f, boolean isInverse) {
-
-        if (!ArithmeticUtils.isPowerOfTwo(f.length)) {
-            throw new MathIllegalArgumentException(
-                    LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
-                    Integer.valueOf(f.length));
-        }
-        Complex[] transformed = new Complex[f.length];
-        if (f.length == 1) {
-            transformed[0] = new Complex(f[0], 0.0);
-            return transformed;
-        }
-
-        // 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[] repacked = new Complex[n];
-        for (int i = 0; i < n; i++) {
-            repacked[i] = new Complex(f[2 * i], f[2 * i + 1]);
-        }
-        roots.computeRoots(isInverse ? n : -n);
-        Complex[] z = fft(repacked);
-
-        // reconstruct the FFT result for the original array
-        roots.computeRoots(isInverse ? 2 * n : -2 * n);
-        transformed[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);
-        transformed[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 = roots.getOmega(i).multiply(Complex.I);
-            Complex d = new Complex(-roots.getImaginary(i),
-                                    roots.getReal(i));
-            transformed[i] = b.subtract(c.multiply(d));
-            transformed[2 * n - i] = transformed[i].conjugate();
-        }
+        final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);
+        final double[] dataR = dataRI[0];
+        final double[] dataI = dataRI[1];
 
-        return TransformUtils.scaleArray(transformed, 0.5);
-    }
-
-    /**
-     * Returns the FFT of the specified complex data set. Performs the base-4
-     * Cooley-Tukey FFT algorithm.
-     *
-     * @param data the complex data array to be transformed
-     * @return the complex transformed array
-     * @throws MathIllegalArgumentException if the length of the data array is
-     * not a power of two
-     */
-    protected Complex[] fft(Complex[] data) {
-
-        if (!ArithmeticUtils.isPowerOfTwo(data.length)) {
-            throw new MathIllegalArgumentException(
-                    LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
-                    Integer.valueOf(data.length));
-        }
+        transformInPlace(dataRI, true);
 
-        final int n = data.length;
-        final Complex[] f = new Complex[n];
-
-        // initial simple cases
-        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;
+        if (unitary) {
+            final double s = FastMath.sqrt(f.length);
+            TransformUtils.scaleArray(dataR, s);
+            TransformUtils.scaleArray(dataI, s);
         }
 
-        // permute original data array in bit-reversal order
-        int ii = 0;
-        for (int i = 0; i < n; i++) {
-            f[i] = data[ii];
-            int k = n >> 1;
-            while (ii >= k && k > 0) {
-                ii -= k; k >>= 1;
-            }
-            ii += k;
-        }
-
-        // the bottom base-4 round
-        for (int i = 0; i < n; i += 4) {
-            final Complex a = f[i].add(f[i + 1]);
-            final Complex b = f[i + 2].add(f[i + 3]);
-            final Complex c = f[i].subtract(f[i + 1]);
-            final Complex d = f[i + 2].subtract(f[i + 3]);
-            final Complex e1 = c.add(d.multiply(Complex.I));
-            final Complex e2 = 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] = roots.isCounterClockWise() ? e1 : e2;
-            f[i + 3] = roots.isCounterClockWise() ? e2 : e1;
-        }
-
-        // iterations from bottom to top take O(N*logN) time
-        for (int i = 4; i < n; i <<= 1) {
-            final int m = n / (i << 1);
-            for (int j = 0; j < n; j += i << 1) {
-                for (int k = 0; k < i; k++) {
-                    //z = f[i+j+k].multiply(roots.getOmega(k*m));
-                    final int km = k * m;
-                    final double omegaKmReal = roots.getReal(km);
-                    final double omegaKmImag = roots.getImaginary(km);
-                    //z = f[i+j+k].multiply(omega[k*m]);
-                    final Complex z = new Complex(
-                        f[i + j + k].getReal() * omegaKmReal -
-                        f[i + j + k].getImaginary() * omegaKmImag,
-                        f[i + j + k].getReal() * omegaKmImag +
-                        f[i + j + k].getImaginary() * omegaKmReal);
-
-                    f[i + j + k] = f[j + k].subtract(z);
-                    f[j + k] = f[j + k].add(z);
-                }
-            }
-        }
-        return f;
+        return TransformUtils.createComplexArray(dataRI);
     }
 
     /**
@@ -395,10 +539,7 @@ public class FastFourierTransformer impl
      * @return transform of {@code mdca} as a Multi-Dimensional Complex Array
      * id est {@code Complex[][][][]}
      * @throws IllegalArgumentException if any dimension is not a power of two
-     * @deprecated see
-     * <a href="https://issues.apache.org/jira/browse/MATH-736">MATH-736</a>
      */
-    @Deprecated
     public Object mdfft(Object mdca, boolean forward) {
         MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)
                 new MultiDimensionalComplexMatrix(mdca).clone();
@@ -467,10 +608,7 @@ public class FastFourierTransformer impl
      * eventually be replaced by jsr-83 of the java community process
      * http://jcp.org/en/jsr/detail?id=83
      * may require additional exception throws for other basic requirements.
-     *
-     * @deprecated see MATH-736
      */
-    @Deprecated
     private static class MultiDimensionalComplexMatrix
         implements Cloneable {
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/TransformUtils.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/TransformUtils.java?rev=1241807&r1=1241806&r2=1241807&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/TransformUtils.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/TransformUtils.java
Wed Feb  8 08:10:56 2012
@@ -16,7 +16,12 @@
  */
 package org.apache.commons.math.transform;
 
+import java.util.Arrays;
+
 import org.apache.commons.math.complex.Complex;
+import org.apache.commons.math.exception.DimensionMismatchException;
+import org.apache.commons.math.exception.MathIllegalArgumentException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
 
 /**
  * Useful functions for the implementation of various transforms.
@@ -25,6 +30,20 @@ import org.apache.commons.math.complex.C
  * @since 3.0
  */
 public class TransformUtils {
+    /**
+     * Table of the powers of 2 to facilitate binary search lookup.
+     *
+     * @see #exactLog2(int)
+     */
+    private static final int[] POWERS_OF_TWO = {
+        0x00000001, 0x00000002, 0x00000004, 0x00000008, 0x00000010, 0x00000020,
+        0x00000040, 0x00000080, 0x00000100, 0x00000200, 0x00000400, 0x00000800,
+        0x00001000, 0x00002000, 0x00004000, 0x00008000, 0x00010000, 0x00020000,
+        0x00040000, 0x00080000, 0x00100000, 0x00200000, 0x00400000, 0x00800000,
+        0x01000000, 0x02000000, 0x04000000, 0x08000000, 0x10000000, 0x20000000,
+        0x40000000
+    };
+
     /** Private constructor. */
     private TransformUtils() {
         super();
@@ -62,4 +81,85 @@ public class TransformUtils {
         return f;
     }
 
+
+    /**
+     * Builds a new two dimensional array of {@code double} filled with the real
+     * and imaginary parts of the specified {@link Complex} numbers. In the
+     * returned array {@code dataRI}, the data is laid out as follows
+     * <ul>
+     * <li>{@code dataRI[0][i] = dataC[i].getReal()},</li>
+     * <li>{@code dataRI[1][i] = dataC[i].getImaginary()}.</li>
+     * </ul>
+     *
+     * @param dataC the array of {@link Complex} data to be transformed
+     * @return a two dimensional array filled with the real and imaginary parts
+     * of the specified complex input
+     */
+    public static double[][] createRealImaginaryArray(final Complex[] dataC) {
+        final double[][] dataRI = new double[2][dataC.length];
+        final double[] dataR = dataRI[0];
+        final double[] dataI = dataRI[1];
+        for (int i = 0; i < dataC.length; i++) {
+            final Complex c = dataC[i];
+            dataR[i] = c.getReal();
+            dataI[i] = c.getImaginary();
+        }
+        return dataRI;
+    }
+
+    /**
+     * Builds a new array of {@link Complex} from the specified two dimensional
+     * array of real and imaginary parts. In the returned array {@code dataC},
+     * the data is laid out as follows
+     * <ul>
+     * <li>{@code dataC[i].getReal() = dataRI[0][i]},</li>
+     * <li>{@code dataC[i].getImaginary() = dataRI[1][i]}.</li>
+     * </ul>
+     *
+     * @param dataRI the array of real and imaginary parts to be transformed
+     * @return an array of {@link Complex} with specified real and imaginary
+     * parts.
+     * @throws DimensionMismatchException if the number of rows of the specified
+     * array is not two, or the array is not rectangular
+     */
+    public static Complex[] createComplexArray(final double[][] dataRI)
+        throws DimensionMismatchException{
+
+        if (dataRI.length != 2) {
+            throw new DimensionMismatchException(dataRI.length, 2);
+        }
+        final double[] dataR = dataRI[0];
+        final double[] dataI = dataRI[1];
+        if (dataR.length != dataI.length) {
+            throw new DimensionMismatchException(dataI.length, dataR.length);
+        }
+
+        final int n = dataR.length;
+        final Complex[] c = new Complex[n];
+        for (int i = 0; i < n; i++) {
+            c[i] = new Complex(dataR[i], dataI[i]);
+        }
+        return c;
+    }
+
+
+    /**
+     * Returns the base-2 logarithm of the specified {@code int}. Throws an
+     * exception if {@code n} is not a power of two.
+     *
+     * @param n the {@code int} whose base-2 logarithm is to be evaluated
+     * @return the base-2 logarithm of {@code n}
+     * @throws MathIllegalArgumentException if {@code n} is not a power of two
+     */
+    public static int exactLog2(final int n)
+        throws MathIllegalArgumentException {
+
+        int index = Arrays.binarySearch(TransformUtils.POWERS_OF_TWO, n);
+        if (index < 0) {
+            throw new MathIllegalArgumentException(
+                    LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
+                    Integer.valueOf(n));
+        }
+        return index;
+    }
 }

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/transform/FastFourierTransformerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/transform/FastFourierTransformerTest.java?rev=1241807&r1=1241806&r2=1241807&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/transform/FastFourierTransformerTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/transform/FastFourierTransformerTest.java
Wed Feb  8 08:10:56 2012
@@ -380,7 +380,7 @@ public final class FastFourierTransforme
         doTestTransformComplex(8, 1.0E-14, forward, standard);
         doTestTransformComplex(16, 1.0E-13, forward, standard);
         doTestTransformComplex(32, 1.0E-13, forward, standard);
-        doTestTransformComplex(64, 1.0E-13, forward, standard);
+        doTestTransformComplex(64, 1.0E-12, forward, standard);
         doTestTransformComplex(128, 1.0E-12, forward, standard);
     }
 
@@ -468,7 +468,7 @@ public final class FastFourierTransforme
         doTestTransformComplex(8, 1.0E-14, forward, standard);
         doTestTransformComplex(16, 1.0E-13, forward, standard);
         doTestTransformComplex(32, 1.0E-13, forward, standard);
-        doTestTransformComplex(64, 1.0E-13, forward, standard);
+        doTestTransformComplex(64, 1.0E-12, forward, standard);
         doTestTransformComplex(128, 1.0E-12, forward, standard);
     }
 
@@ -505,7 +505,7 @@ public final class FastFourierTransforme
     public void testUnitaryInverseTransformComplex() {
         final boolean forward = false;
         final boolean standard = false;
-        doTestTransformComplex(2, 1.0E-15, forward, standard);
+        doTestTransformComplex(2, 1.0E-14, forward, standard);
         doTestTransformComplex(4, 1.0E-14, forward, standard);
         doTestTransformComplex(8, 1.0E-14, forward, standard);
         doTestTransformComplex(16, 1.0E-13, forward, standard);



Mime
View raw message