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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.transform;
23  
24  import java.io.Serializable;
25  
26  import org.hipparchus.analysis.FunctionUtils;
27  import org.hipparchus.analysis.UnivariateFunction;
28  import org.hipparchus.complex.Complex;
29  import org.hipparchus.exception.MathIllegalArgumentException;
30  import org.hipparchus.exception.MathRuntimeException;
31  import org.hipparchus.util.ArithmeticUtils;
32  import org.hipparchus.util.FastMath;
33  import org.hipparchus.util.MathArrays;
34  import org.hipparchus.util.MathUtils;
35  
36  /**
37   * Implements the Fast Fourier Transform for transformation of one-dimensional
38   * real or complex data sets. For reference, see <em>Applied Numerical Linear
39   * Algebra</em>, ISBN 0898713897, chapter 6.
40   * <p>
41   * There are several variants of the discrete Fourier transform, with various
42   * normalization conventions, which are specified by the parameter
43   * {@link DftNormalization}.
44   * <p>
45   * The current implementation of the discrete Fourier transform as a fast
46   * Fourier transform requires the length of the data set to be a power of 2.
47   * This greatly simplifies and speeds up the code. Users can pad the data with
48   * zeros to meet this requirement. There are other flavors of FFT, for
49   * reference, see S. Winograd,
50   * <i>On computing the discrete Fourier transform</i>, Mathematics of
51   * Computation, 32 (1978), 175 - 199.
52   *
53   * @see DftNormalization
54   */
55  public class FastFourierTransformer implements Serializable {
56  
57      /** Serializable version identifier. */
58      private static final long serialVersionUID = 20120210L;
59  
60      /**
61       * {@code W_SUB_N_R[i]} is the real part of
62       * {@code exp(- 2 * i * pi / n)}:
63       * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}.
64       */
65      private static final double[] W_SUB_N_R =
66              {  0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1
67              , 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1
68              , 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1
69              , 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1
70              , 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1
71              , 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1
72              , 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1
73              , 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0
74              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
75              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
76              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
77              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
78              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
79              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
80              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
81              , 0x1.0p0, 0x1.0p0, 0x1.0p0 };
82  
83      /**
84       * {@code W_SUB_N_I[i]} is the imaginary part of
85       * {@code exp(- 2 * i * pi / n)}:
86       * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}.
87       */
88      private static final double[] W_SUB_N_I =
89              {  0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1
90              , -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5
91              , -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9
92              , -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13
93              , -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17
94              , -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21
95              , -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25
96              , -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29
97              , -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33
98              , -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37
99              , -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41
100             , -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45
101             , -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49
102             , -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53
103             , -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57
104             , -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 };
105 
106     /** The type of DFT to be performed. */
107     private final DftNormalization normalization;
108 
109     /**
110      * Creates a new instance of this class, with various normalization
111      * conventions.
112      *
113      * @param normalization the type of normalization to be applied to the
114      * transformed data
115      */
116     public FastFourierTransformer(final DftNormalization normalization) {
117         this.normalization = normalization;
118     }
119 
120     /**
121      * Performs identical index bit reversal shuffles on two arrays of identical
122      * size. Each element in the array is swapped with another element based on
123      * the bit-reversal of the index. For example, in an array with length 16,
124      * item at binary index 0011 (decimal 3) would be swapped with the item at
125      * binary index 1100 (decimal 12).
126      *
127      * @param a the first array to be shuffled
128      * @param b the second array to be shuffled
129      */
130     private static void bitReversalShuffle2(double[] a, double[] b) {
131         final int n = a.length;
132         assert b.length == n;
133         final int halfOfN = n >> 1;
134 
135         int j = 0;
136         for (int i = 0; i < n; i++) {
137             if (i < j) {
138                 // swap indices i & j
139                 double temp = a[i];
140                 a[i] = a[j];
141                 a[j] = temp;
142 
143                 temp = b[i];
144                 b[i] = b[j];
145                 b[j] = temp;
146             }
147 
148             int k = halfOfN;
149             while (k <= j && k > 0) {
150                 j -= k;
151                 k >>= 1;
152             }
153             j += k;
154         }
155     }
156 
157     /**
158      * Applies the proper normalization to the specified transformed data.
159      *
160      * @param dataRI the unscaled transformed data
161      * @param normalization the normalization to be applied
162      * @param type the type of transform (forward, inverse) which resulted in the specified data
163      */
164     private static void normalizeTransformedData(final double[][] dataRI,
165         final DftNormalization normalization, final TransformType type) {
166 
167         final double[] dataR = dataRI[0];
168         final double[] dataI = dataRI[1];
169         final int n = dataR.length;
170         assert dataI.length == n;
171 
172         switch (normalization) {
173             case STANDARD:
174                 if (type == TransformType.INVERSE) {
175                     final double scaleFactor = 1.0 / n;
176                     for (int i = 0; i < n; i++) {
177                         dataR[i] *= scaleFactor;
178                         dataI[i] *= scaleFactor;
179                     }
180                 }
181                 break;
182             case UNITARY:
183                 final double scaleFactor = 1.0 / FastMath.sqrt(n);
184                 for (int i = 0; i < n; i++) {
185                     dataR[i] *= scaleFactor;
186                     dataI[i] *= scaleFactor;
187                 }
188                 break;
189             default:
190                 // This should never occur in normal conditions. However this
191                 // clause has been added as a safeguard if other types of
192                 // normalizations are ever implemented, and the corresponding
193                 // test is forgotten in the present switch.
194                 throw MathRuntimeException.createInternalError();
195         }
196     }
197 
198     /**
199      * Computes the standard transform of the specified complex data. The
200      * computation is done in place. The input data is laid out as follows
201      * <ul>
202      *   <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point,</li>
203      *   <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data point.</li>
204      * </ul>
205      *
206      * @param dataRI the two dimensional array of real and imaginary parts of the data
207      * @param normalization the normalization to be applied to the transformed data
208      * @param type the type of transform (forward, inverse) to be performed
209      * @throws MathIllegalArgumentException if the number of rows of the specified
210      *   array is not two, or the array is not rectangular
211      * @throws MathIllegalArgumentException if the number of data points is not
212      *   a power of two
213      */
214     public static void transformInPlace(final double[][] dataRI,
215         final DftNormalization normalization, final TransformType type) {
216 
217         MathUtils.checkDimension(dataRI.length, 2);
218         final double[] dataR = dataRI[0];
219         final double[] dataI = dataRI[1];
220         MathArrays.checkEqualLength(dataR, dataI);
221 
222         final int n = dataR.length;
223         if (!ArithmeticUtils.isPowerOfTwo(n)) {
224             throw new MathIllegalArgumentException(LocalizedFFTFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
225                                                    Integer.valueOf(n));
226         }
227 
228         if (n == 1) {
229             return;
230         } else if (n == 2) {
231             final double srcR0 = dataR[0];
232             final double srcI0 = dataI[0];
233             final double srcR1 = dataR[1];
234             final double srcI1 = dataI[1];
235 
236             // X_0 = x_0 + x_1
237             dataR[0] = srcR0 + srcR1;
238             dataI[0] = srcI0 + srcI1;
239             // X_1 = x_0 - x_1
240             dataR[1] = srcR0 - srcR1;
241             dataI[1] = srcI0 - srcI1;
242 
243             normalizeTransformedData(dataRI, normalization, type);
244             return;
245         }
246 
247         bitReversalShuffle2(dataR, dataI);
248 
249         // Do 4-term DFT.
250         if (type == TransformType.INVERSE) {
251             for (int i0 = 0; i0 < n; i0 += 4) {
252                 final int i1 = i0 + 1;
253                 final int i2 = i0 + 2;
254                 final int i3 = i0 + 3;
255 
256                 final double srcR0 = dataR[i0];
257                 final double srcI0 = dataI[i0];
258                 final double srcR1 = dataR[i2];
259                 final double srcI1 = dataI[i2];
260                 final double srcR2 = dataR[i1];
261                 final double srcI2 = dataI[i1];
262                 final double srcR3 = dataR[i3];
263                 final double srcI3 = dataI[i3];
264 
265                 // 4-term DFT
266                 // X_0 = x_0 + x_1 + x_2 + x_3
267                 dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
268                 dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
269                 // X_1 = x_0 - x_2 + j * (x_3 - x_1)
270                 dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1);
271                 dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3);
272                 // X_2 = x_0 - x_1 + x_2 - x_3
273                 dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
274                 dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
275                 // X_3 = x_0 - x_2 + j * (x_1 - x_3)
276                 dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3);
277                 dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1);
278             }
279         } else {
280             for (int i0 = 0; i0 < n; i0 += 4) {
281                 final int i1 = i0 + 1;
282                 final int i2 = i0 + 2;
283                 final int i3 = i0 + 3;
284 
285                 final double srcR0 = dataR[i0];
286                 final double srcI0 = dataI[i0];
287                 final double srcR1 = dataR[i2];
288                 final double srcI1 = dataI[i2];
289                 final double srcR2 = dataR[i1];
290                 final double srcI2 = dataI[i1];
291                 final double srcR3 = dataR[i3];
292                 final double srcI3 = dataI[i3];
293 
294                 // 4-term DFT
295                 // X_0 = x_0 + x_1 + x_2 + x_3
296                 dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
297                 dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
298                 // X_1 = x_0 - x_2 + j * (x_3 - x_1)
299                 dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3);
300                 dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1);
301                 // X_2 = x_0 - x_1 + x_2 - x_3
302                 dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
303                 dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
304                 // X_3 = x_0 - x_2 + j * (x_1 - x_3)
305                 dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1);
306                 dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3);
307             }
308         }
309 
310         int lastN0 = 4;
311         int lastLogN0 = 2;
312         while (lastN0 < n) {
313             int n0 = lastN0 << 1;
314             int logN0 = lastLogN0 + 1;
315             double wSubN0R = W_SUB_N_R[logN0];
316             double wSubN0I = W_SUB_N_I[logN0];
317             if (type == TransformType.INVERSE) {
318                 wSubN0I = -wSubN0I;
319             }
320 
321             // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2).
322             for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) {
323                 int destOddStartIndex = destEvenStartIndex + lastN0;
324 
325                 double wSubN0ToRR = 1;
326                 double wSubN0ToRI = 0;
327 
328                 for (int r = 0; r < lastN0; r++) {
329                     double grR = dataR[destEvenStartIndex + r];
330                     double grI = dataI[destEvenStartIndex + r];
331                     double hrR = dataR[destOddStartIndex + r];
332                     double hrI = dataI[destOddStartIndex + r];
333 
334                     // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr
335                     dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI * hrI;
336                     dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI * hrR;
337                     // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr
338                     dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI * hrI);
339                     dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI * hrR);
340 
341                     // WsubN0ToR *= WsubN0R
342                     double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I;
343                     double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R;
344                     wSubN0ToRR = nextWsubN0ToRR;
345                     wSubN0ToRI = nextWsubN0ToRI;
346                 }
347             }
348 
349             lastN0 = n0;
350             lastLogN0 = logN0;
351         }
352 
353         normalizeTransformedData(dataRI, normalization, type);
354     }
355 
356     /**
357      * Returns the (forward, inverse) transform of the specified real data set.
358      *
359      * @param f the real data array to be transformed
360      * @param type the type of transform (forward, inverse) to be performed
361      * @return the complex transformed array
362      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
363      */
364     public Complex[] transform(final double[] f, final TransformType type) {
365         final double[][] dataRI = { f.clone(), new double[f.length] };
366         transformInPlace(dataRI, normalization, type);
367         return TransformUtils.createComplexArray(dataRI);
368     }
369 
370     /**
371      * Returns the (forward, inverse) transform of the specified real function,
372      * sampled on the specified interval.
373      *
374      * @param f the function to be sampled and transformed
375      * @param min the (inclusive) lower bound for the interval
376      * @param max the (exclusive) upper bound for the interval
377      * @param n the number of sample points
378      * @param type the type of transform (forward, inverse) to be performed
379      * @return the complex transformed array
380      * @throws org.hipparchus.exception.MathIllegalArgumentException
381      *   if the lower bound is greater than, or equal to the upper bound
382      * @throws org.hipparchus.exception.MathIllegalArgumentException
383      *   if the number of sample points {@code n} is negative
384      * @throws MathIllegalArgumentException if the number of sample points
385      *   {@code n} is not a power of two
386      */
387     public Complex[] transform(final UnivariateFunction f,
388                                final double min, final double max, final int n,
389                                final TransformType type) {
390 
391         final double[] data = FunctionUtils.sample(f, min, max, n);
392         return transform(data, type);
393     }
394 
395     /**
396      * Returns the (forward, inverse) transform of the specified complex data set.
397      *
398      * @param f the complex data array to be transformed
399      * @param type the type of transform (forward, inverse) to be performed
400      * @return the complex transformed array
401      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
402      */
403     public Complex[] transform(final Complex[] f, final TransformType type) {
404         final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);
405 
406         transformInPlace(dataRI, normalization, type);
407 
408         return TransformUtils.createComplexArray(dataRI);
409     }
410 
411 }