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://mathworld.wolfram.com/FastFourierTransform.html">
26   * Fast Fourier Transform</a> for transformation of one-dimensional data sets.
27   * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897,
28   * chapter 6.
29   * <p>
30   * There are several conventions for the definition of FFT and inverse FFT,
31   * mainly on different coefficient and exponent. Here the equations are listed
32   * in the comments of the corresponding methods.</p>
33   * <p>
34   * We require the length of data set to be power of 2, this greatly simplifies
35   * and speeds up the code. Users can pad the data with zeros to meet this
36   * requirement. There are other flavors of FFT, for reference, see S. Winograd,
37   * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation,
38   * 32 (1978), 175 - 199.</p>
39   *
40   * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
41   * @since 1.2
42   */
43  public class FastFourierTransformer implements Serializable {
44  
45      /** serializable version identifier */
46      static final long serialVersionUID = 5138259215438106000L;
47  
48      /** array of the roots of unity */
49      private Complex omega[] = new Complex[0];
50  
51      /**
52       * |omegaCount| is the length of lasted computed omega[]. omegaCount
53       * is positive for forward transform and negative for inverse transform.
54       */
55      private int omegaCount = 0;
56  
57      /**
58       * Construct a default transformer.
59       */
60      public FastFourierTransformer() {
61          super();
62      }
63  
64      /**
65       * Transform the given real data set.
66       * <p>
67       * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
68       * </p>
69       * 
70       * @param f the real data array to be transformed
71       * @return the complex transformed array
72       * @throws MathException if any math-related errors occur
73       * @throws IllegalArgumentException if any parameters are invalid
74       */
75      public Complex[] transform(double f[]) throws MathException,
76          IllegalArgumentException {
77  
78          return fft(f, false);
79      }
80  
81      /**
82       * Transform the given real function, sampled on the given interval.
83       * <p>
84       * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
85       * </p>
86       * 
87       * @param f the function to be sampled and transformed
88       * @param min the lower bound for the interval
89       * @param max the upper bound for the interval
90       * @param n the number of sample points
91       * @return the complex transformed array
92       * @throws MathException if any math-related errors occur
93       * @throws IllegalArgumentException if any parameters are invalid
94       */
95      public Complex[] transform(
96          UnivariateRealFunction f, double min, double max, int n)
97          throws MathException, IllegalArgumentException {
98  
99          double data[] = sample(f, min, max, n);
100         return fft(data, false);
101     }
102 
103     /**
104      * Transform the given complex data set.
105      * <p>
106      * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
107      * </p>
108      * 
109      * @param f the complex data array to be transformed
110      * @return the complex transformed array
111      * @throws MathException if any math-related errors occur
112      * @throws IllegalArgumentException if any parameters are invalid
113      */
114     public Complex[] transform(Complex f[]) throws MathException,
115         IllegalArgumentException {
116 
117         computeOmega(f.length);
118         return fft(f);
119     }
120 
121     /**
122      * Transform the given real data set.
123      * <p>
124      * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
125      * </p>
126      * 
127      * @param f the real data array to be transformed
128      * @return the complex transformed array
129      * @throws MathException if any math-related errors occur
130      * @throws IllegalArgumentException if any parameters are invalid
131      */
132     public Complex[] transform2(double f[]) throws MathException,
133         IllegalArgumentException {
134 
135         double scaling_coefficient = 1.0 / Math.sqrt(f.length);
136         return scaleArray(fft(f, false), scaling_coefficient);
137     }
138 
139     /**
140      * Transform the given real function, sampled on the given interval.
141      * <p>
142      * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
143      * </p>
144      * 
145      * @param f the function to be sampled and transformed
146      * @param min the lower bound for the interval
147      * @param max the upper bound for the interval
148      * @param n the number of sample points
149      * @return the complex transformed array
150      * @throws MathException if any math-related errors occur
151      * @throws IllegalArgumentException if any parameters are invalid
152      */
153     public Complex[] transform2(
154         UnivariateRealFunction f, double min, double max, int n)
155         throws MathException, IllegalArgumentException {
156 
157         double data[] = sample(f, min, max, n);
158         double scaling_coefficient = 1.0 / Math.sqrt(n);
159         return scaleArray(fft(data, false), scaling_coefficient);
160     }
161 
162     /**
163      * Transform the given complex data set.
164      * <p>
165      * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
166      * </p>
167      * 
168      * @param f the complex data array to be transformed
169      * @return the complex transformed array
170      * @throws MathException if any math-related errors occur
171      * @throws IllegalArgumentException if any parameters are invalid
172      */
173     public Complex[] transform2(Complex f[]) throws MathException,
174         IllegalArgumentException {
175 
176         computeOmega(f.length);
177         double scaling_coefficient = 1.0 / Math.sqrt(f.length);
178         return scaleArray(fft(f), scaling_coefficient);
179     }
180 
181     /**
182      * Inversely transform the given real data set.
183      * <p>
184      * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
185      * </p>
186      * 
187      * @param f the real data array to be inversely transformed
188      * @return the complex inversely transformed array
189      * @throws MathException if any math-related errors occur
190      * @throws IllegalArgumentException if any parameters are invalid
191      */
192     public Complex[] inversetransform(double f[]) throws MathException,
193         IllegalArgumentException {
194 
195         double scaling_coefficient = 1.0 / f.length;
196         return scaleArray(fft(f, true), scaling_coefficient);
197     }
198 
199     /**
200      * Inversely transform the given real function, sampled on the given interval.
201      * <p>
202      * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_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 complex inversely transformed array
210      * @throws MathException if any math-related errors occur
211      * @throws IllegalArgumentException if any parameters are invalid
212      */
213     public Complex[] inversetransform(
214         UnivariateRealFunction f, double min, double max, int n)
215         throws MathException, IllegalArgumentException {
216 
217         double data[] = sample(f, min, max, n);
218         double scaling_coefficient = 1.0 / n;
219         return scaleArray(fft(data, true), scaling_coefficient);
220     }
221 
222     /**
223      * Inversely transform the given complex data set.
224      * <p>
225      * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
226      * </p>
227      * 
228      * @param f the complex data array to be inversely transformed
229      * @return the complex inversely transformed array
230      * @throws MathException if any math-related errors occur
231      * @throws IllegalArgumentException if any parameters are invalid
232      */
233     public Complex[] inversetransform(Complex f[]) throws MathException,
234         IllegalArgumentException {
235 
236         computeOmega(-f.length);    // pass negative argument
237         double scaling_coefficient = 1.0 / f.length;
238         return scaleArray(fft(f), scaling_coefficient);
239     }
240 
241     /**
242      * Inversely transform the given real data set.
243      * <p>
244      * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
245      * </p>
246      * 
247      * @param f the real data array to be inversely transformed
248      * @return the complex inversely transformed array
249      * @throws MathException if any math-related errors occur
250      * @throws IllegalArgumentException if any parameters are invalid
251      */
252     public Complex[] inversetransform2(double f[]) throws MathException,
253         IllegalArgumentException {
254 
255         double scaling_coefficient = 1.0 / Math.sqrt(f.length);
256         return scaleArray(fft(f, true), scaling_coefficient);
257     }
258 
259     /**
260      * Inversely transform the given real function, sampled on the given interval.
261      * <p>
262      * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
263      * </p>
264      * 
265      * @param f the function to be sampled and inversely transformed
266      * @param min the lower bound for the interval
267      * @param max the upper bound for the interval
268      * @param n the number of sample points
269      * @return the complex inversely transformed array
270      * @throws MathException if any math-related errors occur
271      * @throws IllegalArgumentException if any parameters are invalid
272      */
273     public Complex[] inversetransform2(
274         UnivariateRealFunction f, double min, double max, int n)
275         throws MathException, IllegalArgumentException {
276 
277         double data[] = sample(f, min, max, n);
278         double scaling_coefficient = 1.0 / Math.sqrt(n);
279         return scaleArray(fft(data, true), scaling_coefficient);
280     }
281 
282     /**
283      * Inversely transform the given complex data set.
284      * <p>
285      * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
286      * </p>
287      * 
288      * @param f the complex data array to be inversely transformed
289      * @return the complex inversely transformed array
290      * @throws MathException if any math-related errors occur
291      * @throws IllegalArgumentException if any parameters are invalid
292      */
293     public Complex[] inversetransform2(Complex f[]) throws MathException,
294         IllegalArgumentException {
295 
296         computeOmega(-f.length);    // pass negative argument
297         double scaling_coefficient = 1.0 / Math.sqrt(f.length);
298         return scaleArray(fft(f), scaling_coefficient);
299     }
300 
301     /**
302      * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
303      *
304      * @param f the real data array to be transformed
305      * @param isInverse the indicator of forward or inverse transform
306      * @return the complex transformed array
307      * @throws MathException if any math-related errors occur
308      * @throws IllegalArgumentException if any parameters are invalid
309      */
310     protected Complex[] fft(double f[], boolean isInverse) throws
311         MathException, IllegalArgumentException {
312 
313         verifyDataSet(f);
314         Complex F[] = new Complex[f.length];
315         if (f.length == 1) {
316             F[0] = new Complex(f[0], 0.0);
317             return F;
318         }
319 
320         // Rather than the naive real to complex conversion, pack 2N
321         // real numbers into N complex numbers for better performance.
322         int N = f.length >> 1;
323         Complex c[] = new Complex[N];
324         for (int i = 0; i < N; i++) {
325             c[i] = new Complex(f[2*i], f[2*i+1]);
326         }
327         computeOmega(isInverse ? -N : N);
328         Complex z[] = fft(c);
329 
330         // reconstruct the FFT result for the original array
331         computeOmega(isInverse ? -2*N : 2*N);
332         F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);
333         F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0);
334         for (int i = 1; i < N; i++) {
335             Complex A = z[N-i].conjugate();
336             Complex B = z[i].add(A);
337             Complex C = z[i].subtract(A);
338             Complex D = omega[i].multiply(Complex.I);
339             F[i] = B.subtract(C.multiply(D));
340             F[2*N-i] = F[i].conjugate();
341         }
342 
343         return scaleArray(F, 0.5);
344     }
345 
346     /**
347      * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
348      *
349      * @param data the complex data array to be transformed
350      * @return the complex transformed array
351      * @throws MathException if any math-related errors occur
352      * @throws IllegalArgumentException if any parameters are invalid
353      */
354     protected Complex[] fft(Complex data[]) throws MathException,
355         IllegalArgumentException {
356 
357         int i, j, k, m, N = data.length;
358         Complex A, B, C, D, E, F, z, f[] = new Complex[N];
359 
360         // initial simple cases
361         verifyDataSet(data);
362         if (N == 1) {
363             f[0] = data[0];
364             return f;
365         }
366         if (N == 2) {
367             f[0] = data[0].add(data[1]);
368             f[1] = data[0].subtract(data[1]);
369             return f;
370         }
371 
372         // permute original data array in bit-reversal order
373         j = 0;
374         for (i = 0; i < N; i++) {
375             f[i] = data[j];
376             k = N >> 1;
377             while (j >= k && k > 0) {
378                 j -= k; k >>= 1;
379             }
380             j += k;
381         }
382 
383         // the bottom base-4 round
384         for (i = 0; i < N; i += 4) {
385             A = f[i].add(f[i+1]);
386             B = f[i+2].add(f[i+3]);
387             C = f[i].subtract(f[i+1]);
388             D = f[i+2].subtract(f[i+3]);
389             E = C.add(D.multiply(Complex.I));
390             F = C.subtract(D.multiply(Complex.I));
391             f[i] = A.add(B);
392             f[i+2] = A.subtract(B);
393             // omegaCount indicates forward or inverse transform
394             f[i+1] = omegaCount < 0 ? E : F;
395             f[i+3] = omegaCount > 0 ? E : F;
396         }
397 
398         // iterations from bottom to top take O(N*logN) time
399         for (i = 4; i < N; i <<= 1) {
400             m = N / (i<<1);
401             for (j = 0; j < N; j += i<<1) {
402                 for (k = 0; k < i; k++) {
403                     z = f[i+j+k].multiply(omega[k*m]);
404                     f[i+j+k] = f[j+k].subtract(z);
405                     f[j+k] = f[j+k].add(z);
406                 }
407             }
408         }
409         return f;
410     }
411 
412     /**
413      * Calculate the n-th roots of unity.
414      * <p>
415      * The computed omega[] = { 1, w, w^2, ... w^(n-1) } where
416      * w = exp(-2 \pi i / n), i = sqrt(-1). Note n is positive for
417      * forward transform and negative for inverse transform. </p>
418      * 
419      * @param n the integer passed in
420      * @throws IllegalArgumentException if n = 0
421      */
422     protected void computeOmega(int n) throws IllegalArgumentException {
423         if (n == 0) {
424             throw new IllegalArgumentException
425                 ("Cannot compute 0-th root of unity, indefinite result.");
426         }
427         // avoid repetitive calculations
428         if (n == omegaCount) { return; }
429         if (n + omegaCount == 0) {
430             for (int i = 0; i < Math.abs(omegaCount); i++) {
431                 omega[i] = omega[i].conjugate();
432             }
433             omegaCount = n;
434             return;
435         }
436         // calculate everything from scratch
437         omega = new Complex[Math.abs(n)];
438         double t = 2.0 * Math.PI / n;
439         double cost = Math.cos(t);
440         double sint = Math.sin(t);
441         omega[0] = new Complex(1.0, 0.0);
442         for (int i = 1; i < Math.abs(n); i++) {
443             omega[i] = new Complex(
444                 omega[i-1].getReal() * cost + omega[i-1].getImaginary() * sint,
445                 omega[i-1].getImaginary() * cost - omega[i-1].getReal() * sint);
446         }
447         omegaCount = n;
448     }
449 
450     /**
451      * Sample the given univariate real function on the given interval.
452      * <p>
453      * The interval is divided equally into N sections and sample points
454      * are taken from min to max-(max-min)/N. Usually f(x) is periodic
455      * such that f(min) = f(max) (note max is not sampled), but we don't
456      * require that.</p>
457      *
458      * @param f the function to be sampled
459      * @param min the lower bound for the interval
460      * @param max the upper bound for the interval
461      * @param n the number of sample points
462      * @return the samples array
463      * @throws MathException if any math-related errors occur
464      * @throws IllegalArgumentException if any parameters are invalid
465      */
466     public static double[] sample(
467         UnivariateRealFunction f, double min, double max, int n)
468         throws MathException, IllegalArgumentException {
469 
470         if (n <= 0) {
471             throw new IllegalArgumentException("Number of samples not positive.");
472         }
473         verifyInterval(min, max);
474 
475         double s[] = new double[n];
476         double h = (max - min) / n;
477         for (int i = 0; i < n; i++) {
478             s[i] = f.value(min + i * h);
479         }
480         return s;
481     }
482 
483     /**
484      * Multiply every component in the given real array by the
485      * given real number. The change is made in place.
486      *
487      * @param f the real array to be scaled
488      * @param d the real scaling coefficient
489      * @return a reference to the scaled array
490      */
491     public static double[] scaleArray(double f[], double d) {
492         for (int i = 0; i < f.length; i++) {
493             f[i] *= d;
494         }
495         return f;
496     }
497 
498     /**
499      * Multiply every component in the given complex array by the
500      * given real number. The change is made in place.
501      *
502      * @param f the complex array to be scaled
503      * @param d the real scaling coefficient
504      * @return a reference to the scaled array
505      */
506     public static Complex[] scaleArray(Complex f[], double d) {
507         for (int i = 0; i < f.length; i++) {
508             f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary());
509         }
510         return f;
511     }
512 
513     /**
514      * Returns true if the argument is power of 2.
515      * 
516      * @param n the number to test
517      * @return true if the argument is power of 2
518      */
519     public static boolean isPowerOf2(long n) {
520         return (n > 0) && ((n & (n - 1)) == 0);
521     }
522 
523     /**
524      * Verifies that the data set has length of power of 2.
525      * 
526      * @param d the data array
527      * @throws IllegalArgumentException if array length is not power of 2
528      */
529     public static void verifyDataSet(double d[]) throws IllegalArgumentException {
530         if (!isPowerOf2(d.length)) {
531             throw new IllegalArgumentException
532                 ("Number of samples not power of 2, consider padding for fix.");
533         }       
534     }
535 
536     /**
537      * Verifies that the data set has length of power of 2.
538      * 
539      * @param o the data array
540      * @throws IllegalArgumentException if array length is not power of 2
541      */
542     public static void verifyDataSet(Object o[]) throws IllegalArgumentException {
543         if (!isPowerOf2(o.length)) {
544             throw new IllegalArgumentException
545                 ("Number of samples not power of 2, consider padding for fix.");
546         }       
547     }
548 
549     /**
550      * Verifies that the endpoints specify an interval.
551      * 
552      * @param lower lower endpoint
553      * @param upper upper endpoint
554      * @throws IllegalArgumentException if not interval
555      */
556     public static void verifyInterval(double lower, double upper) throws
557         IllegalArgumentException {
558 
559         if (lower >= upper) {
560             throw new IllegalArgumentException
561                 ("Endpoints do not specify an interval: [" + lower +
562                 ", " + upper + "]");
563         }       
564     }
565 }