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 java.io.Serializable;
20  import org.apache.commons.math.analysis.*;
21  import org.apache.commons.math.complex.*;
22  import org.apache.commons.math.MathException;
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   *
37   * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
38   * @since 1.2
39   */
40  public class FastCosineTransformer implements Serializable {
41  
42      /** serializable version identifier */
43      static final long serialVersionUID = -7673941545134707766L;
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_n = (1/2) [f_0 + (-1)^n f_N] +
56       *                        \Sigma_{k=0}^{N-1} f_k \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 MathException if any math-related errors occur
62       * @throws IllegalArgumentException if any parameters are invalid
63       */
64      public double[] transform(double f[]) throws MathException,
65          IllegalArgumentException {
66  
67          return fct(f);
68      }
69  
70      /**
71       * Transform the given real function, sampled on the given interval.
72       * <p>
73       * The formula is $ F_n = (1/2) [f_0 + (-1)^n f_N] +
74       *                        \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
75       * </p>
76       * 
77       * @param f the function to be sampled and transformed
78       * @param min the lower bound for the interval
79       * @param max the upper bound for the interval
80       * @param n the number of sample points
81       * @return the real transformed array
82       * @throws MathException if any math-related errors occur
83       * @throws IllegalArgumentException if any parameters are invalid
84       */
85      public double[] transform(
86          UnivariateRealFunction f, double min, double max, int n)
87          throws MathException, IllegalArgumentException {
88  
89          double data[] = FastFourierTransformer.sample(f, min, max, n);
90          return fct(data);
91      }
92  
93      /**
94       * Transform the given real data set.
95       * <p>
96       * The formula is $ F_n = \sqrt{1/2N} [f_0 + (-1)^n f_N] +
97       *                        \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
98       * </p>
99       * 
100      * @param f the real data array to be transformed
101      * @return the real transformed array
102      * @throws MathException if any math-related errors occur
103      * @throws IllegalArgumentException if any parameters are invalid
104      */
105     public double[] transform2(double f[]) throws MathException,
106         IllegalArgumentException {
107 
108         double scaling_coefficient = Math.sqrt(2.0 / (f.length-1));
109         return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
110     }
111 
112     /**
113      * Transform the given real function, sampled on the given interval.
114      * <p>
115      * The formula is $ F_n = \sqrt{1/2N} [f_0 + (-1)^n f_N] +
116      *                        \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \cos(\pi nk/N) $
117      *
118      * </p>
119      * 
120      * @param f the function to be sampled and transformed
121      * @param min the lower bound for the interval
122      * @param max the upper bound for the interval
123      * @param n the number of sample points
124      * @return the real transformed array
125      * @throws MathException if any math-related errors occur
126      * @throws IllegalArgumentException if any parameters are invalid
127      */
128     public double[] transform2(
129         UnivariateRealFunction f, double min, double max, int n)
130         throws MathException, IllegalArgumentException {
131 
132         double data[] = FastFourierTransformer.sample(f, min, max, n);
133         double scaling_coefficient = Math.sqrt(2.0 / (n-1));
134         return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
135     }
136 
137     /**
138      * Inversely transform the given real data set.
139      * <p>
140      * The formula is $ f_k = (1/N) [F_0 + (-1)^k F_N] +
141      *                        (2/N) \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
142      * </p>
143      * 
144      * @param f the real data array to be inversely transformed
145      * @return the real inversely transformed array
146      * @throws MathException if any math-related errors occur
147      * @throws IllegalArgumentException if any parameters are invalid
148      */
149     public double[] inversetransform(double f[]) throws MathException,
150         IllegalArgumentException {
151 
152         double scaling_coefficient = 2.0 / (f.length - 1);
153         return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient);
154     }
155 
156     /**
157      * Inversely transform the given real function, sampled on the given interval.
158      * <p>
159      * The formula is $ f_k = (1/N) [F_0 + (-1)^k F_N] +
160      *                        (2/N) \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
161      * </p>
162      * 
163      * @param f the function to be sampled and inversely transformed
164      * @param min the lower bound for the interval
165      * @param max the upper bound for the interval
166      * @param n the number of sample points
167      * @return the real inversely transformed array
168      * @throws MathException if any math-related errors occur
169      * @throws IllegalArgumentException if any parameters are invalid
170      */
171     public double[] inversetransform(
172         UnivariateRealFunction f, double min, double max, int n)
173         throws MathException, IllegalArgumentException {
174 
175         double data[] = FastFourierTransformer.sample(f, min, max, n);
176         double scaling_coefficient = 2.0 / (n - 1);
177         return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
178     }
179 
180     /**
181      * Inversely transform the given real data set.
182      * <p>
183      * The formula is $ f_k = \sqrt{1/2N} [F_0 + (-1)^k F_N] +
184      *                        \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
185      * </p>
186      * 
187      * @param f the real data array to be inversely transformed
188      * @return the real inversely transformed array
189      * @throws MathException if any math-related errors occur
190      * @throws IllegalArgumentException if any parameters are invalid
191      */
192     public double[] inversetransform2(double f[]) throws MathException,
193         IllegalArgumentException {
194 
195         return transform2(f);
196     }
197 
198     /**
199      * Inversely transform the given real function, sampled on the given interval.
200      * <p>
201      * The formula is $ f_k = \sqrt{1/2N} [F_0 + (-1)^k F_N] +
202      *                        \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \cos(\pi nk/N) $
203      * </p>
204      * 
205      * @param f the function to be sampled and inversely transformed
206      * @param min the lower bound for the interval
207      * @param max the upper bound for the interval
208      * @param n the number of sample points
209      * @return the real inversely transformed array
210      * @throws MathException if any math-related errors occur
211      * @throws IllegalArgumentException if any parameters are invalid
212      */
213     public double[] inversetransform2(
214         UnivariateRealFunction f, double min, double max, int n)
215         throws MathException, IllegalArgumentException {
216 
217         return transform2(f, min, max, n);
218     }
219 
220     /**
221      * Perform the FCT algorithm (including inverse).
222      *
223      * @param f the real data array to be transformed
224      * @return the real transformed array
225      * @throws MathException if any math-related errors occur
226      * @throws IllegalArgumentException if any parameters are invalid
227      */
228     protected double[] fct(double f[]) throws MathException,
229         IllegalArgumentException {
230 
231         double A, B, C, F1, x[], F[] = new double[f.length];
232 
233         int N = f.length - 1;
234         if (!FastFourierTransformer.isPowerOf2(N)) {
235             throw new IllegalArgumentException
236                 ("Number of samples not power of 2 plus one: " + f.length);
237         }
238         if (N == 1) {       // trivial case
239             F[0] = 0.5 * (f[0] + f[1]);
240             F[1] = 0.5 * (f[0] - f[1]);
241             return F;
242         }
243 
244         // construct a new array and perform FFT on it
245         x = new double[N];
246         x[0] = 0.5 * (f[0] + f[N]);
247         x[N >> 1] = f[N >> 1];
248         F1 = 0.5 * (f[0] - f[N]);   // temporary variable for F[1]
249         for (int i = 1; i < (N >> 1); i++) {
250             A = 0.5 * (f[i] + f[N-i]);
251             B = Math.sin(i * Math.PI / N) * (f[i] - f[N-i]);
252             C = Math.cos(i * Math.PI / N) * (f[i] - f[N-i]);
253             x[i] = A - B;
254             x[N-i] = A + B;
255             F1 += C;
256         }
257         FastFourierTransformer transformer = new FastFourierTransformer();
258         Complex y[] = transformer.transform(x);
259 
260         // reconstruct the FCT result for the original array
261         F[0] = y[0].getReal();
262         F[1] = F1;
263         for (int i = 1; i < (N >> 1); i++) {
264             F[2*i] = y[i].getReal();
265             F[2*i+1] = F[2*i-1] - y[i].getImaginary();
266         }
267         F[N] = y[N >> 1].getReal();
268 
269         return F;
270     }
271 }