    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.math.transform;
    018
    019    import org.apache.commons.math.FunctionEvaluationException;
    020    import org.apache.commons.math.MathRuntimeException;
    021    import org.apache.commons.math.analysis.UnivariateRealFunction;
    022    import org.apache.commons.math.complex.Complex;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.util.FastMath;
    025
    026    /**
    027     * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
    028     * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a>
    029     * for transformation of one-dimensional data sets. For reference, see
    030     * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
    031     * <p>
    032     * FST is its own inverse, up to a multiplier depending on conventions.
    033     * The equations are listed in the comments of the corresponding methods.</p>
    034     * <p>
    035     * Similar to FFT, we also require the length of data set to be power of 2.
    036     * In addition, the first element must be 0 and it's enforced in function
    037     * transformation after sampling.</p>
    038     * <p>As of version 2.0 this no longer implements Serializable</p>
    039     *
    040     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    041     * @since 1.2
    042     */
    043    public class FastSineTransformer implements RealTransformer {
    044
    045        /**
    046         * Construct a default transformer.
    047         */
    048        public FastSineTransformer() {
    049            super();
    050        }
    051
    052        /**
    053         * Transform the given real data set.
    054         * <p>
    055         * The formula is F<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
    056         * </p>
    057         *
    058         * @param f the real data array to be transformed
    059         * @return the real transformed array
    060         * @throws IllegalArgumentException if any parameters are invalid
    061         */
    062        public double[] transform(double f[])
    063            throws IllegalArgumentException {
    064            return fst(f);
    065        }
    066
    067        /**
    068         * Transform the given real function, sampled on the given interval.
    069         * <p>
    070         * The formula is F<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
    071         * </p>
    072         *
    073         * @param f the function to be sampled and transformed
    074         * @param min the lower bound for the interval
    075         * @param max the upper bound for the interval
    076         * @param n the number of sample points
    077         * @return the real transformed array
    078         * @throws FunctionEvaluationException if function cannot be evaluated
    079         * at some point
    080         * @throws IllegalArgumentException if any parameters are invalid
    081         */
    082        public double[] transform(UnivariateRealFunction f,
    083                                  double min, double max, int n)
    084            throws FunctionEvaluationException, IllegalArgumentException {
    085
    086            double data[] = FastFourierTransformer.sample(f, min, max, n);
    087            data[0] = 0.0;
    088            return fst(data);
    089        }
    090
    091        /**
    092         * Transform the given real data set.
    093         * <p>
    094         * The formula is F<sub>n</sub> = &radic;(2/N) &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
    095         * </p>
    096         *
    097         * @param f the real data array to be transformed
    098         * @return the real transformed array
    099         * @throws IllegalArgumentException if any parameters are invalid
    100         */
    101        public double[] transform2(double f[]) throws IllegalArgumentException {
    102
    103            double scaling_coefficient = FastMath.sqrt(2.0 / f.length);
    104            return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
    105        }
    106
    107        /**
    108         * Transform the given real function, sampled on the given interval.
    109         * <p>
    110         * The formula is F<sub>n</sub> = &radic;(2/N) &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
    111         * </p>
    112         *
    113         * @param f the function to be sampled and transformed
    114         * @param min the lower bound for the interval
    115         * @param max the upper bound for the interval
    116         * @param n the number of sample points
    117         * @return the real transformed array
    118         * @throws FunctionEvaluationException if function cannot be evaluated
    119         * at some point
    120         * @throws IllegalArgumentException if any parameters are invalid
    121         */
    122        public double[] transform2(
    123            UnivariateRealFunction f, double min, double max, int n)
    124            throws FunctionEvaluationException, IllegalArgumentException {
    125
    126            double data[] = FastFourierTransformer.sample(f, min, max, n);
    127            data[0] = 0.0;
    128            double scaling_coefficient = FastMath.sqrt(2.0 / n);
    129            return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
    130        }
    131
    132        /**
    133         * Inversely transform the given real data set.
    134         * <p>
    135         * The formula is f<sub>k</sub> = (2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
    136         * </p>
    137         *
    138         * @param f the real data array to be inversely transformed
    139         * @return the real inversely transformed array
    140         * @throws IllegalArgumentException if any parameters are invalid
    141         */
    142        public double[] inversetransform(double f[]) throws IllegalArgumentException {
    143
    144            double scaling_coefficient = 2.0 / f.length;
    145            return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
    146        }
    147
    148        /**
    149         * Inversely transform the given real function, sampled on the given interval.
    150         * <p>
    151         * The formula is f<sub>k</sub> = (2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
    152         * </p>
    153         *
    154         * @param f the function to be sampled and inversely transformed
    155         * @param min the lower bound for the interval
    156         * @param max the upper bound for the interval
    157         * @param n the number of sample points
    158         * @return the real inversely transformed array
    159         * @throws FunctionEvaluationException if function cannot be evaluated at some point
    160         * @throws IllegalArgumentException if any parameters are invalid
    161         */
    162        public double[] inversetransform(UnivariateRealFunction f, double min, double max, int n)
    163            throws FunctionEvaluationException, IllegalArgumentException {
    164
    165            double data[] = FastFourierTransformer.sample(f, min, max, n);
    166            data[0] = 0.0;
    167            double scaling_coefficient = 2.0 / n;
    168            return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
    169        }
    170
    171        /**
    172         * Inversely transform the given real data set.
    173         * <p>
    174         * The formula is f<sub>k</sub> = &radic;(2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
    175         * </p>
    176         *
    177         * @param f the real data array to be inversely transformed
    178         * @return the real inversely transformed array
    179         * @throws IllegalArgumentException if any parameters are invalid
    180         */
    181        public double[] inversetransform2(double f[]) throws IllegalArgumentException {
    182
    183            return transform2(f);
    184        }
    185
    186        /**
    187         * Inversely transform the given real function, sampled on the given interval.
    188         * <p>
    189         * The formula is f<sub>k</sub> = &radic;(2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
    190         * </p>
    191         *
    192         * @param f the function to be sampled and inversely transformed
    193         * @param min the lower bound for the interval
    194         * @param max the upper bound for the interval
    195         * @param n the number of sample points
    196         * @return the real inversely transformed array
    197         * @throws FunctionEvaluationException if function cannot be evaluated at some point
    198         * @throws IllegalArgumentException if any parameters are invalid
    199         */
    200        public double[] inversetransform2(UnivariateRealFunction f, double min, double max, int n)
    201            throws FunctionEvaluationException, IllegalArgumentException {
    202
    203            return transform2(f, min, max, n);
    204        }
    205
    206        /**
    207         * Perform the FST algorithm (including inverse).
    208         *
    209         * @param f the real data array to be transformed
    210         * @return the real transformed array
    211         * @throws IllegalArgumentException if any parameters are invalid
    212         */
    213        protected double[] fst(double f[]) throws IllegalArgumentException {
    214
    215            final double transformed[] = new double[f.length];
    216
    217            FastFourierTransformer.verifyDataSet(f);
    218            if (f[0] != 0.0) {
    219                throw MathRuntimeException.createIllegalArgumentException(
    220                        LocalizedFormats.FIRST_ELEMENT_NOT_ZERO,
    221                        f[0]);
    222            }
    223            final int n = f.length;
    224            if (n == 1) {       // trivial case
    225                transformed[0] = 0.0;
    226                return transformed;
    227            }
    228
    229            // construct a new array and perform FFT on it
    230            final double[] x = new double[n];
    231            x[0] = 0.0;
    232            x[n >> 1] = 2.0 * f[n >> 1];
    233            for (int i = 1; i < (n >> 1); i++) {
    234                final double a = FastMath.sin(i * FastMath.PI / n) * (f[i] + f[n-i]);
    235                final double b = 0.5 * (f[i] - f[n-i]);
    236                x[i]     = a + b;
    237                x[n - i] = a - b;
    238            }
    239            FastFourierTransformer transformer = new FastFourierTransformer();
    240            Complex y[] = transformer.transform(x);
    241
    242            // reconstruct the FST result for the original array
    243            transformed[0] = 0.0;
    244            transformed[1] = 0.5 * y[0].getReal();
    245            for (int i = 1; i < (n >> 1); i++) {
    246                transformed[2 * i]     = -y[i].getImaginary();
    247                transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1];
    248            }
    249
    250            return transformed;
    251        }
    252    }