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.util.Random;
25  
26  import org.hipparchus.analysis.UnivariateFunction;
27  import org.hipparchus.analysis.function.Sin;
28  import org.hipparchus.analysis.function.Sinc;
29  import org.hipparchus.complex.Complex;
30  import org.hipparchus.exception.MathIllegalArgumentException;
31  import org.hipparchus.util.FastMath;
32  import org.junit.Assert;
33  import org.junit.Test;
34  
35  /**
36   * Test case for fast Fourier transformer.
37   * <p>
38   * FFT algorithm is exact, the small tolerance number is used only
39   * to account for round-off errors.
40   *
41   */
42  public final class FastFourierTransformerTest {
43      /** The common seed of all random number generators used in this test. */
44      private final static long SEED = 20110111L;
45  
46      /*
47       * Precondition checks.
48       */
49  
50      @Test
51      public void testTransformComplexSizeNotAPowerOfTwo() {
52          final int n = 127;
53          final Complex[] x = createComplexData(n);
54          final DftNormalization[] norm;
55          norm = DftNormalization.values();
56          final TransformType[] type;
57          type = TransformType.values();
58          for (int i = 0; i < norm.length; i++) {
59              for (int j = 0; j < type.length; j++) {
60                  final FastFourierTransformer fft;
61                  fft = new FastFourierTransformer(norm[i]);
62                  try {
63                      fft.transform(x, type[j]);
64                      Assert.fail(norm[i] + ", " + type[j] +
65                          ": MathIllegalArgumentException was expected");
66                  } catch (MathIllegalArgumentException e) {
67                      // Expected behaviour
68                  }
69              }
70          }
71      }
72  
73      @Test
74      public void testTransformRealSizeNotAPowerOfTwo() {
75          final int n = 127;
76          final double[] x = createRealData(n);
77          final DftNormalization[] norm;
78          norm = DftNormalization.values();
79          final TransformType[] type;
80          type = TransformType.values();
81          for (int i = 0; i < norm.length; i++) {
82              for (int j = 0; j < type.length; j++) {
83                  final FastFourierTransformer fft;
84                  fft = new FastFourierTransformer(norm[i]);
85                  try {
86                      fft.transform(x, type[j]);
87                      Assert.fail(norm[i] + ", " + type[j] +
88                          ": MathIllegalArgumentException was expected");
89                  } catch (MathIllegalArgumentException e) {
90                      // Expected behaviour
91                  }
92              }
93          }
94      }
95  
96      @Test
97      public void testTransformFunctionSizeNotAPowerOfTwo() {
98          final int n = 127;
99          final UnivariateFunction f = new Sin();
100         final DftNormalization[] norm;
101         norm = DftNormalization.values();
102         final TransformType[] type;
103         type = TransformType.values();
104         for (int i = 0; i < norm.length; i++) {
105             for (int j = 0; j < type.length; j++) {
106                 final FastFourierTransformer fft;
107                 fft = new FastFourierTransformer(norm[i]);
108                 try {
109                     fft.transform(f, 0.0, Math.PI, n, type[j]);
110                     Assert.fail(norm[i] + ", " + type[j] +
111                         ": MathIllegalArgumentException was expected");
112                 } catch (MathIllegalArgumentException e) {
113                     // Expected behaviour
114                 }
115             }
116         }
117     }
118 
119     @Test
120     public void testTransformFunctionNotStrictlyPositiveNumberOfSamples() {
121         final int n = -128;
122         final UnivariateFunction f = new Sin();
123         final DftNormalization[] norm;
124         norm = DftNormalization.values();
125         final TransformType[] type;
126         type = TransformType.values();
127         for (int i = 0; i < norm.length; i++) {
128             for (int j = 0; j < type.length; j++) {
129                 final FastFourierTransformer fft;
130                 fft = new FastFourierTransformer(norm[i]);
131                 try {
132                     fft.transform(f, 0.0, Math.PI, n, type[j]);
133                     fft.transform(f, 0.0, Math.PI, n, type[j]);
134                     Assert.fail(norm[i] + ", " + type[j] +
135                         ": MathIllegalArgumentException was expected");
136                 } catch (MathIllegalArgumentException e) {
137                     // Expected behaviour
138                 }
139             }
140         }
141     }
142 
143     @Test
144     public void testTransformFunctionInvalidBounds() {
145         final int n = 128;
146         final UnivariateFunction f = new Sin();
147         final DftNormalization[] norm;
148         norm = DftNormalization.values();
149         final TransformType[] type;
150         type = TransformType.values();
151         for (int i = 0; i < norm.length; i++) {
152             for (int j = 0; j < type.length; j++) {
153                 final FastFourierTransformer fft;
154                 fft = new FastFourierTransformer(norm[i]);
155                 try {
156                     fft.transform(f, Math.PI, 0.0, n, type[j]);
157                     Assert.fail(norm[i] + ", " + type[j] +
158                         ": MathIllegalArgumentException was expected");
159                 } catch (MathIllegalArgumentException e) {
160                     // Expected behaviour
161                 }
162             }
163         }
164     }
165 
166     @Test
167     public void testTransformOnePoint() {
168         double[][] data = new double[][] { { 1.0 }, { 2.0} };
169         FastFourierTransformer.transformInPlace(data,
170                                                 DftNormalization.STANDARD,
171                                                 TransformType.FORWARD);
172         Assert.assertEquals(1.0, data[0][0], 1.0e-15);
173         Assert.assertEquals(2.0, data[1][0], 1.0e-15);
174     }
175 
176     /*
177      * Utility methods for checking (successful) transforms.
178      */
179 
180     private static Complex[] createComplexData(final int n) {
181         final Random random = new Random(SEED);
182         final Complex[] data = new Complex[n];
183         for (int i = 0; i < n; i++) {
184             final double re = 2.0 * random.nextDouble() - 1.0;
185             final double im = 2.0 * random.nextDouble() - 1.0;
186             data[i] = new Complex(re, im);
187         }
188         return data;
189     }
190 
191     private static double[] createRealData(final int n) {
192         final Random random = new Random(SEED);
193         final double[] data = new double[n];
194         for (int i = 0; i < n; i++) {
195             data[i] = 2.0 * random.nextDouble() - 1.0;
196         }
197         return data;
198     }
199 
200     /** Naive implementation of DFT, for reference. */
201     private static Complex[] dft(final Complex[] x, final int sgn) {
202         final int n = x.length;
203         final double[] cos = new double[n];
204         final double[] sin = new double[n];
205         final Complex[] y = new Complex[n];
206         for (int i = 0; i < n; i++) {
207             final double arg = 2.0 * FastMath.PI * i / n;
208             cos[i] = FastMath.cos(arg);
209             sin[i] = FastMath.sin(arg);
210         }
211         for (int i = 0; i < n; i++) {
212             double yr = 0.0;
213             double yi = 0.0;
214             for (int j = 0; j < n; j++) {
215                 final int index = (i * j) % n;
216                 final double c = cos[index];
217                 final double s = sin[index];
218                 final double xr = x[j].getReal();
219                 final double xi = x[j].getImaginary();
220                 yr += c * xr - sgn * s * xi;
221                 yi += sgn * s * xr + c * xi;
222             }
223             y[i] = new Complex(yr, yi);
224         }
225         return y;
226     }
227 
228     private static void doTestTransformComplex(final int n, final double tol,
229         final DftNormalization normalization,
230         final TransformType type) {
231         final FastFourierTransformer fft;
232         fft = new FastFourierTransformer(normalization);
233         final Complex[] x = createComplexData(n);
234         final Complex[] expected;
235         final double s;
236         if (type==TransformType.FORWARD) {
237             expected = dft(x, -1);
238             if (normalization == DftNormalization.STANDARD){
239                 s = 1.0;
240             } else {
241                 s = 1.0 / FastMath.sqrt(n);
242             }
243         } else {
244             expected = dft(x, 1);
245             if (normalization == DftNormalization.STANDARD) {
246                 s = 1.0 / n;
247             } else {
248                 s = 1.0 / FastMath.sqrt(n);
249             }
250         }
251         final Complex[] actual = fft.transform(x, type);
252         for (int i = 0; i < n; i++) {
253             final String msg;
254             msg = String.format("%s, %s, %d, %d", normalization, type, n, i);
255             final double re = s * expected[i].getReal();
256             Assert.assertEquals(msg, re, actual[i].getReal(),
257                 tol * FastMath.abs(re));
258             final double im = s * expected[i].getImaginary();
259             Assert.assertEquals(msg, im, actual[i].getImaginary(), tol *
260                 FastMath.abs(re));
261         }
262     }
263 
264     private static void doTestTransformReal(final int n, final double tol,
265         final DftNormalization normalization,
266         final TransformType type) {
267         final FastFourierTransformer fft;
268         fft = new FastFourierTransformer(normalization);
269         final double[] x = createRealData(n);
270         final Complex[] xc = new Complex[n];
271         for (int i = 0; i < n; i++) {
272             xc[i] = new Complex(x[i], 0.0);
273         }
274         final Complex[] expected;
275         final double s;
276         if (type == TransformType.FORWARD) {
277             expected = dft(xc, -1);
278             if (normalization == DftNormalization.STANDARD) {
279                 s = 1.0;
280             } else {
281                 s = 1.0 / FastMath.sqrt(n);
282             }
283         } else {
284             expected = dft(xc, 1);
285             if (normalization == DftNormalization.STANDARD) {
286                 s = 1.0 / n;
287             } else {
288                 s = 1.0 / FastMath.sqrt(n);
289             }
290         }
291         final Complex[] actual = fft.transform(x, type);
292         for (int i = 0; i < n; i++) {
293             final String msg;
294             msg = String.format("%s, %s, %d, %d", normalization, type, n, i);
295             final double re = s * expected[i].getReal();
296             Assert.assertEquals(msg, re, actual[i].getReal(),
297                 tol * FastMath.abs(re));
298             final double im = s * expected[i].getImaginary();
299             Assert.assertEquals(msg, im, actual[i].getImaginary(), tol *
300                 FastMath.abs(re));
301         }
302     }
303 
304     private static void doTestTransformFunction(final UnivariateFunction f,
305         final double min, final double max, int n, final double tol,
306         final DftNormalization normalization,
307         final TransformType type) {
308         final FastFourierTransformer fft;
309         fft = new FastFourierTransformer(normalization);
310         final Complex[] x = new Complex[n];
311         for (int i = 0; i < n; i++) {
312             final double t = min + i * (max - min) / n;
313             x[i] = new Complex(f.value(t));
314         }
315         final Complex[] expected;
316         final double s;
317         if (type == TransformType.FORWARD) {
318             expected = dft(x, -1);
319             if (normalization == DftNormalization.STANDARD) {
320                 s = 1.0;
321             } else {
322                 s = 1.0 / FastMath.sqrt(n);
323             }
324         } else {
325             expected = dft(x, 1);
326             if (normalization == DftNormalization.STANDARD) {
327                 s = 1.0 / n;
328             } else {
329                 s = 1.0 / FastMath.sqrt(n);
330             }
331         }
332         final Complex[] actual = fft.transform(f, min, max, n, type);
333         for (int i = 0; i < n; i++) {
334             final String msg = String.format("%d, %d", n, i);
335             final double re = s * expected[i].getReal();
336             Assert.assertEquals(msg, re, actual[i].getReal(),
337                 tol * FastMath.abs(re));
338             final double im = s * expected[i].getImaginary();
339             Assert.assertEquals(msg, im, actual[i].getImaginary(), tol *
340                 FastMath.abs(re));
341         }
342     }
343 
344     /*
345      * Tests of standard transform (when data is valid).
346      */
347 
348     @Test
349     public void testTransformComplex() {
350         final DftNormalization[] norm;
351         norm = DftNormalization.values();
352         final TransformType[] type;
353         type = TransformType.values();
354         for (int i = 0; i < norm.length; i++) {
355             for (int j = 0; j < type.length; j++) {
356                 doTestTransformComplex(2, 1.0E-15, norm[i], type[j]);
357                 doTestTransformComplex(4, 1.0E-14, norm[i], type[j]);
358                 doTestTransformComplex(8, 1.0E-14, norm[i], type[j]);
359                 doTestTransformComplex(16, 1.0E-13, norm[i], type[j]);
360                 doTestTransformComplex(32, 1.0E-13, norm[i], type[j]);
361                 doTestTransformComplex(64, 1.0E-12, norm[i], type[j]);
362                 doTestTransformComplex(128, 1.0E-12, norm[i], type[j]);
363             }
364         }
365     }
366 
367     @Test
368     public void testStandardTransformReal() {
369         final DftNormalization[] norm;
370         norm = DftNormalization.values();
371         final TransformType[] type;
372         type = TransformType.values();
373         for (int i = 0; i < norm.length; i++) {
374             for (int j = 0; j < type.length; j++) {
375                 doTestTransformReal(2, 1.0E-15, norm[i], type[j]);
376                 doTestTransformReal(4, 1.0E-14, norm[i], type[j]);
377                 doTestTransformReal(8, 1.0E-14, norm[i], type[j]);
378                 doTestTransformReal(16, 1.0E-13, norm[i], type[j]);
379                 doTestTransformReal(32, 1.0E-13, norm[i], type[j]);
380                 doTestTransformReal(64, 1.0E-13, norm[i], type[j]);
381                 doTestTransformReal(128, 1.0E-11, norm[i], type[j]);
382             }
383         }
384     }
385 
386     @Test
387     public void testStandardTransformFunction() {
388         final UnivariateFunction f = new Sinc();
389         final double min = -FastMath.PI;
390         final double max = FastMath.PI;
391         final DftNormalization[] norm;
392         norm = DftNormalization.values();
393         final TransformType[] type;
394         type = TransformType.values();
395         for (int i = 0; i < norm.length; i++) {
396             for (int j = 0; j < type.length; j++) {
397                 doTestTransformFunction(f, min, max, 2, 1.0E-15, norm[i], type[j]);
398                 doTestTransformFunction(f, min, max, 4, 1.0E-14, norm[i], type[j]);
399                 doTestTransformFunction(f, min, max, 8, 1.0E-14, norm[i], type[j]);
400                 doTestTransformFunction(f, min, max, 16, 1.0E-13, norm[i], type[j]);
401                 doTestTransformFunction(f, min, max, 32, 1.0E-13, norm[i], type[j]);
402                 doTestTransformFunction(f, min, max, 64, 1.0E-12, norm[i], type[j]);
403                 doTestTransformFunction(f, min, max, 128, 1.0E-11, norm[i], type[j]);
404             }
405         }
406     }
407 
408     /*
409      * Additional tests for 1D data.
410      */
411 
412     /**
413      * Test of transformer for the ad hoc data taken from Mathematica.
414      */
415     @Test
416     public void testAdHocData() {
417         FastFourierTransformer transformer;
418         transformer = new FastFourierTransformer(DftNormalization.STANDARD);
419         Complex result[]; double tolerance = 1E-12;
420 
421         double x[] = {1.3, 2.4, 1.7, 4.1, 2.9, 1.7, 5.1, 2.7};
422         Complex y[] = {
423             new Complex(21.9, 0.0),
424             new Complex(-2.09497474683058, 1.91507575950825),
425             new Complex(-2.6, 2.7),
426             new Complex(-1.10502525316942, -4.88492424049175),
427             new Complex(0.1, 0.0),
428             new Complex(-1.10502525316942, 4.88492424049175),
429             new Complex(-2.6, -2.7),
430             new Complex(-2.09497474683058, -1.91507575950825)};
431 
432         result = transformer.transform(x, TransformType.FORWARD);
433         for (int i = 0; i < result.length; i++) {
434             Assert.assertEquals(y[i].getReal(), result[i].getReal(), tolerance);
435             Assert.assertEquals(y[i].getImaginary(), result[i].getImaginary(), tolerance);
436         }
437 
438         result = transformer.transform(y, TransformType.INVERSE);
439         for (int i = 0; i < result.length; i++) {
440             Assert.assertEquals(x[i], result[i].getReal(), tolerance);
441             Assert.assertEquals(0.0, result[i].getImaginary(), tolerance);
442         }
443 
444         double x2[] = {10.4, 21.6, 40.8, 13.6, 23.2, 32.8, 13.6, 19.2};
445         TransformUtils.scaleArray(x2, 1.0 / FastMath.sqrt(x2.length));
446         Complex y2[] = y;
447 
448         transformer = new FastFourierTransformer(DftNormalization.UNITARY);
449         result = transformer.transform(y2, TransformType.FORWARD);
450         for (int i = 0; i < result.length; i++) {
451             Assert.assertEquals(x2[i], result[i].getReal(), tolerance);
452             Assert.assertEquals(0.0, result[i].getImaginary(), tolerance);
453         }
454 
455         result = transformer.transform(x2, TransformType.INVERSE);
456         for (int i = 0; i < result.length; i++) {
457             Assert.assertEquals(y2[i].getReal(), result[i].getReal(), tolerance);
458             Assert.assertEquals(y2[i].getImaginary(), result[i].getImaginary(), tolerance);
459         }
460     }
461 
462     /**
463      * Test of transformer for the sine function.
464      */
465     @Test
466     public void testSinFunction() {
467         UnivariateFunction f = new Sin();
468         FastFourierTransformer transformer;
469         transformer = new FastFourierTransformer(DftNormalization.STANDARD);
470         Complex result[]; int N = 1 << 8;
471         double min, max, tolerance = 1E-12;
472 
473         min = 0.0; max = 2.0 * FastMath.PI;
474         result = transformer.transform(f, min, max, N, TransformType.FORWARD);
475         Assert.assertEquals(0.0, result[1].getReal(), tolerance);
476         Assert.assertEquals(-(N >> 1), result[1].getImaginary(), tolerance);
477         Assert.assertEquals(0.0, result[N-1].getReal(), tolerance);
478         Assert.assertEquals(N >> 1, result[N-1].getImaginary(), tolerance);
479         for (int i = 0; i < N-1; i += (i == 0 ? 2 : 1)) {
480             Assert.assertEquals(0.0, result[i].getReal(), tolerance);
481             Assert.assertEquals(0.0, result[i].getImaginary(), tolerance);
482         }
483 
484         min = -FastMath.PI; max = FastMath.PI;
485         result = transformer.transform(f, min, max, N, TransformType.INVERSE);
486         Assert.assertEquals(0.0, result[1].getReal(), tolerance);
487         Assert.assertEquals(-0.5, result[1].getImaginary(), tolerance);
488         Assert.assertEquals(0.0, result[N-1].getReal(), tolerance);
489         Assert.assertEquals(0.5, result[N-1].getImaginary(), tolerance);
490         for (int i = 0; i < N-1; i += (i == 0 ? 2 : 1)) {
491             Assert.assertEquals(0.0, result[i].getReal(), tolerance);
492             Assert.assertEquals(0.0, result[i].getImaginary(), tolerance);
493         }
494     }
495 
496 }