Added support for multi-dimensional Fourier transform.

JIRA: MATH-152

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@724198 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-12-07 20:48:00 +00:00
parent 71a29150bb
commit 1021f180e5
4 changed files with 226 additions and 1 deletions

View File

@ -99,6 +99,9 @@
<contributor>
<name>C. Scott Ananian</name>
</contributor>
<contributor>
<name>R&#233;mi Arntzen</name>
</contributor>
<contributor>
<name>Paul Cowan</name>
</contributor>

View File

@ -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);
}
}
}
}

View File

@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="add" issue="MATH-152" due-to="Remi Arntzen">
Added support for multi-dimensional Fourier transform.
</action>
<action dev="luc" type="update" issue="MATH-218" >
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

View File

@ -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.
*/