View Javadoc

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   */
17  package org.apache.commons.math.transform;
18  
19  import org.apache.commons.math.FunctionEvaluationException;
20  import org.apache.commons.math.MathRuntimeException;
21  import org.apache.commons.math.analysis.UnivariateRealFunction;
22  import org.apache.commons.math.complex.Complex;
23  
24  /**
25   * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
26   * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine Transform</a>
27   * for transformation of one-dimensional data sets. For reference, see
28   * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
29   * <p>
30   * FCT is its own inverse, up to a multiplier depending on conventions.
31   * The equations are listed in the comments of the corresponding methods.</p>
32   * <p>
33   * Different from FFT and FST, FCT requires the length of data set to be
34   * power of 2 plus one. Users should especially pay attention to the
35   * function transformation on how this affects the sampling.</p>
36   * <p>As of version 2.0 this no longer implements Serializable</p>
37   *
38   * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
39   * @since 1.2
40   */
41  public class FastCosineTransformer implements RealTransformer {
42  
43      /**
44       * Construct a default transformer.
45       */
46      public FastCosineTransformer() {
47          super();
48      }
49  
50      /**
51       * Transform the given real data set.
52       * <p>
53       * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
54       *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
55       * </p>
56       * 
57       * @param f the real data array to be transformed
58       * @return the real transformed array
59       * @throws IllegalArgumentException if any parameters are invalid
60       */
61      public double[] transform(double f[]) throws IllegalArgumentException {
62          return fct(f);
63      }
64  
65      /**
66       * Transform the given real function, sampled on the given interval.
67       * <p>
68       * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
69       *                        &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
70       * </p>
71       * 
72       * @param f the function to be sampled and transformed
73       * @param min the lower bound for the interval
74       * @param max the upper bound for the interval
75       * @param n the number of sample points
76       * @return the real transformed array
77       * @throws FunctionEvaluationException if function cannot be evaluated
78       * at some point
79       * @throws IllegalArgumentException if any parameters are invalid
80       */
81      public double[] transform(UnivariateRealFunction f,
82                                double min, double max, int n)
83          throws FunctionEvaluationException, IllegalArgumentException {
84          double data[] = FastFourierTransformer.sample(f, min, max, n);
85          return fct(data);
86      }
87  
88      /**
89       * Transform the given real data set.
90       * <p>
91       * The formula is F<sub>n</sub> = &radic;(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
92       *                        &radic;(2/N) &sum;<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(&pi; nk/N)
93       * </p>
94       * 
95       * @param f the real data array to be transformed
96       * @return the real transformed array
97       * @throws IllegalArgumentException if any parameters are invalid
98       */
99      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         double A, B, C, F1, x[], F[] = new double[f.length];
222 
223         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             F[0] = 0.5 * (f[0] + f[1]);
231             F[1] = 0.5 * (f[0] - f[1]);
232             return F;
233         }
234 
235         // construct a new array and perform FFT on it
236         x = new double[N];
237         x[0] = 0.5 * (f[0] + f[N]);
238         x[N >> 1] = f[N >> 1];
239         F1 = 0.5 * (f[0] - f[N]);   // temporary variable for F[1]
240         for (int i = 1; i < (N >> 1); i++) {
241             A = 0.5 * (f[i] + f[N-i]);
242             B = Math.sin(i * Math.PI / N) * (f[i] - f[N-i]);
243             C = Math.cos(i * Math.PI / N) * (f[i] - f[N-i]);
244             x[i] = A - B;
245             x[N-i] = A + B;
246             F1 += C;
247         }
248         FastFourierTransformer transformer = new FastFourierTransformer();
249         Complex y[] = transformer.transform(x);
250 
251         // reconstruct the FCT result for the original array
252         F[0] = y[0].getReal();
253         F[1] = F1;
254         for (int i = 1; i < (N >> 1); i++) {
255             F[2*i] = y[i].getReal();
256             F[2*i+1] = F[2*i-1] - y[i].getImaginary();
257         }
258         F[N] = y[N >> 1].getReal();
259 
260         return F;
261     }
262 }