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    }