    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 Cosine 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     * FCT 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     * Different from FFT and FST, FCT requires the length of data set to be
    036     * power of 2 plus one. Users should especially pay attention to the
    037     * function transformation on how this affects the sampling.</p>
    038     * <p>As of version 2.0 this no longer implements Serializable</p>
    039     *
    040     * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
    041     * @since 1.2
    042     */
    043    public class FastCosineTransformer implements RealTransformer {
    044
    045        /**
    046         * Construct a default transformer.
    047         */
    048        public FastCosineTransformer() {
    049            super();
    050        }
    051
    052        /**
    053         * Transform the given real data set.
    054         * <p>
    055         * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
    056         *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
    057         * </p>
    058         *
    059         * @param f the real data array to be transformed
    060         * @return the real transformed array
    061         * @throws IllegalArgumentException if any parameters are invalid
    062         */
    063        public double[] transform(double f[]) throws IllegalArgumentException {
    064            return fct(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> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
    071         *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
    072         * </p>
    073         *
    074         * @param f the function to be sampled and transformed
    075         * @param min the lower bound for the interval
    076         * @param max the upper bound for the interval
    077         * @param n the number of sample points
    078         * @return the real transformed array
    079         * @throws FunctionEvaluationException if function cannot be evaluated
    080         * at some point
    081         * @throws IllegalArgumentException if any parameters are invalid
    082         */
    083        public double[] transform(UnivariateRealFunction f,
    084                                  double min, double max, int n)
    085            throws FunctionEvaluationException, IllegalArgumentException {
    086            double data[] = FastFourierTransformer.sample(f, min, max, n);
    087            return fct(data);
    088        }
    089
    090        /**
    091         * Transform the given real data set.
    092         * <p>
    093         * The formula is F<sub>n</sub> = &radic;(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
    094         *                        &radic;(2/N) &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&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-1));
    104            return FastFourierTransformer.scaleArray(fct(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;(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
    111         *                        &radic;(2/N) &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
    112         *
    113         * </p>
    114         *
    115         * @param f the function to be sampled and transformed
    116         * @param min the lower bound for the interval
    117         * @param max the upper bound for the interval
    118         * @param n the number of sample points
    119         * @return the real transformed array
    120         * @throws FunctionEvaluationException if function cannot be evaluated
    121         * at some point
    122         * @throws IllegalArgumentException if any parameters are invalid
    123         */
    124        public double[] transform2(UnivariateRealFunction f,
    125                                   double min, double max, int n)
    126            throws FunctionEvaluationException, IllegalArgumentException {
    127
    128            double data[] = FastFourierTransformer.sample(f, min, max, n);
    129            double scaling_coefficient = FastMath.sqrt(2.0 / (n-1));
    130            return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
    131        }
    132
    133        /**
    134         * Inversely transform the given real data set.
    135         * <p>
    136         * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
    137         *                        (2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
    138         * </p>
    139         *
    140         * @param f the real data array to be inversely transformed
    141         * @return the real inversely transformed array
    142         * @throws IllegalArgumentException if any parameters are invalid
    143         */
    144        public double[] inversetransform(double f[]) throws IllegalArgumentException {
    145
    146            double scaling_coefficient = 2.0 / (f.length - 1);
    147            return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
    148        }
    149
    150        /**
    151         * Inversely transform the given real function, sampled on the given interval.
    152         * <p>
    153         * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
    154         *                        (2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
    155         * </p>
    156         *
    157         * @param f the function to be sampled and inversely transformed
    158         * @param min the lower bound for the interval
    159         * @param max the upper bound for the interval
    160         * @param n the number of sample points
    161         * @return the real inversely transformed array
    162         * @throws FunctionEvaluationException if function cannot be evaluated at some point
    163         * @throws IllegalArgumentException if any parameters are invalid
    164         */
    165        public double[] inversetransform(UnivariateRealFunction f,
    166                                         double min, double max, int n)
    167            throws FunctionEvaluationException, IllegalArgumentException {
    168
    169            double data[] = FastFourierTransformer.sample(f, min, max, n);
    170            double scaling_coefficient = 2.0 / (n - 1);
    171            return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
    172        }
    173
    174        /**
    175         * Inversely transform the given real data set.
    176         * <p>
    177         * The formula is f<sub>k</sub> = &radic;(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
    178         *                        &radic;(2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
    179         * </p>
    180         *
    181         * @param f the real data array to be inversely transformed
    182         * @return the real inversely transformed array
    183         * @throws IllegalArgumentException if any parameters are invalid
    184         */
    185        public double[] inversetransform2(double f[]) throws IllegalArgumentException {
    186            return transform2(f);
    187        }
    188
    189        /**
    190         * Inversely transform the given real function, sampled on the given interval.
    191         * <p>
    192         * The formula is f<sub>k</sub> = &radic;(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
    193         *                        &radic;(2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
    194         * </p>
    195         *
    196         * @param f the function to be sampled and inversely transformed
    197         * @param min the lower bound for the interval
    198         * @param max the upper bound for the interval
    199         * @param n the number of sample points
    200         * @return the real inversely transformed array
    201         * @throws FunctionEvaluationException if function cannot be evaluated at some point
    202         * @throws IllegalArgumentException if any parameters are invalid
    203         */
    204        public double[] inversetransform2(UnivariateRealFunction f,
    205                                          double min, double max, int n)
    206            throws FunctionEvaluationException, IllegalArgumentException {
    207
    208            return transform2(f, min, max, n);
    209        }
    210
    211        /**
    212         * Perform the FCT algorithm (including inverse).
    213         *
    214         * @param f the real data array to be transformed
    215         * @return the real transformed array
    216         * @throws IllegalArgumentException if any parameters are invalid
    217         */
    218        protected double[] fct(double f[])
    219            throws IllegalArgumentException {
    220
    221            final double transformed[] = new double[f.length];
    222
    223            final int n = f.length - 1;
    224            if (!FastFourierTransformer.isPowerOf2(n)) {
    225                throw MathRuntimeException.createIllegalArgumentException(
    226                        LocalizedFormats.NOT_POWER_OF_TWO_PLUS_ONE,
    227                        f.length);
    228            }
    229            if (n == 1) {       // trivial case
    230                transformed[0] = 0.5 * (f[0] + f[1]);
    231                transformed[1] = 0.5 * (f[0] - f[1]);
    232                return transformed;
    233            }
    234
    235            // construct a new array and perform FFT on it
    236            final double[] x = new double[n];
    237            x[0] = 0.5 * (f[0] + f[n]);
    238            x[n >> 1] = f[n >> 1];
    239            double t1 = 0.5 * (f[0] - f[n]);   // temporary variable for transformed[1]
    240            for (int i = 1; i < (n >> 1); i++) {
    241                final double a = 0.5 * (f[i] + f[n-i]);
    242                final double b = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n-i]);
    243                final double c = FastMath.cos(i * FastMath.PI / n) * (f[i] - f[n-i]);
    244                x[i] = a - b;
    245                x[n-i] = a + b;
    246                t1 += c;
    247            }
    248            FastFourierTransformer transformer = new FastFourierTransformer();
    249            Complex y[] = transformer.transform(x);
    250
    251            // reconstruct the FCT result for the original array
    252            transformed[0] = y[0].getReal();
    253            transformed[1] = t1;
    254            for (int i = 1; i < (n >> 1); i++) {
    255                transformed[2 * i]     = y[i].getReal();
    256                transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary();
    257            }
    258            transformed[n] = y[n >> 1].getReal();
    259
    260            return transformed;
    261        }
    262    }