    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.transform;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math3.analysis.FunctionUtils;
    022    import org.apache.commons.math3.analysis.UnivariateFunction;
    023    import org.apache.commons.math3.complex.Complex;
    024    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    025    import org.apache.commons.math3.exception.util.LocalizedFormats;
    026    import org.apache.commons.math3.util.ArithmeticUtils;
    027    import org.apache.commons.math3.util.FastMath;
    028
    029    /**
    030     * Implements the Fast Sine Transform for transformation of one-dimensional real
    031     * data sets. For reference, see James S. Walker, <em>Fast Fourier
    032     * Transforms</em>, chapter 3 (ISBN 0849371635).
    033     * <p>
    034     * There are several variants of the discrete sine transform. The present
    035     * implementation corresponds to DST-I, with various normalization conventions,
    036     * which are specified by the parameter {@link DstNormalization}.
    037     * <strong>It should be noted that regardless to the convention, the first
    038     * element of the dataset to be transformed must be zero.</strong>
    039     * <p>
    040     * DST-I is equivalent to DFT of an <em>odd extension</em> of the data series.
    041     * More precisely, if x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is the data set
    042     * to be sine transformed, the extended data set x<sub>0</sub><sup>&#35;</sup>,
    043     * &hellip;, x<sub>2N-1</sub><sup>&#35;</sup> is defined as follows
    044     * <ul>
    045     * <li>x<sub>0</sub><sup>&#35;</sup> = x<sub>0</sub> = 0,</li>
    046     * <li>x<sub>k</sub><sup>&#35;</sup> = x<sub>k</sub> if 1 &le; k &lt; N,</li>
    047     * <li>x<sub>N</sub><sup>&#35;</sup> = 0,</li>
    048     * <li>x<sub>k</sub><sup>&#35;</sup> = -x<sub>2N-k</sub> if N + 1 &le; k &lt;
    049     * 2N.</li>
    050     * </ul>
    051     * <p>
    052     * Then, the standard DST-I y<sub>0</sub>, &hellip;, y<sub>N-1</sub> of the real
    053     * data set x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is equal to <em>half</em>
    054     * of i (the pure imaginary number) times the N first elements of the DFT of the
    055     * extended data set x<sub>0</sub><sup>&#35;</sup>, &hellip;,
    056     * x<sub>2N-1</sub><sup>&#35;</sup> <br />
    057     * y<sub>n</sub> = (i / 2) &sum;<sub>k=0</sub><sup>2N-1</sup>
    058     * x<sub>k</sub><sup>&#35;</sup> exp[-2&pi;i nk / (2N)]
    059     * &nbsp;&nbsp;&nbsp;&nbsp;k = 0, &hellip;, N-1.
    060     * <p>
    061     * The present implementation of the discrete sine transform as a fast sine
    062     * transform requires the length of the data to be a power of two. Besides,
    063     * it implicitly assumes that the sampled function is odd. In particular, the
    064     * first element of the data set must be 0, which is enforced in
    065     * {@link #transform(UnivariateFunction, double, double, int, TransformType)},
    066     * after sampling.
    067     *
    068     * @version $Id: FastSineTransformer.java 1385310 2012-09-16 16:32:10Z tn $
    069     * @since 1.2
    070     */
    071    public class FastSineTransformer implements RealTransformer, Serializable {
    072
    073        /** Serializable version identifier. */
    074        static final long serialVersionUID = 20120211L;
    075
    076        /** The type of DST to be performed. */
    077        private final DstNormalization normalization;
    078
    079        /**
    080         * Creates a new instance of this class, with various normalization conventions.
    081         *
    082         * @param normalization the type of normalization to be applied to the transformed data
    083         */
    084        public FastSineTransformer(final DstNormalization normalization) {
    085            this.normalization = normalization;
    086        }
    087
    088        /**
    089         * {@inheritDoc}
    090         *
    091         * The first element of the specified data set is required to be {@code 0}.
    092         *
    093         * @throws MathIllegalArgumentException if the length of the data array is
    094         *   not a power of two, or the first element of the data array is not zero
    095         */
    096        public double[] transform(final double[] f, final TransformType type) {
    097            if (normalization == DstNormalization.ORTHOGONAL_DST_I) {
    098                final double s = FastMath.sqrt(2.0 / f.length);
    099                return TransformUtils.scaleArray(fst(f), s);
    100            }
    101            if (type == TransformType.FORWARD) {
    102                return fst(f);
    103            }
    104            final double s = 2.0 / f.length;
    105            return TransformUtils.scaleArray(fst(f), s);
    106        }
    107
    108        /**
    109         * {@inheritDoc}
    110         *
    111         * This implementation enforces {@code f(x) = 0.0} at {@code x = 0.0}.
    112         *
    113         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    114         *   if the lower bound is greater than, or equal to the upper bound
    115         * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
    116         *   if the number of sample points is negative
    117         * @throws MathIllegalArgumentException if the number of sample points is not a power of two
    118         */
    119        public double[] transform(final UnivariateFunction f,
    120            final double min, final double max, final int n,
    121            final TransformType type) {
    122
    123            final double[] data = FunctionUtils.sample(f, min, max, n);
    124            data[0] = 0.0;
    125            return transform(data, type);
    126        }
    127
    128        /**
    129         * Perform the FST algorithm (including inverse). The first element of the
    130         * data set is required to be {@code 0}.
    131         *
    132         * @param f the real data array to be transformed
    133         * @return the real transformed array
    134         * @throws MathIllegalArgumentException if the length of the data array is
    135         *   not a power of two, or the first element of the data array is not zero
    136         */
    137        protected double[] fst(double[] f) throws MathIllegalArgumentException {
    138
    139            final double[] transformed = new double[f.length];
    140
    141            if (!ArithmeticUtils.isPowerOfTwo(f.length)) {
    142                throw new MathIllegalArgumentException(
    143                        LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
    144                        Integer.valueOf(f.length));
    145            }
    146            if (f[0] != 0.0) {
    147                throw new MathIllegalArgumentException(
    148                        LocalizedFormats.FIRST_ELEMENT_NOT_ZERO,
    149                        Double.valueOf(f[0]));
    150            }
    151            final int n = f.length;
    152            if (n == 1) {       // trivial case
    153                transformed[0] = 0.0;
    154                return transformed;
    155            }
    156
    157            // construct a new array and perform FFT on it
    158            final double[] x = new double[n];
    159            x[0] = 0.0;
    160            x[n >> 1] = 2.0 * f[n >> 1];
    161            for (int i = 1; i < (n >> 1); i++) {
    162                final double a = FastMath.sin(i * FastMath.PI / n) * (f[i] + f[n - i]);
    163                final double b = 0.5 * (f[i] - f[n - i]);
    164                x[i]     = a + b;
    165                x[n - i] = a - b;
    166            }
    167            FastFourierTransformer transformer;
    168            transformer = new FastFourierTransformer(DftNormalization.STANDARD);
    169            Complex[] y = transformer.transform(x, TransformType.FORWARD);
    170
    171            // reconstruct the FST result for the original array
    172            transformed[0] = 0.0;
    173            transformed[1] = 0.5 * y[0].getReal();
    174            for (int i = 1; i < (n >> 1); i++) {
    175                transformed[2 * i]     = -y[i].getImaginary();
    176                transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1];
    177            }
    178
    179            return transformed;
    180        }
    181    }