commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From pste...@apache.org
Subject svn commit: r419098 - in /jakarta/commons/proper/math/trunk: src/java/org/apache/commons/math/transform/ src/test/org/apache/commons/math/transform/ xdocs/
Date Tue, 04 Jul 2006 20:58:54 GMT
Author: psteitz
Date: Tue Jul  4 13:58:54 2006
New Revision: 419098

URL: http://svn.apache.org/viewvc?rev=419098&view=rev
Log:
Added Fast Fourier transform
JIRA: MATH-140
Contributed by Xiaogang Zhang

Added:
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastCosineTransformer.java   (with props)
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java   (with props)
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastSineTransformer.java   (with props)
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html   (with props)
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastCosineTransformerTest.java   (with props)
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java   (with props)
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastSineTransformerTest.java   (with props)
Modified:
    jakarta/commons/proper/math/trunk/xdocs/changes.xml

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastCosineTransformer.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastCosineTransformer.java?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastCosineTransformer.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastCosineTransformer.java Tue Jul  4 13:58:54 2006
@@ -0,0 +1,260 @@
+/*
+ * Copyright 2005 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.transform;
+
+import java.io.Serializable;
+import org.apache.commons.math.analysis.*;
+import org.apache.commons.math.complex.*;
+import org.apache.commons.math.MathException;
+
+/**
+ * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
+ * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine Transform</a>
+ * for transformation of one-dimensional data sets. For reference, see
+ * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
+ * <p>
+ * FCT is its own inverse, up to a multiplier depending on conventions.
+ * The equations are listed in the comments of the corresponding methods.
+ * <p>
+ * Different from FFT and FST, FCT requires the length of data set to be
+ * power of 2 plus one. Users should especially pay attention to the
+ * function transformation on how this affects the sampling.
+ *
+ * @version $Revision$ $Date$
+ */
+public class FastCosineTransformer implements Serializable {
+
+    /** serializable version identifier */
+    static final long serialVersionUID = -7673941545134707766L;
+
+    /**
+     * Construct a default transformer.
+     */
+    FastCosineTransformer() {
+        super();
+    }
+
+    /**
+     * Transform the given real data set.
+     * <p>
+     * The formula is $ F_n = (1/2) [f_0 + (-1)^n f_N] +
+     *                        \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
+     *
+     * @param f the real data array to be transformed
+     * @return the real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        return fct(f);
+    }
+
+    /**
+     * Transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ F_n = (1/2) [f_0 + (-1)^n f_N] +
+     *                        \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
+     *
+     * @param f the function to be sampled and transformed
+     * @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 real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = FastFourierTransformer.sample(f, min, max, n);
+        return fct(data);
+    }
+
+    /**
+     * Transform the given real data set.
+     * <p>
+     * The formula is $ F_n = \sqrt{1/2N} [f_0 + (-1)^n f_N] +
+     *                        \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
+     *
+     * @param f the real data array to be transformed
+     * @return the real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform2(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = Math.sqrt(2.0 / (f.length-1));
+        return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
+    }
+
+    /**
+     * Transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ F_n = \sqrt{1/2N} [f_0 + (-1)^n f_N] +
+     *                        \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
+     *
+     * @param f the function to be sampled and transformed
+     * @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 real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform2(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = FastFourierTransformer.sample(f, min, max, n);
+        double scaling_coefficient = Math.sqrt(2.0 / (n-1));
+        return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real data set.
+     * <p>
+     * The formula is $ f_k = (1/N) [F_0 + (-1)^k F_N] +
+     *                        (2/N) \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
+     *
+     * @param f the real data array to be inversely transformed
+     * @return the real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = 2.0 / (f.length - 1);
+        return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ f_k = (1/N) [F_0 + (-1)^k F_N] +
+     *                        (2/N) \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
+     *
+     * @param f the function to be sampled and inversely transformed
+     * @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 real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = FastFourierTransformer.sample(f, min, max, n);
+        double scaling_coefficient = 2.0 / (n - 1);
+        return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real data set.
+     * <p>
+     * The formula is $ f_k = \sqrt{1/2N} [F_0 + (-1)^k F_N] +
+     *                        \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
+     *
+     * @param f the real data array to be inversely transformed
+     * @return the real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform2(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        return transform2(f);
+    }
+
+    /**
+     * Inversely transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ f_k = \sqrt{1/2N} [F_0 + (-1)^k F_N] +
+     *                        \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
+     *
+     * @param f the function to be sampled and inversely transformed
+     * @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 real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform2(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        return transform2(f, min, max, n);
+    }
+
+    /**
+     * Perform the FCT algorithm (including inverse).
+     *
+     * @param f the real data array to be transformed
+     * @return the real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    protected double[] fct(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double A, B, C, F1, x[], F[] = new double[f.length];
+
+        int N = f.length - 1;
+        if (!FastFourierTransformer.isPowerOf2(N)) {
+            throw new IllegalArgumentException
+                ("Number of samples not power of 2 plus one: " + f.length);
+        }
+        if (N == 1) {       // trivial case
+            F[0] = 0.5 * (f[0] + f[1]);
+            F[1] = 0.5 * (f[0] - f[1]);
+            return F;
+        }
+
+        // construct a new array and perform FFT on it
+        x = new double[N];
+        x[0] = 0.5 * (f[0] + f[N]);
+        x[N >> 1] = f[N >> 1];
+        F1 = 0.5 * (f[0] - f[N]);   // temporary variable for F[1]
+        for (int i = 1; i < (N >> 1); i++) {
+            A = 0.5 * (f[i] + f[N-i]);
+            B = Math.sin(i * Math.PI / N) * (f[i] - f[N-i]);
+            C = Math.cos(i * Math.PI / N) * (f[i] - f[N-i]);
+            x[i] = A - B;
+            x[N-i] = A + B;
+            F1 += C;
+        }
+        FastFourierTransformer transformer = new FastFourierTransformer();
+        Complex y[] = transformer.transform(x);
+
+        // reconstruct the FCT result for the original array
+        F[0] = y[0].getReal();
+        F[1] = F1;
+        for (int i = 1; i < (N >> 1); i++) {
+            F[2*i] = y[i].getReal();
+            F[2*i+1] = F[2*i-1] - y[i].getImaginary();
+        }
+        F[N] = y[N >> 1].getReal();
+
+        return F;
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastCosineTransformer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java Tue Jul  4 13:58:54 2006
@@ -0,0 +1,551 @@
+/*
+ * Copyright 2005 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.transform;
+
+import java.io.Serializable;
+import org.apache.commons.math.analysis.*;
+import org.apache.commons.math.complex.*;
+import org.apache.commons.math.MathException;
+
+/**
+ * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html">
+ * Fast Fourier Transform</a> for transformation of one-dimensional data sets.
+ * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897,
+ * chapter 6.
+ * <p>
+ * There are several conventions for the definition of FFT and inverse FFT,
+ * mainly on different coefficient and exponent. Here the equations are listed
+ * in the comments of the corresponding methods.
+ * <p>
+ * We require the length of data set to be power of 2, this greatly simplifies
+ * and speeds up the code. Users can pad the data with zeros to meet this
+ * requirement. There are other flavors of FFT, for reference, see S. Winograd,
+ * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation,
+ * 32 (1978), 175 - 199.
+ *
+ * @version $Revision$ $Date$
+ */
+public class FastFourierTransformer implements Serializable {
+
+    /** 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;
+
+    /**
+     * Construct a default transformer.
+     */
+    FastFourierTransformer() {
+        super();
+    }
+
+    /**
+     * Transform the given real data set.
+     * <p>
+     * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
+     *
+     * @param f the real 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
+     */
+    public Complex[] transform(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        return fft(f, false);
+    }
+
+    /**
+     * Transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
+     *
+     * @param f the function to be sampled and transformed
+     * @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 complex transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public Complex[] transform(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = sample(f, min, max, n);
+        return fft(data, false);
+    }
+
+    /**
+     * Transform the given complex data set.
+     * <p>
+     * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
+     *
+     * @param f 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
+     */
+    public Complex[] transform(Complex f[]) throws MathException,
+        IllegalArgumentException {
+
+        computeOmega(f.length);
+        return fft(f);
+    }
+
+    /**
+     * Transform the given real data set.
+     * <p>
+     * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
+     *
+     * @param f the real 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
+     */
+    public Complex[] transform2(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = 1.0 / Math.sqrt(f.length);
+        return scaleArray(fft(f, false), scaling_coefficient);
+    }
+
+    /**
+     * Transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
+     *
+     * @param f the function to be sampled and transformed
+     * @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 complex transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public Complex[] transform2(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = sample(f, min, max, n);
+        double scaling_coefficient = 1.0 / Math.sqrt(n);
+        return scaleArray(fft(data, false), scaling_coefficient);
+    }
+
+    /**
+     * Transform the given complex data set.
+     * <p>
+     * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
+     *
+     * @param f 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
+     */
+    public Complex[] transform2(Complex f[]) throws MathException,
+        IllegalArgumentException {
+
+        computeOmega(f.length);
+        double scaling_coefficient = 1.0 / Math.sqrt(f.length);
+        return scaleArray(fft(f), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real data set.
+     * <p>
+     * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
+     *
+     * @param f the real 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[] inversetransform(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = 1.0 / f.length;
+        return scaleArray(fft(f, true), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
+     *
+     * @param f the function to be sampled and inversely transformed
+     * @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 complex inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public Complex[] inversetransform(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = sample(f, min, max, n);
+        double scaling_coefficient = 1.0 / n;
+        return scaleArray(fft(data, true), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given complex data set.
+     * <p>
+     * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
+     *
+     * @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[] inversetransform(Complex f[]) throws MathException,
+        IllegalArgumentException {
+
+        computeOmega(-f.length);    // pass negative argument
+        double scaling_coefficient = 1.0 / f.length;
+        return scaleArray(fft(f), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real data set.
+     * <p>
+     * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
+     *
+     * @param f the real 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(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = 1.0 / Math.sqrt(f.length);
+        return scaleArray(fft(f, true), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
+     *
+     * @param f the function to be sampled and inversely transformed
+     * @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 complex inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public Complex[] inversetransform2(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = sample(f, min, max, n);
+        double scaling_coefficient = 1.0 / Math.sqrt(n);
+        return scaleArray(fft(data, true), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given complex data set.
+     * <p>
+     * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
+     *
+     * @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.
+     * 
+     * @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.
+     *
+     * @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 + "]");
+        }       
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastFourierTransformer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastSineTransformer.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastSineTransformer.java?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastSineTransformer.java (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastSineTransformer.java Tue Jul  4 13:58:54 2006
@@ -0,0 +1,251 @@
+/*
+ * Copyright 2005 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.transform;
+
+import java.io.Serializable;
+import org.apache.commons.math.analysis.*;
+import org.apache.commons.math.complex.*;
+import org.apache.commons.math.MathException;
+
+/**
+ * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
+ * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a>
+ * for transformation of one-dimensional data sets. For reference, see
+ * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
+ * <p>
+ * FST is its own inverse, up to a multiplier depending on conventions.
+ * The equations are listed in the comments of the corresponding methods.
+ * <p>
+ * Similar to FFT, we also require the length of data set to be power of 2.
+ * In addition, the first element must be 0 and it's enforced in function
+ * transformation after sampling.
+ *
+ * @version $Revision$ $Date$
+ */
+public class FastSineTransformer implements Serializable {
+
+    /** serializable version identifier */
+    static final long serialVersionUID = -478002039949390854L;
+
+    /**
+     * Construct a default transformer.
+     */
+    FastSineTransformer() {
+        super();
+    }
+
+    /**
+     * Transform the given real data set.
+     * <p>
+     * The formula is $ F_n = \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
+     *
+     * @param f the real data array to be transformed
+     * @return the real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        return fst(f);
+    }
+
+    /**
+     * Transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ F_n = \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
+     *
+     * @param f the function to be sampled and transformed
+     * @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 real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = FastFourierTransformer.sample(f, min, max, n);
+        data[0] = 0.0;
+        return fst(data);
+    }
+
+    /**
+     * Transform the given real data set.
+     * <p>
+     * The formula is $ F_n = \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
+     *
+     * @param f the real data array to be transformed
+     * @return the real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform2(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = Math.sqrt(2.0 / f.length);
+        return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
+    }
+
+    /**
+     * Transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ F_n = \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
+     *
+     * @param f the function to be sampled and transformed
+     * @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 real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] transform2(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = FastFourierTransformer.sample(f, min, max, n);
+        data[0] = 0.0;
+        double scaling_coefficient = Math.sqrt(2.0 / n);
+        return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real data set.
+     * <p>
+     * The formula is $ f_k = (2/N) \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
+     *
+     * @param f the real data array to be inversely transformed
+     * @return the real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double scaling_coefficient = 2.0 / f.length;
+        return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ f_k = (2/N) \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
+     *
+     * @param f the function to be sampled and inversely transformed
+     * @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 real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        double data[] = FastFourierTransformer.sample(f, min, max, n);
+        data[0] = 0.0;
+        double scaling_coefficient = 2.0 / n;
+        return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
+    }
+
+    /**
+     * Inversely transform the given real data set.
+     * <p>
+     * The formula is $ f_k = \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
+     *
+     * @param f the real data array to be inversely transformed
+     * @return the real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform2(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        return transform2(f);
+    }
+
+    /**
+     * Inversely transform the given real function, sampled on the given interval.
+     * <p>
+     * The formula is $ f_k = \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
+     *
+     * @param f the function to be sampled and inversely transformed
+     * @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 real inversely transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    public double[] inversetransform2(
+        UnivariateRealFunction f, double min, double max, int n)
+        throws MathException, IllegalArgumentException {
+
+        return transform2(f, min, max, n);
+    }
+
+    /**
+     * Perform the FST algorithm (including inverse).
+     *
+     * @param f the real data array to be transformed
+     * @return the real transformed array
+     * @throws MathException if any math-related errors occur
+     * @throws IllegalArgumentException if any parameters are invalid
+     */
+    protected double[] fst(double f[]) throws MathException,
+        IllegalArgumentException {
+
+        double A, B, x[], F[] = new double[f.length];
+
+        FastFourierTransformer.verifyDataSet(f);
+        if (f[0] != 0.0) {
+            throw new IllegalArgumentException
+                ("The first element is not zero: " + f[0]);
+        }
+        int N = f.length;
+        if (N == 1) {       // trivial case
+            F[0] = 0.0;
+            return F;
+        }
+
+        // construct a new array and perform FFT on it
+        x = new double[N];
+        x[0] = 0.0;
+        x[N >> 1] = 2.0 * f[N >> 1];
+        for (int i = 1; i < (N >> 1); i++) {
+            A = Math.sin(i * Math.PI / N) * (f[i] + f[N-i]);
+            B = 0.5 * (f[i] - f[N-i]);
+            x[i] = A + B;
+            x[N-i] = A - B;
+        }
+        FastFourierTransformer transformer = new FastFourierTransformer();
+        Complex y[] = transformer.transform(x);
+
+        // reconstruct the FST result for the original array
+        F[0] = 0.0;
+        F[1] = 0.5 * y[0].getReal();
+        for (int i = 1; i < (N >> 1); i++) {
+            F[2*i] = -y[i].getImaginary();
+            F[2*i+1] = y[i].getReal() + F[2*i-1];
+        }
+
+        return F;
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/FastSineTransformer.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html (added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html Tue Jul  4 13:58:54 2006
@@ -0,0 +1,21 @@
+<html>
+<!--
+   Copyright 2006 The Apache Software Foundation
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+  -->
+    <!-- $Revision$ $Date$ -->
+    <body>
+     Implementations of transform methods, including Fast Fourier transforms.
+    </body>
+</html>

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/transform/package.html
------------------------------------------------------------------------------
    svn:executable = *

Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastCosineTransformerTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastCosineTransformerTest.java?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastCosineTransformerTest.java (added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastCosineTransformerTest.java Tue Jul  4 13:58:54 2006
@@ -0,0 +1,120 @@
+/*
+ * Copyright 2005 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.transform;
+
+import org.apache.commons.math.analysis.*;
+import org.apache.commons.math.MathException;
+import junit.framework.TestCase;
+
+/**
+ * Testcase for fast cosine transformer.
+ * <p>
+ * FCT algorithm is exact, the small tolerance number is used only
+ * to account for round-off errors.
+ * 
+ * @version $Revision$ $Date$ 
+ */
+public final class FastCosineTransformerTest extends TestCase {
+
+    /**
+     * Test of transformer for the ad hoc data.
+     */
+    public void testAdHocData() throws MathException {
+        FastCosineTransformer transformer = new FastCosineTransformer();
+        double result[], tolerance = 1E-12;
+
+        double x[] = { 0.0, 1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0 };
+        double y[] = { 172.0, -105.096569476353, 27.3137084989848,
+                      -12.9593152353742, 8.0, -5.78585076868676,
+                       4.68629150101524, -4.15826451958632, 4.0 };
+
+        result = transformer.transform(x);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(y[i], result[i], tolerance);
+        }
+
+        result = transformer.inversetransform(y);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(x[i], result[i], tolerance);
+        }
+
+        FastFourierTransformer.scaleArray(x, Math.sqrt(0.5 * (x.length-1)));
+
+        result = transformer.transform2(y);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(x[i], result[i], tolerance);
+        }
+
+        result = transformer.inversetransform2(x);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(y[i], result[i], tolerance);
+        }
+    }
+
+    /**
+     * Test of transformer for the sine function.
+     */
+    public void testSinFunction() throws MathException {
+        UnivariateRealFunction f = new SinFunction();
+        FastCosineTransformer transformer = new FastCosineTransformer();
+        double min, max, result[], tolerance = 1E-12; int N = 9;
+
+        double expected[] = { 0.0, 3.26197262739567, 0.0,
+                             -2.17958042710327, 0.0, -0.648846697642915,
+                              0.0, -0.433545502649478, 0.0 };
+        min = 0.0; max = 2.0 * Math.PI * N / (N-1);
+        result = transformer.transform(f, min, max, N);
+        for (int i = 0; i < N; i++) {
+            assertEquals(expected[i], result[i], tolerance);
+        }
+
+        min = -Math.PI; max = Math.PI * (N+1) / (N-1);
+        result = transformer.transform(f, min, max, N);
+        for (int i = 0; i < N; i++) {
+            assertEquals(-expected[i], result[i], tolerance);
+        }
+    }
+
+    /**
+     * Test of parameters for the transformer.
+     */
+    public void testParameters() throws Exception {
+        UnivariateRealFunction f = new SinFunction();
+        FastCosineTransformer transformer = new FastCosineTransformer();
+
+        try {
+            // bad interval
+            transformer.transform(f, 1, -1, 65);
+            fail("Expecting IllegalArgumentException - bad interval");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+        try {
+            // bad samples number
+            transformer.transform(f, -1, 1, 1);
+            fail("Expecting IllegalArgumentException - bad samples number");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+        try {
+            // bad samples number
+            transformer.transform(f, -1, 1, 64);
+            fail("Expecting IllegalArgumentException - bad samples number");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastCosineTransformerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java (added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java Tue Jul  4 13:58:54 2006
@@ -0,0 +1,141 @@
+/*
+ * Copyright 2005 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.transform;
+
+import org.apache.commons.math.analysis.*;
+import org.apache.commons.math.complex.*;
+import org.apache.commons.math.MathException;
+import junit.framework.TestCase;
+
+/**
+ * Testcase for fast Fourier transformer.
+ * <p>
+ * FFT algorithm is exact, the small tolerance number is used only
+ * to account for round-off errors.
+ * 
+ * @version $Revision$ $Date$ 
+ */
+public final class FastFourierTransformerTest extends TestCase {
+
+    /**
+     * Test of transformer for the ad hoc data taken from Mathematica.
+     */
+    public void testAdHocData() throws MathException {
+        FastFourierTransformer transformer = new FastFourierTransformer();
+        Complex result[]; double tolerance = 1E-12;
+
+        double x[] = {1.3, 2.4, 1.7, 4.1, 2.9, 1.7, 5.1, 2.7};
+        Complex y[] = {
+            new Complex(21.9, 0.0),
+            new Complex(-2.09497474683058, 1.91507575950825),
+            new Complex(-2.6, 2.7),
+            new Complex(-1.10502525316942, -4.88492424049175),
+            new Complex(0.1, 0.0),
+            new Complex(-1.10502525316942, 4.88492424049175),
+            new Complex(-2.6, -2.7),
+            new Complex(-2.09497474683058, -1.91507575950825)};
+
+        result = transformer.transform(x);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(y[i].getReal(), result[i].getReal(), tolerance);
+            assertEquals(y[i].getImaginary(), result[i].getImaginary(), tolerance);
+        }
+
+        result = transformer.inversetransform(y);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(x[i], result[i].getReal(), tolerance);
+            assertEquals(0.0, result[i].getImaginary(), tolerance);
+        }
+
+        double x2[] = {10.4, 21.6, 40.8, 13.6, 23.2, 32.8, 13.6, 19.2};
+        transformer.scaleArray(x2, 1.0 / Math.sqrt(x2.length));
+        Complex y2[] = y;
+
+        result = transformer.transform2(y2);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(x2[i], result[i].getReal(), tolerance);
+            assertEquals(0.0, result[i].getImaginary(), tolerance);
+        }
+
+        result = transformer.inversetransform2(x2);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(y2[i].getReal(), result[i].getReal(), tolerance);
+            assertEquals(y2[i].getImaginary(), result[i].getImaginary(), tolerance);
+        }
+    }
+
+    /**
+     * Test of transformer for the sine function.
+     */
+    public void testSinFunction() throws MathException {
+        UnivariateRealFunction f = new SinFunction();
+        FastFourierTransformer transformer = new FastFourierTransformer();
+        Complex result[]; int N = 1 << 8;
+        double min, max, tolerance = 1E-12;
+
+        min = 0.0; max = 2.0 * Math.PI;
+        result = transformer.transform(f, min, max, N);
+        assertEquals(0.0, result[1].getReal(), tolerance);
+        assertEquals(-(N >> 1), result[1].getImaginary(), tolerance);
+        assertEquals(0.0, result[N-1].getReal(), tolerance);
+        assertEquals(N >> 1, result[N-1].getImaginary(), tolerance);
+        for (int i = 0; i < N-1; i += (i == 0 ? 2 : 1)) {
+            assertEquals(0.0, result[i].getReal(), tolerance);
+            assertEquals(0.0, result[i].getImaginary(), tolerance);
+        }
+
+        min = -Math.PI; max = Math.PI;
+        result = transformer.inversetransform(f, min, max, N);
+        assertEquals(0.0, result[1].getReal(), tolerance);
+        assertEquals(-0.5, result[1].getImaginary(), tolerance);
+        assertEquals(0.0, result[N-1].getReal(), tolerance);
+        assertEquals(0.5, result[N-1].getImaginary(), tolerance);
+        for (int i = 0; i < N-1; i += (i == 0 ? 2 : 1)) {
+            assertEquals(0.0, result[i].getReal(), tolerance);
+            assertEquals(0.0, result[i].getImaginary(), tolerance);
+        }
+    }
+
+    /**
+     * Test of parameters for the transformer.
+     */
+    public void testParameters() throws Exception {
+        UnivariateRealFunction f = new SinFunction();
+        FastFourierTransformer transformer = new FastFourierTransformer();
+
+        try {
+            // bad interval
+            transformer.transform(f, 1, -1, 64);
+            fail("Expecting IllegalArgumentException - bad interval");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+        try {
+            // bad samples number
+            transformer.transform(f, -1, 1, 0);
+            fail("Expecting IllegalArgumentException - bad samples number");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+        try {
+            // bad samples number
+            transformer.transform(f, -1, 1, 100);
+            fail("Expecting IllegalArgumentException - bad samples number");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastFourierTransformerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastSineTransformerTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastSineTransformerTest.java?rev=419098&view=auto
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastSineTransformerTest.java (added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastSineTransformerTest.java Tue Jul  4 13:58:54 2006
@@ -0,0 +1,119 @@
+/*
+ * Copyright 2005 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.transform;
+
+import org.apache.commons.math.analysis.*;
+import org.apache.commons.math.MathException;
+import junit.framework.TestCase;
+
+/**
+ * Testcase for fast sine transformer.
+ * <p>
+ * FST algorithm is exact, the small tolerance number is used only
+ * to account for round-off errors.
+ * 
+ * @version $Revision$ $Date$ 
+ */
+public final class FastSineTransformerTest extends TestCase {
+
+    /**
+     * Test of transformer for the ad hoc data.
+     */
+    public void testAdHocData() throws MathException {
+        FastSineTransformer transformer = new FastSineTransformer();
+        double result[], tolerance = 1E-12;
+
+        double x[] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
+        double y[] = { 0.0, 20.1093579685034, -9.65685424949238,
+                       5.98642305066196, -4.0, 2.67271455167720,
+                      -1.65685424949238, 0.795649469518633 };
+
+        result = transformer.transform(x);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(y[i], result[i], tolerance);
+        }
+
+        result = transformer.inversetransform(y);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(x[i], result[i], tolerance);
+        }
+
+        FastFourierTransformer.scaleArray(x, Math.sqrt(x.length / 2.0));
+
+        result = transformer.transform2(y);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(x[i], result[i], tolerance);
+        }
+
+        result = transformer.inversetransform2(x);
+        for (int i = 0; i < result.length; i++) {
+            assertEquals(y[i], result[i], tolerance);
+        }
+    }
+
+    /**
+     * Test of transformer for the sine function.
+     */
+    public void testSinFunction() throws MathException {
+        UnivariateRealFunction f = new SinFunction();
+        FastSineTransformer transformer = new FastSineTransformer();
+        double min, max, result[], tolerance = 1E-12; int N = 1 << 8;
+
+        min = 0.0; max = 2.0 * Math.PI;
+        result = transformer.transform(f, min, max, N);
+        assertEquals(N >> 1, result[2], tolerance);
+        for (int i = 0; i < N; i += (i == 1 ? 2 : 1)) {
+            assertEquals(0.0, result[i], tolerance);
+        }
+
+        min = -Math.PI; max = Math.PI;
+        result = transformer.transform(f, min, max, N);
+        assertEquals(-(N >> 1), result[2], tolerance);
+        for (int i = 0; i < N; i += (i == 1 ? 2 : 1)) {
+            assertEquals(0.0, result[i], tolerance);
+        }
+    }
+
+    /**
+     * Test of parameters for the transformer.
+     */
+    public void testParameters() throws Exception {
+        UnivariateRealFunction f = new SinFunction();
+        FastSineTransformer transformer = new FastSineTransformer();
+
+        try {
+            // bad interval
+            transformer.transform(f, 1, -1, 64);
+            fail("Expecting IllegalArgumentException - bad interval");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+        try {
+            // bad samples number
+            transformer.transform(f, -1, 1, 0);
+            fail("Expecting IllegalArgumentException - bad samples number");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+        try {
+            // bad samples number
+            transformer.transform(f, -1, 1, 100);
+            fail("Expecting IllegalArgumentException - bad samples number");
+        } catch (IllegalArgumentException ex) {
+            // expected
+        }
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/transform/FastSineTransformerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Modified: jakarta/commons/proper/math/trunk/xdocs/changes.xml
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/xdocs/changes.xml?rev=419098&r1=419097&r2=419098&view=diff
==============================================================================
--- jakarta/commons/proper/math/trunk/xdocs/changes.xml (original)
+++ jakarta/commons/proper/math/trunk/xdocs/changes.xml Tue Jul  4 13:58:54 2006
@@ -59,6 +59,9 @@
         not return incorrect results for numbers with bad IEEE754 
         representations.
       </action>
+      <action dev="psteitz" type="update" issue=MATH-140 due-to="Xiaogang Zhang">
+        Added Fast Fourier Transform implementation.
+      </action>
     </release>
     <release version="1.1" date="2005-12-17"  
  description="This is a maintenance release containing bug fixes and enhancements.



---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message