diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/transform')
11 files changed, 1872 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/transform/DctNormalization.java b/src/main/java/org/apache/commons/math3/transform/DctNormalization.java new file mode 100644 index 0000000..2f38d04 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/DctNormalization.java @@ -0,0 +1,62 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +/** + * This enumeration defines the various types of normalizations that can be applied to discrete + * cosine transforms (DCT). The exact definition of these normalizations is detailed below. + * + * @see FastCosineTransformer + * @since 3.0 + */ +public enum DctNormalization { + /** + * Should be passed to the constructor of {@link FastCosineTransformer} to use the + * <em>standard</em> normalization convention. The standard DCT-I normalization convention is + * defined as follows + * + * <ul> + * <li>forward transform: y<sub>n</sub> = (1/2) [x<sub>0</sub> + + * (-1)<sup>n</sup>x<sub>N-1</sub>] + ∑<sub>k=1</sub><sup>N-2</sup> x<sub>k</sub> + * cos[π nk / (N - 1)], + * <li>inverse transform: x<sub>k</sub> = [1 / (N - 1)] [y<sub>0</sub> + + * (-1)<sup>k</sup>y<sub>N-1</sub>] + [2 / (N - 1)] ∑<sub>n=1</sub><sup>N-2</sup> + * y<sub>n</sub> cos[π nk / (N - 1)], + * </ul> + * + * where N is the size of the data sample. + */ + STANDARD_DCT_I, + + /** + * Should be passed to the constructor of {@link FastCosineTransformer} to use the + * <em>orthogonal</em> normalization convention. The orthogonal DCT-I normalization convention + * is defined as follows + * + * <ul> + * <li>forward transform: y<sub>n</sub> = [2(N - 1)]<sup>-1/2</sup> [x<sub>0</sub> + + * (-1)<sup>n</sup>x<sub>N-1</sub>] + [2 / (N - 1)]<sup>1/2</sup> + * ∑<sub>k=1</sub><sup>N-2</sup> x<sub>k</sub> cos[π nk / (N - 1)], + * <li>inverse transform: x<sub>k</sub> = [2(N - 1)]<sup>-1/2</sup> [y<sub>0</sub> + + * (-1)<sup>k</sup>y<sub>N-1</sub>] + [2 / (N - 1)]<sup>1/2</sup> + * ∑<sub>n=1</sub><sup>N-2</sup> y<sub>n</sub> cos[π nk / (N - 1)], + * </ul> + * + * which makes the transform orthogonal. N is the size of the data sample. + */ + ORTHOGONAL_DCT_I; +} diff --git a/src/main/java/org/apache/commons/math3/transform/DftNormalization.java b/src/main/java/org/apache/commons/math3/transform/DftNormalization.java new file mode 100644 index 0000000..db0eddb --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/DftNormalization.java @@ -0,0 +1,58 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +/** + * This enumeration defines the various types of normalizations that can be applied to discrete + * Fourier transforms (DFT). The exact definition of these normalizations is detailed below. + * + * @see FastFourierTransformer + * @since 3.0 + */ +public enum DftNormalization { + /** + * Should be passed to the constructor of {@link FastFourierTransformer} to use the + * <em>standard</em> normalization convention. This normalization convention is defined as + * follows + * + * <ul> + * <li>forward transform: y<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> x<sub>k</sub> + * exp(-2πi n k / N), + * <li>inverse transform: x<sub>k</sub> = N<sup>-1</sup> ∑<sub>n=0</sub><sup>N-1</sup> + * y<sub>n</sub> exp(2πi n k / N), + * </ul> + * + * where N is the size of the data sample. + */ + STANDARD, + + /** + * Should be passed to the constructor of {@link FastFourierTransformer} to use the + * <em>unitary</em> normalization convention. This normalization convention is defined as + * follows + * + * <ul> + * <li>forward transform: y<sub>n</sub> = (1 / √N) ∑<sub>k=0</sub><sup>N-1</sup> + * x<sub>k</sub> exp(-2πi n k / N), + * <li>inverse transform: x<sub>k</sub> = (1 / √N) ∑<sub>n=0</sub><sup>N-1</sup> + * y<sub>n</sub> exp(2πi n k / N), + * </ul> + * + * which makes the transform unitary. N is the size of the data sample. + */ + UNITARY; +} diff --git a/src/main/java/org/apache/commons/math3/transform/DstNormalization.java b/src/main/java/org/apache/commons/math3/transform/DstNormalization.java new file mode 100644 index 0000000..0aba09e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/DstNormalization.java @@ -0,0 +1,59 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +/** + * This enumeration defines the various types of normalizations that can be applied to discrete sine + * transforms (DST). The exact definition of these normalizations is detailed below. + * + * @see FastSineTransformer + * @since 3.0 + */ +public enum DstNormalization { + /** + * Should be passed to the constructor of {@link FastSineTransformer} to use the + * <em>standard</em> normalization convention. The standard DST-I normalization convention is + * defined as follows + * + * <ul> + * <li>forward transform: y<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> x<sub>k</sub> + * sin(π nk / N), + * <li>inverse transform: x<sub>k</sub> = (2 / N) ∑<sub>n=0</sub><sup>N-1</sup> + * y<sub>n</sub> sin(π nk / N), + * </ul> + * + * where N is the size of the data sample, and x<sub>0</sub> = 0. + */ + STANDARD_DST_I, + + /** + * Should be passed to the constructor of {@link FastSineTransformer} to use the + * <em>orthogonal</em> normalization convention. The orthogonal DCT-I normalization convention + * is defined as follows + * + * <ul> + * <li>Forward transform: y<sub>n</sub> = √(2 / N) ∑<sub>k=0</sub><sup>N-1</sup> + * x<sub>k</sub> sin(π nk / N), + * <li>Inverse transform: x<sub>k</sub> = √(2 / N) ∑<sub>n=0</sub><sup>N-1</sup> + * y<sub>n</sub> sin(π nk / N), + * </ul> + * + * which makes the transform orthogonal. N is the size of the data sample, and x<sub>0</sub> = + * 0. + */ + ORTHOGONAL_DST_I +} diff --git a/src/main/java/org/apache/commons/math3/transform/FastCosineTransformer.java b/src/main/java/org/apache/commons/math3/transform/FastCosineTransformer.java new file mode 100644 index 0000000..1e73187 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/FastCosineTransformer.java @@ -0,0 +1,177 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +import org.apache.commons.math3.analysis.FunctionUtils; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.complex.Complex; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.ArithmeticUtils; +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * Implements the Fast Cosine Transform for transformation of one-dimensional real data sets. For + * reference, see James S. Walker, <em>Fast Fourier Transforms</em>, chapter 3 (ISBN 0849371635). + * + * <p>There are several variants of the discrete cosine transform. The present implementation + * corresponds to DCT-I, with various normalization conventions, which are specified by the + * parameter {@link DctNormalization}. + * + * <p>DCT-I is equivalent to DFT of an <em>even extension</em> of the data series. More precisely, + * if x<sub>0</sub>, …, x<sub>N-1</sub> is the data set to be cosine transformed, the + * extended data set x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup> is + * defined as follows + * + * <ul> + * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 0 ≤ k < N, + * <li>x<sub>k</sub><sup>#</sup> = x<sub>2N-2-k</sub> if N ≤ k < 2N - 2. + * </ul> + * + * <p>Then, the standard DCT-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real data set + * x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em> of the N first elements of the + * DFT of the extended data set x<sub>0</sub><sup>#</sup>, …, + * x<sub>2N-3</sub><sup>#</sup> <br> + * y<sub>n</sub> = (1 / 2) ∑<sub>k=0</sub><sup>2N-3</sup> x<sub>k</sub><sup>#</sup> + * exp[-2πi nk / (2N - 2)] k = 0, …, N-1. + * + * <p>The present implementation of the discrete cosine transform as a fast cosine transform + * requires the length of the data set to be a power of two plus one + * (N = 2<sup>n</sup> + 1). Besides, it implicitly assumes that the sampled + * function is even. + * + * @since 1.2 + */ +public class FastCosineTransformer implements RealTransformer, Serializable { + + /** Serializable version identifier. */ + static final long serialVersionUID = 20120212L; + + /** The type of DCT to be performed. */ + private final DctNormalization normalization; + + /** + * Creates a new instance of this class, with various normalization conventions. + * + * @param normalization the type of normalization to be applied to the transformed data + */ + public FastCosineTransformer(final DctNormalization normalization) { + this.normalization = normalization; + } + + /** + * {@inheritDoc} + * + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + * plus one + */ + public double[] transform(final double[] f, final TransformType type) + throws MathIllegalArgumentException { + if (type == TransformType.FORWARD) { + if (normalization == DctNormalization.ORTHOGONAL_DCT_I) { + final double s = FastMath.sqrt(2.0 / (f.length - 1)); + return TransformUtils.scaleArray(fct(f), s); + } + return fct(f); + } + final double s2 = 2.0 / (f.length - 1); + final double s1; + if (normalization == DctNormalization.ORTHOGONAL_DCT_I) { + s1 = FastMath.sqrt(s2); + } else { + s1 = s2; + } + return TransformUtils.scaleArray(fct(f), s1); + } + + /** + * {@inheritDoc} + * + * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException if the lower bound + * is greater than, or equal to the upper bound + * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException if the number of + * sample points is negative + * @throws MathIllegalArgumentException if the number of sample points is not a power of two + * plus one + */ + public double[] transform( + final UnivariateFunction f, + final double min, + final double max, + final int n, + final TransformType type) + throws MathIllegalArgumentException { + + final double[] data = FunctionUtils.sample(f, min, max, n); + return transform(data, type); + } + + /** + * Perform the FCT algorithm (including inverse). + * + * @param f the real data array to be transformed + * @return the real transformed array + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + * plus one + */ + protected double[] fct(double[] f) throws MathIllegalArgumentException { + + final double[] transformed = new double[f.length]; + + final int n = f.length - 1; + if (!ArithmeticUtils.isPowerOfTwo(n)) { + throw new MathIllegalArgumentException( + LocalizedFormats.NOT_POWER_OF_TWO_PLUS_ONE, Integer.valueOf(f.length)); + } + if (n == 1) { // trivial case + transformed[0] = 0.5 * (f[0] + f[1]); + transformed[1] = 0.5 * (f[0] - f[1]); + return transformed; + } + + // construct a new array and perform FFT on it + final double[] x = new double[n]; + x[0] = 0.5 * (f[0] + f[n]); + x[n >> 1] = f[n >> 1]; + // temporary variable for transformed[1] + double t1 = 0.5 * (f[0] - f[n]); + for (int i = 1; i < (n >> 1); i++) { + final double a = 0.5 * (f[i] + f[n - i]); + final double b = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n - i]); + final double c = FastMath.cos(i * FastMath.PI / n) * (f[i] - f[n - i]); + x[i] = a - b; + x[n - i] = a + b; + t1 += c; + } + FastFourierTransformer transformer; + transformer = new FastFourierTransformer(DftNormalization.STANDARD); + Complex[] y = transformer.transform(x, TransformType.FORWARD); + + // reconstruct the FCT result for the original array + transformed[0] = y[0].getReal(); + transformed[1] = t1; + for (int i = 1; i < (n >> 1); i++) { + transformed[2 * i] = y[i].getReal(); + transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary(); + } + transformed[n] = y[n >> 1].getReal(); + + return transformed; + } +} diff --git a/src/main/java/org/apache/commons/math3/transform/FastFourierTransformer.java b/src/main/java/org/apache/commons/math3/transform/FastFourierTransformer.java new file mode 100644 index 0000000..116b21d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/FastFourierTransformer.java @@ -0,0 +1,749 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +import org.apache.commons.math3.analysis.FunctionUtils; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.complex.Complex; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.ArithmeticUtils; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; + +import java.io.Serializable; +import java.lang.reflect.Array; + +/** + * Implements the Fast Fourier Transform for transformation of one-dimensional real or complex data + * sets. For reference, see <em>Applied Numerical Linear Algebra</em>, ISBN 0898713897, chapter 6. + * + * <p>There are several variants of the discrete Fourier transform, with various normalization + * conventions, which are specified by the parameter {@link DftNormalization}. + * + * <p>The current implementation of the discrete Fourier transform as a fast Fourier transform + * requires the length of the data set to be a 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. + * + * @see DftNormalization + * @since 1.2 + */ +public class FastFourierTransformer implements Serializable { + + /** Serializable version identifier. */ + static final long serialVersionUID = 20120210L; + + /** + * {@code W_SUB_N_R[i]} is the real part of {@code exp(- 2 * i * pi / n)}: {@code W_SUB_N_R[i] = + * cos(2 * pi/ n)}, where {@code n = 2^i}. + */ + private static final double[] W_SUB_N_R = { + 0x1.0p0, + -0x1.0p0, + 0x1.1a62633145c07p-54, + 0x1.6a09e667f3bcdp-1, + 0x1.d906bcf328d46p-1, + 0x1.f6297cff75cbp-1, + 0x1.fd88da3d12526p-1, + 0x1.ff621e3796d7ep-1, + 0x1.ffd886084cd0dp-1, + 0x1.fff62169b92dbp-1, + 0x1.fffd8858e8a92p-1, + 0x1.ffff621621d02p-1, + 0x1.ffffd88586ee6p-1, + 0x1.fffff62161a34p-1, + 0x1.fffffd8858675p-1, + 0x1.ffffff621619cp-1, + 0x1.ffffffd885867p-1, + 0x1.fffffff62161ap-1, + 0x1.fffffffd88586p-1, + 0x1.ffffffff62162p-1, + 0x1.ffffffffd8858p-1, + 0x1.fffffffff6216p-1, + 0x1.fffffffffd886p-1, + 0x1.ffffffffff621p-1, + 0x1.ffffffffffd88p-1, + 0x1.fffffffffff62p-1, + 0x1.fffffffffffd9p-1, + 0x1.ffffffffffff6p-1, + 0x1.ffffffffffffep-1, + 0x1.fffffffffffffp-1, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0, + 0x1.0p0 + }; + + /** + * {@code W_SUB_N_I[i]} is the imaginary part of {@code exp(- 2 * i * pi / n)}: {@code + * W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}. + */ + private static final double[] W_SUB_N_I = { + 0x1.1a62633145c07p-52, + -0x1.1a62633145c07p-53, + -0x1.0p0, + -0x1.6a09e667f3bccp-1, + -0x1.87de2a6aea963p-2, + -0x1.8f8b83c69a60ap-3, + -0x1.917a6bc29b42cp-4, + -0x1.91f65f10dd814p-5, + -0x1.92155f7a3667ep-6, + -0x1.921d1fcdec784p-7, + -0x1.921f0fe670071p-8, + -0x1.921f8becca4bap-9, + -0x1.921faaee6472dp-10, + -0x1.921fb2aecb36p-11, + -0x1.921fb49ee4ea6p-12, + -0x1.921fb51aeb57bp-13, + -0x1.921fb539ecf31p-14, + -0x1.921fb541ad59ep-15, + -0x1.921fb5439d73ap-16, + -0x1.921fb544197ap-17, + -0x1.921fb544387bap-18, + -0x1.921fb544403c1p-19, + -0x1.921fb544422c2p-20, + -0x1.921fb54442a83p-21, + -0x1.921fb54442c73p-22, + -0x1.921fb54442cefp-23, + -0x1.921fb54442d0ep-24, + -0x1.921fb54442d15p-25, + -0x1.921fb54442d17p-26, + -0x1.921fb54442d18p-27, + -0x1.921fb54442d18p-28, + -0x1.921fb54442d18p-29, + -0x1.921fb54442d18p-30, + -0x1.921fb54442d18p-31, + -0x1.921fb54442d18p-32, + -0x1.921fb54442d18p-33, + -0x1.921fb54442d18p-34, + -0x1.921fb54442d18p-35, + -0x1.921fb54442d18p-36, + -0x1.921fb54442d18p-37, + -0x1.921fb54442d18p-38, + -0x1.921fb54442d18p-39, + -0x1.921fb54442d18p-40, + -0x1.921fb54442d18p-41, + -0x1.921fb54442d18p-42, + -0x1.921fb54442d18p-43, + -0x1.921fb54442d18p-44, + -0x1.921fb54442d18p-45, + -0x1.921fb54442d18p-46, + -0x1.921fb54442d18p-47, + -0x1.921fb54442d18p-48, + -0x1.921fb54442d18p-49, + -0x1.921fb54442d18p-50, + -0x1.921fb54442d18p-51, + -0x1.921fb54442d18p-52, + -0x1.921fb54442d18p-53, + -0x1.921fb54442d18p-54, + -0x1.921fb54442d18p-55, + -0x1.921fb54442d18p-56, + -0x1.921fb54442d18p-57, + -0x1.921fb54442d18p-58, + -0x1.921fb54442d18p-59, + -0x1.921fb54442d18p-60 + }; + + /** The type of DFT to be performed. */ + private final DftNormalization normalization; + + /** + * Creates a new instance of this class, with various normalization conventions. + * + * @param normalization the type of normalization to be applied to the transformed data + */ + public FastFourierTransformer(final DftNormalization normalization) { + this.normalization = normalization; + } + + /** + * Performs identical index bit reversal shuffles on two arrays of identical size. Each element + * in the array is swapped with another element based on the bit-reversal of the index. For + * example, in an array with length 16, item at binary index 0011 (decimal 3) would be swapped + * with the item at binary index 1100 (decimal 12). + * + * @param a the first array to be shuffled + * @param b the second array to be shuffled + */ + private static void bitReversalShuffle2(double[] a, double[] b) { + final int n = a.length; + assert b.length == n; + final int halfOfN = n >> 1; + + int j = 0; + for (int i = 0; i < n; i++) { + if (i < j) { + // swap indices i & j + double temp = a[i]; + a[i] = a[j]; + a[j] = temp; + + temp = b[i]; + b[i] = b[j]; + b[j] = temp; + } + + int k = halfOfN; + while (k <= j && k > 0) { + j -= k; + k >>= 1; + } + j += k; + } + } + + /** + * Applies the proper normalization to the specified transformed data. + * + * @param dataRI the unscaled transformed data + * @param normalization the normalization to be applied + * @param type the type of transform (forward, inverse) which resulted in the specified data + */ + private static void normalizeTransformedData( + final double[][] dataRI, + final DftNormalization normalization, + final TransformType type) { + + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + final int n = dataR.length; + assert dataI.length == n; + + switch (normalization) { + case STANDARD: + if (type == TransformType.INVERSE) { + final double scaleFactor = 1.0 / ((double) n); + for (int i = 0; i < n; i++) { + dataR[i] *= scaleFactor; + dataI[i] *= scaleFactor; + } + } + break; + case UNITARY: + final double scaleFactor = 1.0 / FastMath.sqrt(n); + for (int i = 0; i < n; i++) { + dataR[i] *= scaleFactor; + dataI[i] *= scaleFactor; + } + break; + default: + /* + * This should never occur in normal conditions. However this + * clause has been added as a safeguard if other types of + * normalizations are ever implemented, and the corresponding + * test is forgotten in the present switch. + */ + throw new MathIllegalStateException(); + } + } + + /** + * Computes the standard transform of the specified complex data. The computation is done in + * place. The input data is laid out as follows + * + * <ul> + * <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point, + * <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data point. + * </ul> + * + * @param dataRI the two dimensional array of real and imaginary parts of the data + * @param normalization the normalization to be applied to the transformed data + * @param type the type of transform (forward, inverse) to be performed + * @throws DimensionMismatchException if the number of rows of the specified array is not two, + * or the array is not rectangular + * @throws MathIllegalArgumentException if the number of data points is not a power of two + */ + public static void transformInPlace( + final double[][] dataRI, + final DftNormalization normalization, + final TransformType type) { + + if (dataRI.length != 2) { + throw new DimensionMismatchException(dataRI.length, 2); + } + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + if (dataR.length != dataI.length) { + throw new DimensionMismatchException(dataI.length, dataR.length); + } + + final int n = dataR.length; + if (!ArithmeticUtils.isPowerOfTwo(n)) { + throw new MathIllegalArgumentException( + LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, Integer.valueOf(n)); + } + + if (n == 1) { + return; + } else if (n == 2) { + final double srcR0 = dataR[0]; + final double srcI0 = dataI[0]; + final double srcR1 = dataR[1]; + final double srcI1 = dataI[1]; + + // X_0 = x_0 + x_1 + dataR[0] = srcR0 + srcR1; + dataI[0] = srcI0 + srcI1; + // X_1 = x_0 - x_1 + dataR[1] = srcR0 - srcR1; + dataI[1] = srcI0 - srcI1; + + normalizeTransformedData(dataRI, normalization, type); + return; + } + + bitReversalShuffle2(dataR, dataI); + + // Do 4-term DFT. + if (type == TransformType.INVERSE) { + for (int i0 = 0; i0 < n; i0 += 4) { + final int i1 = i0 + 1; + final int i2 = i0 + 2; + final int i3 = i0 + 3; + + final double srcR0 = dataR[i0]; + final double srcI0 = dataI[i0]; + final double srcR1 = dataR[i2]; + final double srcI1 = dataI[i2]; + final double srcR2 = dataR[i1]; + final double srcI2 = dataI[i1]; + final double srcR3 = dataR[i3]; + final double srcI3 = dataI[i3]; + + // 4-term DFT + // X_0 = x_0 + x_1 + x_2 + x_3 + dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3; + dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3; + // X_1 = x_0 - x_2 + j * (x_3 - x_1) + dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1); + dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3); + // X_2 = x_0 - x_1 + x_2 - x_3 + dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3; + dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3; + // X_3 = x_0 - x_2 + j * (x_1 - x_3) + dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3); + dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1); + } + } else { + for (int i0 = 0; i0 < n; i0 += 4) { + final int i1 = i0 + 1; + final int i2 = i0 + 2; + final int i3 = i0 + 3; + + final double srcR0 = dataR[i0]; + final double srcI0 = dataI[i0]; + final double srcR1 = dataR[i2]; + final double srcI1 = dataI[i2]; + final double srcR2 = dataR[i1]; + final double srcI2 = dataI[i1]; + final double srcR3 = dataR[i3]; + final double srcI3 = dataI[i3]; + + // 4-term DFT + // X_0 = x_0 + x_1 + x_2 + x_3 + dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3; + dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3; + // X_1 = x_0 - x_2 + j * (x_3 - x_1) + dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3); + dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1); + // X_2 = x_0 - x_1 + x_2 - x_3 + dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3; + dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3; + // X_3 = x_0 - x_2 + j * (x_1 - x_3) + dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1); + dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3); + } + } + + int lastN0 = 4; + int lastLogN0 = 2; + while (lastN0 < n) { + int n0 = lastN0 << 1; + int logN0 = lastLogN0 + 1; + double wSubN0R = W_SUB_N_R[logN0]; + double wSubN0I = W_SUB_N_I[logN0]; + if (type == TransformType.INVERSE) { + wSubN0I = -wSubN0I; + } + + // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2). + for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) { + int destOddStartIndex = destEvenStartIndex + lastN0; + + double wSubN0ToRR = 1; + double wSubN0ToRI = 0; + + for (int r = 0; r < lastN0; r++) { + double grR = dataR[destEvenStartIndex + r]; + double grI = dataI[destEvenStartIndex + r]; + double hrR = dataR[destOddStartIndex + r]; + double hrI = dataI[destOddStartIndex + r]; + + // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr + dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI * hrI; + dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI * hrR; + // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr + dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI * hrI); + dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI * hrR); + + // WsubN0ToR *= WsubN0R + double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I; + double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R; + wSubN0ToRR = nextWsubN0ToRR; + wSubN0ToRI = nextWsubN0ToRI; + } + } + + lastN0 = n0; + lastLogN0 = logN0; + } + + normalizeTransformedData(dataRI, normalization, type); + } + + /** + * Returns the (forward, inverse) transform of the specified real data set. + * + * @param f the real data array to be transformed + * @param type the type of transform (forward, inverse) to be performed + * @return the complex transformed array + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + */ + public Complex[] transform(final double[] f, final TransformType type) { + final double[][] dataRI = + new double[][] {MathArrays.copyOf(f, f.length), new double[f.length]}; + + transformInPlace(dataRI, normalization, type); + + return TransformUtils.createComplexArray(dataRI); + } + + /** + * Returns the (forward, inverse) transform of the specified real function, sampled on the + * specified interval. + * + * @param f the function to be sampled and transformed + * @param min the (inclusive) lower bound for the interval + * @param max the (exclusive) upper bound for the interval + * @param n the number of sample points + * @param type the type of transform (forward, inverse) to be performed + * @return the complex transformed array + * @throws org.apache.commons.math3.exception.NumberIsTooLargeException if the lower bound is + * greater than, or equal to the upper bound + * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException if the number of + * sample points {@code n} is negative + * @throws MathIllegalArgumentException if the number of sample points {@code n} is not a power + * of two + */ + public Complex[] transform( + final UnivariateFunction f, + final double min, + final double max, + final int n, + final TransformType type) { + + final double[] data = FunctionUtils.sample(f, min, max, n); + return transform(data, type); + } + + /** + * Returns the (forward, inverse) transform of the specified complex data set. + * + * @param f the complex data array to be transformed + * @param type the type of transform (forward, inverse) to be performed + * @return the complex transformed array + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + */ + public Complex[] transform(final Complex[] f, final TransformType type) { + final double[][] dataRI = TransformUtils.createRealImaginaryArray(f); + + transformInPlace(dataRI, normalization, type); + + return TransformUtils.createComplexArray(dataRI); + } + + /** + * Performs a multi-dimensional Fourier transform on a given array. Use {@link + * #transform(Complex[], TransformType)} in a row-column implementation in any number of + * dimensions with O(N×log(N)) complexity with N = n<sub>1</sub> × n<sub>2</sub> + * ×n<sub>3</sub> × ... × n<sub>d</sub>, where n<sub>k</sub> is the number of + * elements in dimension k, and d is the total number of dimensions. + * + * @param mdca Multi-Dimensional Complex Array, i.e. {@code Complex[][][][]} + * @param type the type of transform (forward, inverse) to be performed + * @return transform of {@code mdca} as a Multi-Dimensional Complex Array, i.e. {@code + * Complex[][][][]} + * @throws IllegalArgumentException if any dimension is not a power of two + * @deprecated see MATH-736 + */ + @Deprecated + public Object mdfft(Object mdca, TransformType type) { + 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, type, i, new int[0]); + } + return mdcm.getArray(); + } + + /** + * Performs one dimension of a multi-dimensional Fourier transform. + * + * @param mdcm input matrix + * @param type the type of transform (forward, inverse) to be performed + * @param d index of the dimension to process + * @param subVector recursion subvector + * @throws IllegalArgumentException if any dimension is not a power of two + * @deprecated see MATH-736 + */ + @Deprecated + private void mdfft( + MultiDimensionalComplexMatrix mdcm, TransformType type, int d, int[] subVector) { + + 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); + } + + temp = transform(temp, type); + + 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, type, d, vector); + } else { + for (int i = 0; i < dimensionSize[subVector.length]; i++) { + vector[subVector.length] = i; + // further split along the next dimension + mdfft(mdcm, type, d, vector); + } + } + } + } + + /** + * Complex matrix implementation. 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. + * + * @deprecated see MATH-736 + */ + @Deprecated + private static class MultiDimensionalComplexMatrix implements Cloneable { + + /** Size in all dimensions. */ + protected int[] dimensionSize; + + /** Storage array. */ + protected Object multiDimensionalComplexArray; + + /** + * Simple constructor. + * + * @param multiDimensionalComplexArray array containing the matrix elements + */ + MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) { + + this.multiDimensionalComplexArray = multiDimensionalComplexArray; + + // count dimensions + int numOfDimensions = 0; + for (Object lastDimension = multiDimensionalComplexArray; + lastDimension instanceof Object[]; ) { + final Object[] array = (Object[]) lastDimension; + numOfDimensions++; + lastDimension = array[0]; + } + + // allocate array with exact count + dimensionSize = new int[numOfDimensions]; + + // fill array + numOfDimensions = 0; + for (Object lastDimension = multiDimensionalComplexArray; + lastDimension instanceof Object[]; ) { + final Object[] array = (Object[]) lastDimension; + dimensionSize[numOfDimensions++] = array.length; + lastDimension = array[0]; + } + } + + /** + * Get a matrix element. + * + * @param vector indices of the element + * @return matrix element + * @exception DimensionMismatchException if dimensions do not match + */ + public Complex get(int... vector) throws DimensionMismatchException { + + if (vector == null) { + if (dimensionSize.length > 0) { + throw new DimensionMismatchException(0, dimensionSize.length); + } + return null; + } + if (vector.length != dimensionSize.length) { + throw new DimensionMismatchException(vector.length, dimensionSize.length); + } + + Object lastDimension = multiDimensionalComplexArray; + + for (int i = 0; i < dimensionSize.length; i++) { + lastDimension = ((Object[]) lastDimension)[vector[i]]; + } + return (Complex) lastDimension; + } + + /** + * Set a matrix element. + * + * @param magnitude magnitude of the element + * @param vector indices of the element + * @return the previous value + * @exception DimensionMismatchException if dimensions do not match + */ + public Complex set(Complex magnitude, int... vector) throws DimensionMismatchException { + + if (vector == null) { + if (dimensionSize.length > 0) { + throw new DimensionMismatchException(0, dimensionSize.length); + } + return null; + } + if (vector.length != dimensionSize.length) { + throw new DimensionMismatchException(vector.length, dimensionSize.length); + } + + Object[] lastDimension = (Object[]) multiDimensionalComplexArray; + for (int i = 0; i < dimensionSize.length - 1; i++) { + lastDimension = (Object[]) lastDimension[vector[i]]; + } + + Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]]; + lastDimension[vector[dimensionSize.length - 1]] = magnitude; + + return lastValue; + } + + /** + * Get the size in all dimensions. + * + * @return size in all dimensions + */ + public int[] getDimensionSizes() { + return dimensionSize.clone(); + } + + /** + * Get the underlying storage array. + * + * @return underlying storage array + */ + public Object getArray() { + return multiDimensionalComplexArray; + } + + /** {@inheritDoc} */ + @Override + public Object clone() { + MultiDimensionalComplexMatrix mdcm = + new MultiDimensionalComplexMatrix( + Array.newInstance(Complex.class, dimensionSize)); + clone(mdcm); + return mdcm; + } + + /** + * Copy contents of current array into mdcm. + * + * @param mdcm array where to copy data + */ + 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/main/java/org/apache/commons/math3/transform/FastHadamardTransformer.java b/src/main/java/org/apache/commons/math3/transform/FastHadamardTransformer.java new file mode 100644 index 0000000..5f7e6a2 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/FastHadamardTransformer.java @@ -0,0 +1,316 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +import org.apache.commons.math3.analysis.FunctionUtils; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.ArithmeticUtils; + +import java.io.Serializable; + +/** + * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard + * Transform</a> (FHT). Transformation of an input vector x to the output vector y. + * + * <p>In addition to transformation of real vectors, the Hadamard transform can transform integer + * vectors into integer vectors. However, this integer transform cannot be inverted directly. Due to + * a scaling factor it may lead to rational results. As an example, the inverse transform of integer + * vector (0, 1, 0, 1) is rational vector (1/2, -1/2, 0, 0). + * + * @since 2.0 + */ +public class FastHadamardTransformer implements RealTransformer, Serializable { + + /** Serializable version identifier. */ + static final long serialVersionUID = 20120211L; + + /** + * {@inheritDoc} + * + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + */ + public double[] transform(final double[] f, final TransformType type) { + if (type == TransformType.FORWARD) { + return fht(f); + } + return TransformUtils.scaleArray(fht(f), 1.0 / f.length); + } + + /** + * {@inheritDoc} + * + * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException if the lower bound + * is greater than, or equal to the upper bound + * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException if the number of + * sample points is negative + * @throws MathIllegalArgumentException if the number of sample points is not a power of two + */ + public double[] transform( + final UnivariateFunction f, + final double min, + final double max, + final int n, + final TransformType type) { + + return transform(FunctionUtils.sample(f, min, max, n), type); + } + + /** + * Returns the forward transform of the specified integer data set.The integer transform cannot + * be inverted directly, due to a scaling factor which may lead to double results. + * + * @param f the integer data array to be transformed (signal) + * @return the integer transformed array (spectrum) + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + */ + public int[] transform(final int[] f) { + return fht(f); + } + + /** + * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. Requires + * {@code N * log2(N) = n * 2^n} additions. + * + * <h3>Short Table of manual calculation for N=8</h3> + * + * <ol> + * <li><b>x</b> is the input vector to be transformed, + * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>), + * <li>a and b are helper rows. + * </ol> + * + * <table align="center" border="1" cellpadding="3"> + * <tbody align="center"> + * <tr> + * <th>x</th> + * <th>a</th> + * <th>b</th> + * <th>y</th> + * </tr> + * <tr> + * <th>x<sub>0</sub></th> + * <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td> + * <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td> + * <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td> + * </tr> + * <tr> + * <th>x<sub>1</sub></th> + * <td>a<sub>1</sub> = x<sub>2</sub> + x<sub>3</sub></td> + * <td>b<sub>0</sub> = a<sub>2</sub> + a<sub>3</sub></td> + * <td>y<sub>0</sub> = b<sub>2</sub> + b<sub>3</sub></td> + * </tr> + * <tr> + * <th>x<sub>2</sub></th> + * <td>a<sub>2</sub> = x<sub>4</sub> + x<sub>5</sub></td> + * <td>b<sub>0</sub> = a<sub>4</sub> + a<sub>5</sub></td> + * <td>y<sub>0</sub> = b<sub>4</sub> + b<sub>5</sub></td> + * </tr> + * <tr> + * <th>x<sub>3</sub></th> + * <td>a<sub>3</sub> = x<sub>6</sub> + x<sub>7</sub></td> + * <td>b<sub>0</sub> = a<sub>6</sub> + a<sub>7</sub></td> + * <td>y<sub>0</sub> = b<sub>6</sub> + b<sub>7</sub></td> + * </tr> + * <tr> + * <th>x<sub>4</sub></th> + * <td>a<sub>0</sub> = x<sub>0</sub> - x<sub>1</sub></td> + * <td>b<sub>0</sub> = a<sub>0</sub> - a<sub>1</sub></td> + * <td>y<sub>0</sub> = b<sub>0</sub> - b<sub>1</sub></td> + * </tr> + * <tr> + * <th>x<sub>5</sub></th> + * <td>a<sub>1</sub> = x<sub>2</sub> - x<sub>3</sub></td> + * <td>b<sub>0</sub> = a<sub>2</sub> - a<sub>3</sub></td> + * <td>y<sub>0</sub> = b<sub>2</sub> - b<sub>3</sub></td> + * </tr> + * <tr> + * <th>x<sub>6</sub></th> + * <td>a<sub>2</sub> = x<sub>4</sub> - x<sub>5</sub></td> + * <td>b<sub>0</sub> = a<sub>4</sub> - a<sub>5</sub></td> + * <td>y<sub>0</sub> = b<sub>4</sub> - b<sub>5</sub></td> + * </tr> + * <tr> + * <th>x<sub>7</sub></th> + * <td>a<sub>3</sub> = x<sub>6</sub> - x<sub>7</sub></td> + * <td>b<sub>0</sub> = a<sub>6</sub> - a<sub>7</sub></td> + * <td>y<sub>0</sub> = b<sub>6</sub> - b<sub>7</sub></td> + * </tr> + * </tbody> + * </table> + * + * <h3>How it works</h3> + * + * <ol> + * <li>Construct a matrix with {@code N} rows and {@code n + 1} columns, {@code hadm[n+1][N]}. + * <br> + * <em>(If I use [x][y] it always means [row-offset][column-offset] of a Matrix with n + * rows and m columns. Its entries go from M[0][0] to M[n][N])</em> + * <li>Place the input vector {@code x[N]} in the first column of the matrix {@code hadm}. + * <li>The entries of the submatrix {@code D_top} are calculated as follows + * <ul> + * <li>{@code D_top} goes from entry {@code [0][1]} to {@code [N / 2 - 1][n + 1]}, + * <li>the columns of {@code D_top} are the pairwise mutually exclusive sums of the + * previous column. + * </ul> + * <li>The entries of the submatrix {@code D_bottom} are calculated as follows + * <ul> + * <li>{@code D_bottom} goes from entry {@code [N / 2][1]} to {@code [N][n + 1]}, + * <li>the columns of {@code D_bottom} are the pairwise differences of the previous + * column. + * </ul> + * <li>The consputation of {@code D_top} and {@code D_bottom} are best understood with the + * above example (for {@code N = 8}). + * <li>The output vector {@code y} is now in the last column of {@code hadm}. + * <li><em>Algorithm from <a + * href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">chipcenter</a>.</em> + * </ol> + * + * <h3>Visually</h3> + * + * <table border="1" align="center" cellpadding="3"> + * <tbody align="center"> + * <tr> + * <td></td><th>0</th><th>1</th><th>2</th><th>3</th> + * <th>…</th> + * <th>n + 1</th> + * </tr> + * <tr> + * <th>0</th> + * <td>x<sub>0</sub></td> + * <td colspan="5" rowspan="5" align="center" valign="middle"> + * ↑<br/> + * ← D<sub>top</sub> →<br/> + * ↓ + * </td> + * </tr> + * <tr><th>1</th><td>x<sub>1</sub></td></tr> + * <tr><th>2</th><td>x<sub>2</sub></td></tr> + * <tr><th>…</th><td>…</td></tr> + * <tr><th>N / 2 - 1</th><td>x<sub>N/2-1</sub></td></tr> + * <tr> + * <th>N / 2</th> + * <td>x<sub>N/2</sub></td> + * <td colspan="5" rowspan="5" align="center" valign="middle"> + * ↑<br/> + * ← D<sub>bottom</sub> →<br/> + * ↓ + * </td> + * </tr> + * <tr><th>N / 2 + 1</th><td>x<sub>N/2+1</sub></td></tr> + * <tr><th>N / 2 + 2</th><td>x<sub>N/2+2</sub></td></tr> + * <tr><th>…</th><td>…</td></tr> + * <tr><th>N</th><td>x<sub>N</sub></td></tr> + * </tbody> + * </table> + * + * @param x the real data array to be transformed + * @return the real transformed array, {@code y} + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + */ + protected double[] fht(double[] x) throws MathIllegalArgumentException { + + final int n = x.length; + final int halfN = n / 2; + + if (!ArithmeticUtils.isPowerOfTwo(n)) { + throw new MathIllegalArgumentException( + LocalizedFormats.NOT_POWER_OF_TWO, Integer.valueOf(n)); + } + + /* + * Instead of creating a matrix with p+1 columns and n rows, we use two + * one dimension arrays which we are used in an alternating way. + */ + double[] yPrevious = new double[n]; + double[] yCurrent = x.clone(); + + // iterate from left to right (column) + for (int j = 1; j < n; j <<= 1) { + + // switch columns + final double[] yTmp = yCurrent; + yCurrent = yPrevious; + yPrevious = yTmp; + + // iterate from top to bottom (row) + for (int i = 0; i < halfN; ++i) { + // Dtop: the top part works with addition + final int twoI = 2 * i; + yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; + } + for (int i = halfN; i < n; ++i) { + // Dbottom: the bottom part works with subtraction + final int twoI = 2 * i; + yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; + } + } + + return yCurrent; + } + + /** + * Returns the forward transform of the specified integer data set. The FHT (Fast Hadamard + * Transform) uses only subtraction and addition. + * + * @param x the integer data array to be transformed + * @return the integer transformed array, {@code y} + * @throws MathIllegalArgumentException if the length of the data array is not a power of two + */ + protected int[] fht(int[] x) throws MathIllegalArgumentException { + + final int n = x.length; + final int halfN = n / 2; + + if (!ArithmeticUtils.isPowerOfTwo(n)) { + throw new MathIllegalArgumentException( + LocalizedFormats.NOT_POWER_OF_TWO, Integer.valueOf(n)); + } + + /* + * Instead of creating a matrix with p+1 columns and n rows, we use two + * one dimension arrays which we are used in an alternating way. + */ + int[] yPrevious = new int[n]; + int[] yCurrent = x.clone(); + + // iterate from left to right (column) + for (int j = 1; j < n; j <<= 1) { + + // switch columns + final int[] yTmp = yCurrent; + yCurrent = yPrevious; + yPrevious = yTmp; + + // iterate from top to bottom (row) + for (int i = 0; i < halfN; ++i) { + // Dtop: the top part works with addition + final int twoI = 2 * i; + yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; + } + for (int i = halfN; i < n; ++i) { + // Dbottom: the bottom part works with subtraction + final int twoI = 2 * i; + yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; + } + } + + // return the last computed output vector y + return yCurrent; + } +} diff --git a/src/main/java/org/apache/commons/math3/transform/FastSineTransformer.java b/src/main/java/org/apache/commons/math3/transform/FastSineTransformer.java new file mode 100644 index 0000000..b33b226 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/FastSineTransformer.java @@ -0,0 +1,175 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +import org.apache.commons.math3.analysis.FunctionUtils; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.complex.Complex; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.ArithmeticUtils; +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * Implements the Fast Sine Transform for transformation of one-dimensional real data sets. For + * reference, see James S. Walker, <em>Fast Fourier Transforms</em>, chapter 3 (ISBN 0849371635). + * + * <p>There are several variants of the discrete sine transform. The present implementation + * corresponds to DST-I, with various normalization conventions, which are specified by the + * parameter {@link DstNormalization}. <strong>It should be noted that regardless to the convention, + * the first element of the dataset to be transformed must be zero.</strong> + * + * <p>DST-I is equivalent to DFT of an <em>odd extension</em> of the data series. More precisely, if + * x<sub>0</sub>, …, x<sub>N-1</sub> is the data set to be sine transformed, the extended + * data set x<sub>0</sub><sup>#</sup>, …, x<sub>2N-1</sub><sup>#</sup> is defined as + * follows + * + * <ul> + * <li>x<sub>0</sub><sup>#</sup> = x<sub>0</sub> = 0, + * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 1 ≤ k < N, + * <li>x<sub>N</sub><sup>#</sup> = 0, + * <li>x<sub>k</sub><sup>#</sup> = -x<sub>2N-k</sub> if N + 1 ≤ k < 2N. + * </ul> + * + * <p>Then, the standard DST-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real data set + * x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em> of i (the pure imaginary + * number) times the N first elements of the DFT of the extended data set + * x<sub>0</sub><sup>#</sup>, …, x<sub>2N-1</sub><sup>#</sup> <br> + * y<sub>n</sub> = (i / 2) ∑<sub>k=0</sub><sup>2N-1</sup> x<sub>k</sub><sup>#</sup> + * exp[-2πi nk / (2N)] k = 0, …, N-1. + * + * <p>The present implementation of the discrete sine transform as a fast sine transform requires + * the length of the data to be a power of two. Besides, it implicitly assumes that the sampled + * function is odd. In particular, the first element of the data set must be 0, which is enforced in + * {@link #transform(UnivariateFunction, double, double, int, TransformType)}, after sampling. + * + * @since 1.2 + */ +public class FastSineTransformer implements RealTransformer, Serializable { + + /** Serializable version identifier. */ + static final long serialVersionUID = 20120211L; + + /** The type of DST to be performed. */ + private final DstNormalization normalization; + + /** + * Creates a new instance of this class, with various normalization conventions. + * + * @param normalization the type of normalization to be applied to the transformed data + */ + public FastSineTransformer(final DstNormalization normalization) { + this.normalization = normalization; + } + + /** + * {@inheritDoc} + * + * <p>The first element of the specified data set is required to be {@code 0}. + * + * @throws MathIllegalArgumentException if the length of the data array is not a power of two, + * or the first element of the data array is not zero + */ + public double[] transform(final double[] f, final TransformType type) { + if (normalization == DstNormalization.ORTHOGONAL_DST_I) { + final double s = FastMath.sqrt(2.0 / f.length); + return TransformUtils.scaleArray(fst(f), s); + } + if (type == TransformType.FORWARD) { + return fst(f); + } + final double s = 2.0 / f.length; + return TransformUtils.scaleArray(fst(f), s); + } + + /** + * {@inheritDoc} + * + * <p>This implementation enforces {@code f(x) = 0.0} at {@code x = 0.0}. + * + * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException if the lower bound + * is greater than, or equal to the upper bound + * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException if the number of + * sample points is negative + * @throws MathIllegalArgumentException if the number of sample points is not a power of two + */ + public double[] transform( + final UnivariateFunction f, + final double min, + final double max, + final int n, + final TransformType type) { + + final double[] data = FunctionUtils.sample(f, min, max, n); + data[0] = 0.0; + return transform(data, type); + } + + /** + * Perform the FST algorithm (including inverse). The first element of the data set is required + * to be {@code 0}. + * + * @param f the real data array to be transformed + * @return the real transformed array + * @throws MathIllegalArgumentException if the length of the data array is not a power of two, + * or the first element of the data array is not zero + */ + protected double[] fst(double[] f) throws MathIllegalArgumentException { + + final double[] transformed = new double[f.length]; + + if (!ArithmeticUtils.isPowerOfTwo(f.length)) { + throw new MathIllegalArgumentException( + LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, Integer.valueOf(f.length)); + } + if (f[0] != 0.0) { + throw new MathIllegalArgumentException( + LocalizedFormats.FIRST_ELEMENT_NOT_ZERO, Double.valueOf(f[0])); + } + final int n = f.length; + if (n == 1) { // trivial case + transformed[0] = 0.0; + return transformed; + } + + // construct a new array and perform FFT on it + final double[] x = new double[n]; + x[0] = 0.0; + x[n >> 1] = 2.0 * f[n >> 1]; + for (int i = 1; i < (n >> 1); i++) { + final double a = FastMath.sin(i * FastMath.PI / n) * (f[i] + f[n - i]); + final double b = 0.5 * (f[i] - f[n - i]); + x[i] = a + b; + x[n - i] = a - b; + } + FastFourierTransformer transformer; + transformer = new FastFourierTransformer(DftNormalization.STANDARD); + Complex[] y = transformer.transform(x, TransformType.FORWARD); + + // reconstruct the FST result for the original array + transformed[0] = 0.0; + transformed[1] = 0.5 * y[0].getReal(); + for (int i = 1; i < (n >> 1); i++) { + transformed[2 * i] = -y[i].getImaginary(); + transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1]; + } + + return transformed; + } +} diff --git a/src/main/java/org/apache/commons/math3/transform/RealTransformer.java b/src/main/java/org/apache/commons/math3/transform/RealTransformer.java new file mode 100644 index 0000000..b5a105b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/RealTransformer.java @@ -0,0 +1,68 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; + +/** + * Interface for one-dimensional data sets transformations producing real results. + * + * <p>Such transforms include {@link FastSineTransformer sine transform}, {@link + * FastCosineTransformer cosine transform} or {@link FastHadamardTransformer Hadamard transform}. + * {@link FastFourierTransformer Fourier transform} is of a different kind and does not implement + * this interface since it produces {@link org.apache.commons.math3.complex.Complex} results instead + * of real ones. + * + * @since 2.0 + */ +public interface RealTransformer { + + /** + * Returns the (forward, inverse) transform of the specified real data set. + * + * @param f the real data array to be transformed (signal) + * @param type the type of transform (forward, inverse) to be performed + * @return the real transformed array (spectrum) + * @throws MathIllegalArgumentException if the array cannot be transformed with the given type + * (this may be for example due to array size, which is constrained in some transforms) + */ + double[] transform(double[] f, TransformType type) throws MathIllegalArgumentException; + + /** + * Returns the (forward, inverse) transform of the specified real function, sampled on the + * specified interval. + * + * @param f the function to be sampled and transformed + * @param min the (inclusive) lower bound for the interval + * @param max the (exclusive) upper bound for the interval + * @param n the number of sample points + * @param type the type of transform (forward, inverse) to be performed + * @return the real transformed array + * @throws NonMonotonicSequenceException if the lower bound is greater than, or equal to the + * upper bound + * @throws NotStrictlyPositiveException if the number of sample points is negative + * @throws MathIllegalArgumentException if the sample cannot be transformed with the given type + * (this may be for example due to sample size, which is constrained in some transforms) + */ + double[] transform(UnivariateFunction f, double min, double max, int n, TransformType type) + throws NonMonotonicSequenceException, + NotStrictlyPositiveException, + MathIllegalArgumentException; +} diff --git a/src/main/java/org/apache/commons/math3/transform/TransformType.java b/src/main/java/org/apache/commons/math3/transform/TransformType.java new file mode 100644 index 0000000..ae1a6c5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/TransformType.java @@ -0,0 +1,30 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +/** + * This enumeration defines the type of transform which is to be computed. + * + * @since 3.0 + */ +public enum TransformType { + /** The type to be specified for forward transforms. */ + FORWARD, + + /** The type to be specified for inverse transforms. */ + INVERSE; +} diff --git a/src/main/java/org/apache/commons/math3/transform/TransformUtils.java b/src/main/java/org/apache/commons/math3/transform/TransformUtils.java new file mode 100644 index 0000000..5cf83de --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/TransformUtils.java @@ -0,0 +1,160 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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.math3.transform; + +import org.apache.commons.math3.complex.Complex; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.util.LocalizedFormats; + +import java.util.Arrays; + +/** + * Useful functions for the implementation of various transforms. + * + * @since 3.0 + */ +public class TransformUtils { + /** + * Table of the powers of 2 to facilitate binary search lookup. + * + * @see #exactLog2(int) + */ + private static final int[] POWERS_OF_TWO = { + 0x00000001, 0x00000002, 0x00000004, 0x00000008, 0x00000010, 0x00000020, + 0x00000040, 0x00000080, 0x00000100, 0x00000200, 0x00000400, 0x00000800, + 0x00001000, 0x00002000, 0x00004000, 0x00008000, 0x00010000, 0x00020000, + 0x00040000, 0x00080000, 0x00100000, 0x00200000, 0x00400000, 0x00800000, + 0x01000000, 0x02000000, 0x04000000, 0x08000000, 0x10000000, 0x20000000, + 0x40000000 + }; + + /** Private constructor. */ + private TransformUtils() { + super(); + } + + /** + * 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; + } + + /** + * Builds a new two dimensional array of {@code double} filled with the real and imaginary parts + * of the specified {@link Complex} numbers. In the returned array {@code dataRI}, the data is + * laid out as follows + * + * <ul> + * <li>{@code dataRI[0][i] = dataC[i].getReal()}, + * <li>{@code dataRI[1][i] = dataC[i].getImaginary()}. + * </ul> + * + * @param dataC the array of {@link Complex} data to be transformed + * @return a two dimensional array filled with the real and imaginary parts of the specified + * complex input + */ + public static double[][] createRealImaginaryArray(final Complex[] dataC) { + final double[][] dataRI = new double[2][dataC.length]; + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + for (int i = 0; i < dataC.length; i++) { + final Complex c = dataC[i]; + dataR[i] = c.getReal(); + dataI[i] = c.getImaginary(); + } + return dataRI; + } + + /** + * Builds a new array of {@link Complex} from the specified two dimensional array of real and + * imaginary parts. In the returned array {@code dataC}, the data is laid out as follows + * + * <ul> + * <li>{@code dataC[i].getReal() = dataRI[0][i]}, + * <li>{@code dataC[i].getImaginary() = dataRI[1][i]}. + * </ul> + * + * @param dataRI the array of real and imaginary parts to be transformed + * @return an array of {@link Complex} with specified real and imaginary parts. + * @throws DimensionMismatchException if the number of rows of the specified array is not two, + * or the array is not rectangular + */ + public static Complex[] createComplexArray(final double[][] dataRI) + throws DimensionMismatchException { + + if (dataRI.length != 2) { + throw new DimensionMismatchException(dataRI.length, 2); + } + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + if (dataR.length != dataI.length) { + throw new DimensionMismatchException(dataI.length, dataR.length); + } + + final int n = dataR.length; + final Complex[] c = new Complex[n]; + for (int i = 0; i < n; i++) { + c[i] = new Complex(dataR[i], dataI[i]); + } + return c; + } + + /** + * Returns the base-2 logarithm of the specified {@code int}. Throws an exception if {@code n} + * is not a power of two. + * + * @param n the {@code int} whose base-2 logarithm is to be evaluated + * @return the base-2 logarithm of {@code n} + * @throws MathIllegalArgumentException if {@code n} is not a power of two + */ + public static int exactLog2(final int n) throws MathIllegalArgumentException { + + int index = Arrays.binarySearch(TransformUtils.POWERS_OF_TWO, n); + if (index < 0) { + throw new MathIllegalArgumentException( + LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, Integer.valueOf(n)); + } + return index; + } +} diff --git a/src/main/java/org/apache/commons/math3/transform/package-info.java b/src/main/java/org/apache/commons/math3/transform/package-info.java new file mode 100644 index 0000000..728c7c4 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/transform/package-info.java @@ -0,0 +1,18 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ +/** Implementations of transform methods, including Fast Fourier transforms. */ +package org.apache.commons.math3.transform; |