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    
024    /**
025     * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
026     * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine Transform</a>
027     * for transformation of one-dimensional data sets. For reference, see
028     * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
029     * <p>
030     * FCT is its own inverse, up to a multiplier depending on conventions.
031     * The equations are listed in the comments of the corresponding methods.</p>
032     * <p>
033     * Different from FFT and FST, FCT requires the length of data set to be
034     * power of 2 plus one. Users should especially pay attention to the
035     * function transformation on how this affects the sampling.</p>
036     * <p>As of version 2.0 this no longer implements Serializable</p>
037     *
038     * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
039     * @since 1.2
040     */
041    public class FastCosineTransformer implements RealTransformer {
042    
043        /**
044         * Construct a default transformer.
045         */
046        public FastCosineTransformer() {
047            super();
048        }
049    
050        /**
051         * Transform the given real data set.
052         * <p>
053         * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
054         *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
055         * </p>
056         *
057         * @param f the real data array to be transformed
058         * @return the real transformed array
059         * @throws IllegalArgumentException if any parameters are invalid
060         */
061        public double[] transform(double f[]) throws IllegalArgumentException {
062            return fct(f);
063        }
064    
065        /**
066         * Transform the given real function, sampled on the given interval.
067         * <p>
068         * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
069         *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
070         * </p>
071         *
072         * @param f the function to be sampled and transformed
073         * @param min the lower bound for the interval
074         * @param max the upper bound for the interval
075         * @param n the number of sample points
076         * @return the real transformed array
077         * @throws FunctionEvaluationException if function cannot be evaluated
078         * at some point
079         * @throws IllegalArgumentException if any parameters are invalid
080         */
081        public double[] transform(UnivariateRealFunction f,
082                                  double min, double max, int n)
083            throws FunctionEvaluationException, IllegalArgumentException {
084            double data[] = FastFourierTransformer.sample(f, min, max, n);
085            return fct(data);
086        }
087    
088        /**
089         * Transform the given real data set.
090         * <p>
091         * The formula is F<sub>n</sub> = &radic;(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
092         *                        &radic;(2/N) &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
093         * </p>
094         *
095         * @param f the real data array to be transformed
096         * @return the real transformed array
097         * @throws IllegalArgumentException if any parameters are invalid
098         */
099        public double[] transform2(double f[]) throws IllegalArgumentException {
100    
101            double scaling_coefficient = Math.sqrt(2.0 / (f.length-1));
102            return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
103        }
104    
105        /**
106         * Transform the given real function, sampled on the given interval.
107         * <p>
108         * The formula is F<sub>n</sub> = &radic;(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
109         *                        &radic;(2/N) &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
110         *
111         * </p>
112         *
113         * @param f the function to be sampled and transformed
114         * @param min the lower bound for the interval
115         * @param max the upper bound for the interval
116         * @param n the number of sample points
117         * @return the real transformed array
118         * @throws FunctionEvaluationException if function cannot be evaluated
119         * at some point
120         * @throws IllegalArgumentException if any parameters are invalid
121         */
122        public double[] transform2(UnivariateRealFunction f,
123                                   double min, double max, int n)
124            throws FunctionEvaluationException, IllegalArgumentException {
125    
126            double data[] = FastFourierTransformer.sample(f, min, max, n);
127            double scaling_coefficient = Math.sqrt(2.0 / (n-1));
128            return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
129        }
130    
131        /**
132         * Inversely transform the given real data set.
133         * <p>
134         * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
135         *                        (2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
136         * </p>
137         *
138         * @param f the real data array to be inversely transformed
139         * @return the real inversely transformed array
140         * @throws IllegalArgumentException if any parameters are invalid
141         */
142        public double[] inversetransform(double f[]) throws IllegalArgumentException {
143    
144            double scaling_coefficient = 2.0 / (f.length - 1);
145            return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
146        }
147    
148        /**
149         * Inversely transform the given real function, sampled on the given interval.
150         * <p>
151         * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
152         *                        (2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
153         * </p>
154         *
155         * @param f the function to be sampled and inversely transformed
156         * @param min the lower bound for the interval
157         * @param max the upper bound for the interval
158         * @param n the number of sample points
159         * @return the real inversely transformed array
160         * @throws FunctionEvaluationException if function cannot be evaluated
161         * at some point
162         * @throws IllegalArgumentException if any parameters are invalid
163         */
164        public double[] inversetransform(UnivariateRealFunction f,
165                                         double min, double max, int n)
166            throws FunctionEvaluationException, IllegalArgumentException {
167    
168            double data[] = FastFourierTransformer.sample(f, min, max, n);
169            double scaling_coefficient = 2.0 / (n - 1);
170            return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
171        }
172    
173        /**
174         * Inversely transform the given real data set.
175         * <p>
176         * The formula is f<sub>k</sub> = &radic;(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
177         *                        &radic;(2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
178         * </p>
179         *
180         * @param f the real data array to be inversely transformed
181         * @return the real inversely transformed array
182         * @throws IllegalArgumentException if any parameters are invalid
183         */
184        public double[] inversetransform2(double f[]) throws IllegalArgumentException {
185            return transform2(f);
186        }
187    
188        /**
189         * Inversely transform the given real function, sampled on the given interval.
190         * <p>
191         * The formula is f<sub>k</sub> = &radic;(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
192         *                        &radic;(2/N) &sum;<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(&pi; nk/N)
193         * </p>
194         *
195         * @param f the function to be sampled and inversely transformed
196         * @param min the lower bound for the interval
197         * @param max the upper bound for the interval
198         * @param n the number of sample points
199         * @return the real inversely transformed array
200         * @throws FunctionEvaluationException if function cannot be evaluated
201         * 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                        "{0} is not a power of 2 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 = Math.sin(i * Math.PI / n) * (f[i] - f[n-i]);
243                final double c = Math.cos(i * Math.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    }