diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/transform/FastCosineTransformer.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/transform/FastCosineTransformer.java | 177 |
1 files changed, 177 insertions, 0 deletions
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; + } +} |