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    }