    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 Cosine Transform for transformation of one-dimensional
    031     * real 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 cosine transform. The present
    035     * implementation corresponds to DCT-I, with various normalization conventions,
    036     * which are specified by the parameter {@link DctNormalization}.
    037     * <p>
    038     * DCT-I is equivalent to DFT of an <em>even extension</em> of the data series.
    039     * More precisely, if x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is the data set
    040     * to be cosine transformed, the extended data set
    041     * x<sub>0</sub><sup>&#35;</sup>, &hellip;, x<sub>2N-3</sub><sup>&#35;</sup>
    042     * is defined as follows
    043     * <ul>
    044     * <li>x<sub>k</sub><sup>&#35;</sup> = x<sub>k</sub> if 0 &le; k &lt; N,</li>
    045     * <li>x<sub>k</sub><sup>&#35;</sup> = x<sub>2N-2-k</sub>
    046     * if N &le; k &lt; 2N - 2.</li>
    047     * </ul>
    048     * <p>
    049     * Then, the standard DCT-I y<sub>0</sub>, &hellip;, y<sub>N-1</sub> of the real
    050     * data set x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is equal to <em>half</em>
    051     * of the N first elements of the DFT of the extended data set
    052     * x<sub>0</sub><sup>&#35;</sup>, &hellip;, x<sub>2N-3</sub><sup>&#35;</sup>
    053     * <br/>
    054     * y<sub>n</sub> = (1 / 2) &sum;<sub>k=0</sub><sup>2N-3</sup>
    055     * x<sub>k</sub><sup>&#35;</sup> exp[-2&pi;i nk / (2N - 2)]
    056     * &nbsp;&nbsp;&nbsp;&nbsp;k = 0, &hellip;, N-1.
    057     * <p>
    058     * The present implementation of the discrete cosine transform as a fast cosine
    059     * transform requires the length of the data set to be a power of two plus one
    060     * (N&nbsp;=&nbsp;2<sup>n</sup>&nbsp;+&nbsp;1). Besides, it implicitly assumes
    061     * that the sampled function is even.
    062     *
    063     * @version $Id: FastCosineTransformer.java 1385310 2012-09-16 16:32:10Z tn $
    064     * @since 1.2
    065     */
    066    public class FastCosineTransformer implements RealTransformer, Serializable {
    067
    068        /** Serializable version identifier. */
    069        static final long serialVersionUID = 20120212L;
    070
    071        /** The type of DCT to be performed. */
    072        private final DctNormalization normalization;
    073
    074        /**
    075         * Creates a new instance of this class, with various normalization
    076         * conventions.
    077         *
    078         * @param normalization the type of normalization to be applied to the
    079         * transformed data
    080         */
    081        public FastCosineTransformer(final DctNormalization normalization) {
    082            this.normalization = normalization;
    083        }
    084
    085        /**
    086         * {@inheritDoc}
    087         *
    088         * @throws MathIllegalArgumentException if the length of the data array is
    089         * not a power of two plus one
    090         */
    091        public double[] transform(final double[] f, final TransformType type)
    092          throws MathIllegalArgumentException {
    093            if (type == TransformType.FORWARD) {
    094                if (normalization == DctNormalization.ORTHOGONAL_DCT_I) {
    095                    final double s = FastMath.sqrt(2.0 / (f.length - 1));
    096                    return TransformUtils.scaleArray(fct(f), s);
    097                }
    098                return fct(f);
    099            }
    100            final double s2 = 2.0 / (f.length - 1);
    101            final double s1;
    102            if (normalization == DctNormalization.ORTHOGONAL_DCT_I) {
    103                s1 = FastMath.sqrt(s2);
    104            } else {
    105                s1 = s2;
    106            }
    107            return TransformUtils.scaleArray(fct(f), s1);
    108        }
    109
    110        /**
    111         * {@inheritDoc}
    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
    118         * not a power of two plus one
    119         */
    120        public double[] transform(final UnivariateFunction f,
    121            final double min, final double max, final int n,
    122            final TransformType type) throws MathIllegalArgumentException {
    123
    124            final double[] data = FunctionUtils.sample(f, min, max, n);
    125            return transform(data, type);
    126        }
    127
    128        /**
    129         * Perform the FCT algorithm (including inverse).
    130         *
    131         * @param f the real data array to be transformed
    132         * @return the real transformed array
    133         * @throws MathIllegalArgumentException if the length of the data array is
    134         * not a power of two plus one
    135         */
    136        protected double[] fct(double[] f)
    137            throws MathIllegalArgumentException {
    138
    139            final double[] transformed = new double[f.length];
    140
    141            final int n = f.length - 1;
    142            if (!ArithmeticUtils.isPowerOfTwo(n)) {
    143                throw new MathIllegalArgumentException(
    144                    LocalizedFormats.NOT_POWER_OF_TWO_PLUS_ONE,
    145                    Integer.valueOf(f.length));
    146            }
    147            if (n == 1) {       // trivial case
    148                transformed[0] = 0.5 * (f[0] + f[1]);
    149                transformed[1] = 0.5 * (f[0] - f[1]);
    150                return transformed;
    151            }
    152
    153            // construct a new array and perform FFT on it
    154            final double[] x = new double[n];
    155            x[0] = 0.5 * (f[0] + f[n]);
    156            x[n >> 1] = f[n >> 1];
    157            // temporary variable for transformed[1]
    158            double t1 = 0.5 * (f[0] - f[n]);
    159            for (int i = 1; i < (n >> 1); i++) {
    160                final double a = 0.5 * (f[i] + f[n - i]);
    161                final double b = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n - i]);
    162                final double c = FastMath.cos(i * FastMath.PI / n) * (f[i] - f[n - i]);
    163                x[i] = a - b;
    164                x[n - i] = a + b;
    165                t1 += c;
    166            }
    167            FastFourierTransformer transformer;
    168            transformer = new FastFourierTransformer(DftNormalization.STANDARD);
    169            Complex[] y = transformer.transform(x, TransformType.FORWARD);
    170
    171            // reconstruct the FCT result for the original array
    172            transformed[0] = y[0].getReal();
    173            transformed[1] = t1;
    174            for (int i = 1; i < (n >> 1); i++) {
    175                transformed[2 * i]     = y[i].getReal();
    176                transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary();
    177            }
    178            transformed[n] = y[n >> 1].getReal();
    179
    180            return transformed;
    181        }
    182    }