org.apache.commons.math3.transform.FastSineTransformer Maven / Gradle / Ivy
Show all versions of virtdata-lib-realer Show documentation
/*
* 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 java.io.Serializable;
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;
/**
* Implements the Fast Sine Transform for transformation of one-dimensional real
* data sets. For reference, see James S. Walker, Fast Fourier
* Transforms, chapter 3 (ISBN 0849371635).
*
* 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}.
* It should be noted that regardless to the convention, the first
* element of the dataset to be transformed must be zero.
*
* DST-I is equivalent to DFT of an odd extension of the data series.
* More precisely, if x0, …, xN-1 is the data set
* to be sine transformed, the extended data set x0#,
* …, x2N-1# is defined as follows
*
* - x0# = x0 = 0,
* - xk# = xk if 1 ≤ k < N,
* - xN# = 0,
* - xk# = -x2N-k if N + 1 ≤ k <
* 2N.
*
*
* Then, the standard DST-I y0, …, yN-1 of the real
* data set x0, …, xN-1 is equal to half
* of i (the pure imaginary number) times the N first elements of the DFT of the
* extended data set x0#, …,
* x2N-1#
* yn = (i / 2) ∑k=02N-1
* xk# exp[-2πi nk / (2N)]
* k = 0, …, N-1.
*
* 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}
*
* 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}
*
* 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;
}
}