1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.transform;
18
19import org.apache.commons.math.FunctionEvaluationException;
20import org.apache.commons.math.MathRuntimeException;
21import org.apache.commons.math.analysis.UnivariateRealFunction;
22import org.apache.commons.math.complex.Complex;
23import org.apache.commons.math.exception.util.LocalizedFormats;
24import org.apache.commons.math.util.FastMath;
25
26/**
27 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
28 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine Transform</a>
29 * for transformation of one-dimensional data sets. For reference, see
30 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
31 * <p>
32 * FCT is its own inverse, up to a multiplier depending on conventions.
33 * The equations are listed in the comments of the corresponding methods.</p>
34 * <p>
35 * Different from FFT and FST, FCT requires the length of data set to be
36 * power of 2 plus one. Users should especially pay attention to the
37 * function transformation on how this affects the sampling.</p>
38 * <p>As of version 2.0 this no longer implements Serializable</p>
39 *
40 * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
41 * @since 1.2
42 */
43public class FastCosineTransformer implements RealTransformer {
44
45    /**
46     * Construct a default transformer.
47     */
48    public FastCosineTransformer() {
49        super();
50    }
51
52    /**
53     * Transform the given real data set.
54     * <p>
55     * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
56     *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
57     * </p>
58     *
59     * @param f the real data array to be transformed
60     * @return the real transformed array
61     * @throws IllegalArgumentException if any parameters are invalid
62     */
63    public double[] transform(double f[]) throws IllegalArgumentException {
64        return fct(f);
65    }
66
67    /**
68     * Transform the given real function, sampled on the given interval.
69     * <p>
70     * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
71     *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
72     * </p>
73     *
74     * @param f the function to be sampled and transformed
75     * @param min the lower bound for the interval
76     * @param max the upper bound for the interval
77     * @param n the number of sample points
78     * @return the real transformed array
79     * @throws FunctionEvaluationException if function cannot be evaluated
80     * at some point
81     * @throws IllegalArgumentException if any parameters are invalid
82     */
83    public double[] transform(UnivariateRealFunction f,
84                              double min, double max, int n)
85        throws FunctionEvaluationException, IllegalArgumentException {
86        double data[] = FastFourierTransformer.sample(f, min, max, n);
87        return fct(data);
88    }
89
90    /**
91     * Transform the given real data set.
92     * <p>
93     * The formula is F<sub>n</sub> = &radic;(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
94     *                        &radic;(2/N) &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
95     * </p>
96     *
97     * @param f the real data array to be transformed
98     * @return the real transformed array
99     * @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}
263