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