| /* |
| * 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.math4.transform; |
| |
| import java.io.Serializable; |
| |
| import org.apache.commons.numbers.complex.Complex; |
| import org.apache.commons.numbers.core.ArithmeticUtils; |
| import org.apache.commons.math4.analysis.FunctionUtils; |
| import org.apache.commons.math4.analysis.UnivariateFunction; |
| import org.apache.commons.math4.exception.MathIllegalArgumentException; |
| import org.apache.commons.math4.exception.util.LocalizedFormats; |
| import org.apache.commons.math4.util.FastMath; |
| |
| /** |
| * 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> |
| * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 1 ≤ k < N,</li> |
| * <li>x<sub>N</sub><sup>#</sup> = 0,</li> |
| * <li>x<sub>k</sub><sup>#</sup> = -x<sub>2N-k</sub> if N + 1 ≤ k < |
| * 2N.</li> |
| * </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} |
| * |
| * 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 |
| */ |
| @Override |
| 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.math4.exception.NonMonotonicSequenceException |
| * if the lower bound is greater than, or equal to the upper bound |
| * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException |
| * if the number of sample points is negative |
| * @throws MathIllegalArgumentException if the number of sample points is not a power of two |
| */ |
| @Override |
| 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; |
| } |
| } |
