You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/02/03 20:59:20 UTC
svn commit: r740400 - in /commons/proper/math/trunk: pom.xml
src/java/org/apache/commons/math/MessagesResources_fr.java
src/java/org/apache/commons/math/transform/FastFourierTransformer.java
src/site/xdoc/changes.xml
Author: luc
Date: Tue Feb 3 19:59:20 2009
New Revision: 740400
URL: http://svn.apache.org/viewvc?rev=740400&view=rev
Log:
applied Cyril Briquet's patch (with slight changes) to improve FastFourierTransform efficiency
JIRA: MATH-216
Modified:
commons/proper/math/trunk/pom.xml
commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java
commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
Modified: commons/proper/math/trunk/pom.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=740400&r1=740399&r2=740400&view=diff
==============================================================================
--- commons/proper/math/trunk/pom.xml (original)
+++ commons/proper/math/trunk/pom.xml Tue Feb 3 19:59:20 2009
@@ -106,6 +106,9 @@
<name>Rémi Arntzen</name>
</contributor>
<contributor>
+ <name>Cyril Briquet</name>
+ </contributor>
+ <contributor>
<name>Paul Cowan</name>
</contributor>
<contributor>
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java?rev=740400&r1=740399&r2=740400&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java Tue Feb 3 19:59:20 2009
@@ -372,6 +372,10 @@
// org.apache.commons.math.transform.FastFourierTransformer
{ "cannot compute 0-th root of unity, indefinite result",
"impossible de calculer la racine z\u00e9roi\u00e8me de l''unit\u00e9, r\u00e9sultat ind\u00e9fini" },
+ { "roots of unity have not been computed yet",
+ "les racines de l''unit\u00e9 n''ont pas encore \u00e9t\u00e9 calcul\u00e9es" },
+ { "out of range root of unity index {0} (must be in [{1};{2}])",
+ "index de racine de l''unit\u00e9 hors domaine (devrait \u00eatre dans [{1}; {2}])" },
{ "number of sample is not positive: {0}",
"le nombre d''\u00e9chantillons n''est pas positif : {0}" },
{ "{0} is not a power of 2, consider padding for fix",
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java?rev=740400&r1=740399&r2=740400&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java Tue Feb 3 19:59:20 2009
@@ -48,14 +48,8 @@
/** Serializable version identifier. */
static final long serialVersionUID = 5138259215438106000L;
- /** array of the roots of unity */
- private Complex omega[] = new Complex[0];
-
- /**
- * |omegaCount| is the length of lasted computed omega[]. omegaCount
- * is positive for forward transform and negative for inverse transform.
- */
- private int omegaCount = 0;
+ /** roots of unity */
+ private RootsOfUnity roots = new RootsOfUnity();
/**
* Construct a default transformer.
@@ -113,7 +107,7 @@
*/
public Complex[] transform(Complex f[])
throws IllegalArgumentException {
- computeOmega(f.length);
+ roots.computeOmega(f.length);
return fft(f);
}
@@ -171,7 +165,7 @@
public Complex[] transform2(Complex f[])
throws IllegalArgumentException {
- computeOmega(f.length);
+ roots.computeOmega(f.length);
double scaling_coefficient = 1.0 / Math.sqrt(f.length);
return scaleArray(fft(f), scaling_coefficient);
}
@@ -230,7 +224,7 @@
public Complex[] inversetransform(Complex f[])
throws IllegalArgumentException {
- computeOmega(-f.length); // pass negative argument
+ roots.computeOmega(-f.length); // pass negative argument
double scaling_coefficient = 1.0 / f.length;
return scaleArray(fft(f), scaling_coefficient);
}
@@ -289,7 +283,7 @@
public Complex[] inversetransform2(Complex f[])
throws IllegalArgumentException {
- computeOmega(-f.length); // pass negative argument
+ roots.computeOmega(-f.length); // pass negative argument
double scaling_coefficient = 1.0 / Math.sqrt(f.length);
return scaleArray(fft(f), scaling_coefficient);
}
@@ -319,18 +313,20 @@
for (int i = 0; i < N; i++) {
c[i] = new Complex(f[2*i], f[2*i+1]);
}
- computeOmega(isInverse ? -N : N);
+ roots.computeOmega(isInverse ? -N : N);
Complex z[] = fft(c);
// reconstruct the FFT result for the original array
- computeOmega(isInverse ? -2*N : 2*N);
+ roots.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);
+ //Complex D = roots.getOmega(i).multiply(Complex.I);
+ Complex D = new Complex(-roots.getOmegaImaginary(i),
+ roots.getOmegaReal(i));
F[i] = B.subtract(C.multiply(D));
F[2*N-i] = F[i].conjugate();
}
@@ -385,8 +381,8 @@
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;
+ f[i+1] = roots.isForward() ? F : E;
+ f[i+3] = roots.isForward() ? E : F;
}
// iterations from bottom to top take O(N*logN) time
@@ -394,7 +390,17 @@
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]);
+ //z = f[i+j+k].multiply(roots.getOmega(k*m));
+ final int k_times_m = k*m;
+ final double omega_k_times_m_real = roots.getOmegaReal(k_times_m);
+ final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m);
+ //z = f[i+j+k].multiply(omega[k*m]);
+ z = new Complex(
+ f[i+j+k].getReal() * omega_k_times_m_real -
+ f[i+j+k].getImaginary() * omega_k_times_m_imaginary,
+ f[i+j+k].getReal() * omega_k_times_m_imaginary +
+ f[i+j+k].getImaginary() * omega_k_times_m_real);
+
f[i+j+k] = f[j+k].subtract(z);
f[j+k] = f[j+k].add(z);
}
@@ -404,45 +410,6 @@
}
/**
- * 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 MathRuntimeException.createIllegalArgumentException("cannot compute 0-th root of unity, indefinite result",
- null);
- }
- // 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
@@ -803,4 +770,143 @@
}
}
}
+
+
+ /** Computes the n<sup>th</sup> roots of unity.
+ * A cache of already computed values is maintained.
+ */
+ private static class RootsOfUnity {
+
+ private int omegaCount;
+ private double[] omegaReal;
+ private double[] omegaImaginaryForward;
+ private double[] omegaImaginaryInverse;
+ private boolean isForward;
+
+ public RootsOfUnity() {
+
+ omegaCount = 0;
+ omegaReal = null;
+ omegaImaginaryForward = null;
+ omegaImaginaryInverse = null;
+ isForward = true;
+
+ }
+
+ /**
+ * Check if computation has been done for forward or reverse transform.
+ * @return true if computation has been done for forward transform
+ * @throws IllegalStateException if no roots of unity have been computed yet
+ */
+ public synchronized boolean isForward() throws IllegalStateException {
+
+ if (omegaCount == 0) {
+ throw MathRuntimeException.createIllegalStateException(
+ "roots of unity have not been computed yet",
+ null);
+ }
+ return isForward;
+
+ }
+
+ /** Computes the n<sup>th</sup> roots of unity.
+ * <p>The computed omega[] = { 1, w, w<sup>2</sup>, ... w<sup>(n-1)</sup> } where
+ * w = exp(-2 π i / n), i = &sqrt;(-1).</p>
+ * <p>Note that n is positive for
+ * forward transform and negative for inverse transform.</p>
+ * @param n number of roots of unity to compute,
+ * positive for forward transform, negative for inverse transform
+ * @throws IllegalArgumentException if n = 0
+ */
+ public synchronized void computeOmega(int n) throws IllegalArgumentException {
+
+ if (n == 0) {
+ throw MathRuntimeException.createIllegalArgumentException(
+ "cannot compute 0-th root of unity, indefinite result",
+ null);
+ }
+
+ isForward = (n > 0);
+
+ // avoid repetitive calculations
+ final int absN = Math.abs(n);
+
+ if (absN == omegaCount) {
+ return;
+ }
+
+ // calculate everything from scratch, for both forward and inverse versions
+ final double t = 2.0 * Math.PI / absN;
+ final double cosT = Math.cos(t);
+ final double sinT = Math.sin(t);
+ omegaReal = new double[absN];
+ omegaImaginaryForward = new double[absN];
+ omegaImaginaryInverse = new double[absN];
+ omegaReal[0] = 1.0;
+ omegaImaginaryForward[0] = 0.0;
+ omegaImaginaryInverse[0] = 0.0;
+ for (int i = 1; i < absN; i++) {
+ omegaReal[i] =
+ omegaReal[i-1] * cosT + omegaImaginaryForward[i-1] * sinT;
+ omegaImaginaryForward[i] =
+ omegaImaginaryForward[i-1] * cosT - omegaReal[i-1] * sinT;
+ omegaImaginaryInverse[i] = -omegaImaginaryForward[i];
+ }
+ omegaCount = absN;
+
+ }
+
+ /**
+ * Get the real part of the k<sup>th</sup> n<sup>th</sup> root of unity
+ * @param k index of the n<sup>th</sup> root of unity
+ * @return real part of the k<sup>th</sup> n<sup>th</sup> root of unity
+ * @throws IllegalStateException if no roots of unity have been computed yet
+ * @throws IllegalArgumentException if k is out of range
+ */
+ public synchronized double getOmegaReal(int k)
+ throws IllegalStateException, IllegalArgumentException {
+
+ if (omegaCount == 0) {
+ throw MathRuntimeException.createIllegalStateException(
+ "roots of unity have not been computed yet",
+ null);
+ }
+ if ((k < 0) || (k >= omegaCount)) {
+ throw MathRuntimeException.createIllegalArgumentException(
+ "out of range root of unity index {0} (must be in [{1};{2}])",
+ new Object[] { k, 0, omegaCount - 1 });
+ }
+
+ return omegaReal[k];
+
+ }
+
+ /**
+ * Get the imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity
+ * @param k index of the n<sup>th</sup> root of unity
+ * @return imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity
+ * @throws IllegalStateException if no roots of unity have been computed yet
+ * @throws IllegalArgumentException if k is out of range
+ */
+ public synchronized double getOmegaImaginary(int k)
+ throws IllegalStateException, IllegalArgumentException {
+
+ if (omegaCount == 0) {
+ throw MathRuntimeException.createIllegalStateException(
+ "roots of unity have not been computed yet",
+ null);
+ }
+ if ((k < 0) || (k >= omegaCount)) {
+ throw MathRuntimeException.createIllegalArgumentException(
+ "out of range root of unity index {0} (must be in [{1};{2}])",
+ new Object[] { k, 0, omegaCount - 1 });
+ }
+
+ return (isForward) ?
+ omegaImaginaryForward[k] : omegaImaginaryInverse[k];
+
+ }
+
+ }
+
}
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=740400&r1=740399&r2=740400&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Tue Feb 3 19:59:20 2009
@@ -39,6 +39,9 @@
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
+ <action dev="luc" type="fix" issue="MATH-216" due-to="Cyril Briquet">
+ Improved fast Fourier transform efficiency.
+ </action>
<action dev="luc" type="add" >
Added factory methods to create Chebyshev, Hermite, Laguerre and Legendre polynomials.
</action>