diff --git a/pom.xml b/pom.xml
index 50e8a9cd4..8c7b1a96a 100644
--- a/pom.xml
+++ b/pom.xml
@@ -99,6 +99,9 @@
C. Scott Ananian
+
+ Rémi Arntzen
+
Paul Cowan
diff --git a/src/java/org/apache/commons/math/transform/FastFourierTransformer.java b/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
index c2baa5dc5..59757a772 100644
--- a/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
+++ b/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
@@ -16,7 +16,9 @@
*/
package org.apache.commons.math.transform;
+import java.lang.reflect.Array;
import java.io.Serializable;
+import java.util.Arrays;
import org.apache.commons.math.analysis.*;
import org.apache.commons.math.complex.*;
import org.apache.commons.math.MathException;
@@ -562,4 +564,189 @@ public class FastFourierTransformer implements Serializable {
", " + upper + "]");
}
}
+
+ /**
+ * Performs a multi-dimensional Fourier transform on a given
+ * array, using {@link #inversetransform2(Complex[])} and
+ * {@link #transform2(Complex[])} in a row-column implementation
+ * in any number of dimensions with Θ(N×log(N)) complexity with
+ * N=n_1×n_2×n_3×⋯×n_d, n_x=number of elements in dimension x,
+ * and d=total number of dimensions.
+ *
+ * @param forward inverseTransform2 is preformed if this is false
+ * @param mdca Multi-Dimensional Complex Array id est Complex[][][][]
+ * @throws MathException if any dimension is not a power of two
+ */
+ public Object mdfft(Object mdca, boolean forward) throws MathException {
+ MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)
+ new MultiDimensionalComplexMatrix(mdca).clone();
+ int[] dimensionSize = mdcm.getDimensionSizes();
+ //cycle through each dimension
+ for (int i = 0; i < dimensionSize.length; i++) {
+ mdfft(mdcm, forward, i, new int[0]);
+ }
+ return mdcm.getArray();
+ }
+
+ private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward,
+ int d, int[] subVector) throws MathException {
+ int[] dimensionSize = mdcm.getDimensionSizes();
+ //if done
+ if (subVector.length == dimensionSize.length) {
+ Complex[] temp = new Complex[dimensionSize[d]];
+ for (int i = 0; i < dimensionSize[d]; i++) {
+ //fft along dimension d
+ subVector[d] = i;
+ temp[i] = mdcm.get(subVector);
+ }
+
+ if (forward)
+ temp = transform2(temp);
+ else
+ temp = inversetransform2(temp);
+
+ for (int i = 0; i < dimensionSize[d]; i++) {
+ subVector[d] = i;
+ mdcm.set(temp[i], subVector);
+ }
+ } else {
+ int[] vector = new int[subVector.length + 1];
+ System.arraycopy(subVector, 0, vector, 0, subVector.length);
+ if (subVector.length == d) {
+ //value is not important once the recursion is done.
+ //then an fft will be applied along the dimension d.
+ vector[d] = 0;
+ mdfft(mdcm, forward, d, vector);
+ } else {
+ for (int i = 0; i < dimensionSize[subVector.length]; i++) {
+ vector[subVector.length] = i;
+ //further split along the next dimension
+ mdfft(mdcm, forward, d, vector);
+ }
+ }
+ }
+ return;
+ }
+
+ /*
+ * not designed for synchronized access
+ * may 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.
+ */
+ private class MultiDimensionalComplexMatrix implements Serializable,
+ Cloneable {
+ private static final long serialVersionUID = 0x564FCD47EBA8169BL;
+
+ protected int[] dimensionSize = new int[1];
+ protected Object multiDimensionalComplexArray;
+
+ public MultiDimensionalComplexMatrix(Object
+ multiDimensionalComplexArray) {
+ this.multiDimensionalComplexArray = multiDimensionalComplexArray;
+ int numOfDimensions = 0;
+
+ Object lastDimension = multiDimensionalComplexArray;
+ while(lastDimension instanceof Object[]) {
+ numOfDimensions++;
+ //manually implement variable size int[]
+ if (dimensionSize.length < numOfDimensions) {
+ int[] newDimensionSize = new int[(int) Math.ceil(
+ dimensionSize.length*1.6)];
+ System.arraycopy(dimensionSize, 0, newDimensionSize, 0,
+ dimensionSize.length);
+ dimensionSize = newDimensionSize;
+ }
+ dimensionSize[numOfDimensions - 1] = ((Object[])
+ lastDimension).length;
+ lastDimension = ((Object[]) lastDimension)[0];
+ }
+ if (dimensionSize.length > numOfDimensions) {
+ int[] newDimensionSize = new int[numOfDimensions];
+ System.arraycopy(dimensionSize, 0, newDimensionSize, 0,
+ numOfDimensions);
+ dimensionSize = newDimensionSize;
+ }
+ }
+
+ public Complex get(int... vector) {
+ if ((vector == null && dimensionSize.length > 1) ||
+ vector.length != dimensionSize.length) {
+ throw new IllegalArgumentException("Number of dimensions must "
+ + "match");
+ }
+
+ Object lastDimension = multiDimensionalComplexArray;
+
+ for (int i = 0; i < dimensionSize.length; i++) {
+ lastDimension = ((Object[]) lastDimension)[vector[i]];
+ }
+ return (Complex) lastDimension;
+ }
+
+ public Complex set(Complex magnitude, int... vector) {
+ if ((vector == null && dimensionSize.length > 1) ||
+ vector.length != dimensionSize.length) {
+ throw new IllegalArgumentException("Number of dimensions must "
+ + "match");
+ }
+
+ Object lastDimension = multiDimensionalComplexArray;
+
+ for (int i = 0; i < dimensionSize.length - 1; i++) {
+ lastDimension = ((Object[]) lastDimension)[vector[i]];
+ }
+
+ Complex lastValue = (Complex) ((Object[])
+ lastDimension)[vector[dimensionSize.length - 1]];
+ ((Object[]) lastDimension)[vector[dimensionSize.length - 1]] =
+ magnitude;
+ return lastValue;
+ }
+
+ public int[] getDimensionSizes() {
+ return dimensionSize.clone();
+ }
+
+ public Object getArray() {
+ return multiDimensionalComplexArray;
+ }
+
+ @Override
+ public Object clone() {
+ MultiDimensionalComplexMatrix mdcm =
+ new MultiDimensionalComplexMatrix(Array.newInstance(
+ Complex.class, dimensionSize));
+ clone(mdcm);
+ return mdcm;
+ }
+
+ /*
+ * Copy contents of current array into mdcm.
+ */
+ private void clone(MultiDimensionalComplexMatrix mdcm) {
+ int[] vector = new int[dimensionSize.length];
+ int size = 1;
+ for (int i = 0; i < dimensionSize.length; i++) {
+ size *= dimensionSize[i];
+ }
+ int[][] vectorList = new int[size][dimensionSize.length];
+ for (int[] nextVector: vectorList) {
+ System.arraycopy(vector, 0, nextVector, 0,
+ dimensionSize.length);
+ for (int i = 0; i < dimensionSize.length; i++) {
+ vector[i]++;
+ if (vector[i] < dimensionSize[i]) {
+ break;
+ } else {
+ vector[i] = 0;
+ }
+ }
+ }
+
+ for (int[] nextVector: vectorList) {
+ mdcm.set(get(nextVector), nextVector);
+ }
+ }
+ }
}
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index c4c88e11c..0467ca6e7 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -39,6 +39,9 @@ The type attribute can be add,update,fix,remove.
+
+ Added support for multi-dimensional Fourier transform.
+
The root solvers now take the function to solve as a parameter to
the solve methods, thus allowing to reuse the same solver for different
diff --git a/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java b/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java
index 94636118e..49561dfce 100644
--- a/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java
+++ b/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java
@@ -77,7 +77,39 @@ public final class FastFourierTransformerTest extends TestCase {
assertEquals(y2[i].getImaginary(), result[i].getImaginary(), tolerance);
}
}
-
+
+ public void test2DData() throws MathException {
+ FastFourierTransformer transformer = new FastFourierTransformer();
+ double tolerance = 1E-12;
+ Complex[][] input = new Complex[][] {new Complex[] {new Complex(1, 0),
+ new Complex(2, 0)},
+ new Complex[] {new Complex(3, 1),
+ new Complex(4, 2)}};
+ Complex[][] goodOutput = new Complex[][] {new Complex[] {new Complex(5,
+ 1.5), new Complex(-1, -.5)}, new Complex[] {new Complex(-2,
+ -1.5), new Complex(0, .5)}};
+ Complex[][] output = (Complex[][])transformer.mdfft(input, true);
+ Complex[][] output2 = (Complex[][])transformer.mdfft(output, false);
+
+ assertEquals(input.length, output.length);
+ assertEquals(input.length, output2.length);
+ assertEquals(input[0].length, output[0].length);
+ assertEquals(input[0].length, output2[0].length);
+ assertEquals(input[1].length, output[1].length);
+ assertEquals(input[1].length, output2[1].length);
+
+ for (int i = 0; i < input.length; i++) {
+ for (int j = 0; j < input[0].length; j++) {
+ assertEquals(input[i][j].getImaginary(), output2[i][j].getImaginary(),
+ tolerance);
+ assertEquals(input[i][j].getReal(), output2[i][j].getReal(), tolerance);
+ assertEquals(goodOutput[i][j].getImaginary(), output[i][j].getImaginary(),
+ tolerance);
+ assertEquals(goodOutput[i][j].getReal(), output[i][j].getReal(), tolerance);
+ }
+ }
+ }
+
/**
* Test of transformer for the sine function.
*/