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