1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.transform;
18
19 import java.io.Serializable;
20 import java.lang.reflect.Array;
21
22 import org.apache.commons.math.FunctionEvaluationException;
23 import org.apache.commons.math.MathRuntimeException;
24 import org.apache.commons.math.analysis.UnivariateRealFunction;
25 import org.apache.commons.math.complex.Complex;
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46 public class FastFourierTransformer implements Serializable {
47
48
49 static final long serialVersionUID = 5138259215438106000L;
50
51
52 private RootsOfUnity roots = new RootsOfUnity();
53
54
55
56
57 public FastFourierTransformer() {
58 super();
59 }
60
61
62
63
64
65
66
67
68
69
70
71 public Complex[] transform(double f[])
72 throws IllegalArgumentException {
73 return fft(f, false);
74 }
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91 public Complex[] transform(UnivariateRealFunction f,
92 double min, double max, int n)
93 throws FunctionEvaluationException, IllegalArgumentException {
94 double data[] = sample(f, min, max, n);
95 return fft(data, false);
96 }
97
98
99
100
101
102
103
104
105
106
107
108 public Complex[] transform(Complex f[])
109 throws IllegalArgumentException {
110 roots.computeOmega(f.length);
111 return fft(f);
112 }
113
114
115
116
117
118
119
120
121
122
123
124 public Complex[] transform2(double f[])
125 throws IllegalArgumentException {
126
127 double scaling_coefficient = 1.0 / Math.sqrt(f.length);
128 return scaleArray(fft(f, false), scaling_coefficient);
129 }
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146 public Complex[] transform2(UnivariateRealFunction f,
147 double min, double max, int n)
148 throws FunctionEvaluationException, IllegalArgumentException {
149
150 double data[] = sample(f, min, max, n);
151 double scaling_coefficient = 1.0 / Math.sqrt(n);
152 return scaleArray(fft(data, false), scaling_coefficient);
153 }
154
155
156
157
158
159
160
161
162
163
164
165 public Complex[] transform2(Complex f[])
166 throws IllegalArgumentException {
167
168 roots.computeOmega(f.length);
169 double scaling_coefficient = 1.0 / Math.sqrt(f.length);
170 return scaleArray(fft(f), scaling_coefficient);
171 }
172
173
174
175
176
177
178
179
180
181
182
183 public Complex[] inversetransform(double f[])
184 throws IllegalArgumentException {
185
186 double scaling_coefficient = 1.0 / f.length;
187 return scaleArray(fft(f, true), scaling_coefficient);
188 }
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205 public Complex[] inversetransform(UnivariateRealFunction f,
206 double min, double max, int n)
207 throws FunctionEvaluationException, IllegalArgumentException {
208
209 double data[] = sample(f, min, max, n);
210 double scaling_coefficient = 1.0 / n;
211 return scaleArray(fft(data, true), scaling_coefficient);
212 }
213
214
215
216
217
218
219
220
221
222
223
224 public Complex[] inversetransform(Complex f[])
225 throws IllegalArgumentException {
226
227 roots.computeOmega(-f.length);
228 double scaling_coefficient = 1.0 / f.length;
229 return scaleArray(fft(f), scaling_coefficient);
230 }
231
232
233
234
235
236
237
238
239
240
241
242 public Complex[] inversetransform2(double f[])
243 throws IllegalArgumentException {
244
245 double scaling_coefficient = 1.0 / Math.sqrt(f.length);
246 return scaleArray(fft(f, true), scaling_coefficient);
247 }
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264 public Complex[] inversetransform2(UnivariateRealFunction f,
265 double min, double max, int n)
266 throws FunctionEvaluationException, IllegalArgumentException {
267
268 double data[] = sample(f, min, max, n);
269 double scaling_coefficient = 1.0 / Math.sqrt(n);
270 return scaleArray(fft(data, true), scaling_coefficient);
271 }
272
273
274
275
276
277
278
279
280
281
282
283 public Complex[] inversetransform2(Complex f[])
284 throws IllegalArgumentException {
285
286 roots.computeOmega(-f.length);
287 double scaling_coefficient = 1.0 / Math.sqrt(f.length);
288 return scaleArray(fft(f), scaling_coefficient);
289 }
290
291
292
293
294
295
296
297
298
299 protected Complex[] fft(double f[], boolean isInverse)
300 throws IllegalArgumentException {
301
302 verifyDataSet(f);
303 Complex F[] = new Complex[f.length];
304 if (f.length == 1) {
305 F[0] = new Complex(f[0], 0.0);
306 return F;
307 }
308
309
310
311 int N = f.length >> 1;
312 Complex c[] = new Complex[N];
313 for (int i = 0; i < N; i++) {
314 c[i] = new Complex(f[2*i], f[2*i+1]);
315 }
316 roots.computeOmega(isInverse ? -N : N);
317 Complex z[] = fft(c);
318
319
320 roots.computeOmega(isInverse ? -2*N : 2*N);
321 F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);
322 F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0);
323 for (int i = 1; i < N; i++) {
324 Complex A = z[N-i].conjugate();
325 Complex B = z[i].add(A);
326 Complex C = z[i].subtract(A);
327
328 Complex D = new Complex(-roots.getOmegaImaginary(i),
329 roots.getOmegaReal(i));
330 F[i] = B.subtract(C.multiply(D));
331 F[2*N-i] = F[i].conjugate();
332 }
333
334 return scaleArray(F, 0.5);
335 }
336
337
338
339
340
341
342
343
344 protected Complex[] fft(Complex data[])
345 throws IllegalArgumentException {
346
347 int i, j, k, m, N = data.length;
348 Complex A, B, C, D, E, F, z, f[] = new Complex[N];
349
350
351 verifyDataSet(data);
352 if (N == 1) {
353 f[0] = data[0];
354 return f;
355 }
356 if (N == 2) {
357 f[0] = data[0].add(data[1]);
358 f[1] = data[0].subtract(data[1]);
359 return f;
360 }
361
362
363 j = 0;
364 for (i = 0; i < N; i++) {
365 f[i] = data[j];
366 k = N >> 1;
367 while (j >= k && k > 0) {
368 j -= k; k >>= 1;
369 }
370 j += k;
371 }
372
373
374 for (i = 0; i < N; i += 4) {
375 A = f[i].add(f[i+1]);
376 B = f[i+2].add(f[i+3]);
377 C = f[i].subtract(f[i+1]);
378 D = f[i+2].subtract(f[i+3]);
379 E = C.add(D.multiply(Complex.I));
380 F = C.subtract(D.multiply(Complex.I));
381 f[i] = A.add(B);
382 f[i+2] = A.subtract(B);
383
384 f[i+1] = roots.isForward() ? F : E;
385 f[i+3] = roots.isForward() ? E : F;
386 }
387
388
389 for (i = 4; i < N; i <<= 1) {
390 m = N / (i<<1);
391 for (j = 0; j < N; j += i<<1) {
392 for (k = 0; k < i; k++) {
393
394 final int k_times_m = k*m;
395 final double omega_k_times_m_real = roots.getOmegaReal(k_times_m);
396 final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m);
397
398 z = new Complex(
399 f[i+j+k].getReal() * omega_k_times_m_real -
400 f[i+j+k].getImaginary() * omega_k_times_m_imaginary,
401 f[i+j+k].getReal() * omega_k_times_m_imaginary +
402 f[i+j+k].getImaginary() * omega_k_times_m_real);
403
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
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429 public static double[] sample(UnivariateRealFunction f,
430 double min, double max, int n)
431 throws FunctionEvaluationException, IllegalArgumentException {
432
433 if (n <= 0) {
434 throw MathRuntimeException.createIllegalArgumentException(
435 "number of sample is not positive: {0}",
436 n);
437 }
438 verifyInterval(min, max);
439
440 double s[] = new double[n];
441 double h = (max - min) / n;
442 for (int i = 0; i < n; i++) {
443 s[i] = f.value(min + i * h);
444 }
445 return s;
446 }
447
448
449
450
451
452
453
454
455
456 public static double[] scaleArray(double f[], double d) {
457 for (int i = 0; i < f.length; i++) {
458 f[i] *= d;
459 }
460 return f;
461 }
462
463
464
465
466
467
468
469
470
471 public static Complex[] scaleArray(Complex f[], double d) {
472 for (int i = 0; i < f.length; i++) {
473 f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary());
474 }
475 return f;
476 }
477
478
479
480
481
482
483
484 public static boolean isPowerOf2(long n) {
485 return (n > 0) && ((n & (n - 1)) == 0);
486 }
487
488
489
490
491
492
493
494 public static void verifyDataSet(double d[]) throws IllegalArgumentException {
495 if (!isPowerOf2(d.length)) {
496 throw MathRuntimeException.createIllegalArgumentException(
497 "{0} is not a power of 2, consider padding for fix",
498 d.length);
499 }
500 }
501
502
503
504
505
506
507
508 public static void verifyDataSet(Object o[]) throws IllegalArgumentException {
509 if (!isPowerOf2(o.length)) {
510 throw MathRuntimeException.createIllegalArgumentException(
511 "{0} is not a power of 2, consider padding for fix",
512 o.length);
513 }
514 }
515
516
517
518
519
520
521
522
523 public static void verifyInterval(double lower, double upper)
524 throws IllegalArgumentException {
525
526 if (lower >= upper) {
527 throw MathRuntimeException.createIllegalArgumentException(
528 "endpoints do not specify an interval: [{0}, {1}]",
529 lower, upper);
530 }
531 }
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547 public Object mdfft(Object mdca, boolean forward)
548 throws IllegalArgumentException {
549 MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)
550 new MultiDimensionalComplexMatrix(mdca).clone();
551 int[] dimensionSize = mdcm.getDimensionSizes();
552
553 for (int i = 0; i < dimensionSize.length; i++) {
554 mdfft(mdcm, forward, i, new int[0]);
555 }
556 return mdcm.getArray();
557 }
558
559
560
561
562
563
564
565
566
567
568 private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward,
569 int d, int[] subVector)
570 throws IllegalArgumentException {
571 int[] dimensionSize = mdcm.getDimensionSizes();
572
573 if (subVector.length == dimensionSize.length) {
574 Complex[] temp = new Complex[dimensionSize[d]];
575 for (int i = 0; i < dimensionSize[d]; i++) {
576
577 subVector[d] = i;
578 temp[i] = mdcm.get(subVector);
579 }
580
581 if (forward)
582 temp = transform2(temp);
583 else
584 temp = inversetransform2(temp);
585
586 for (int i = 0; i < dimensionSize[d]; i++) {
587 subVector[d] = i;
588 mdcm.set(temp[i], subVector);
589 }
590 } else {
591 int[] vector = new int[subVector.length + 1];
592 System.arraycopy(subVector, 0, vector, 0, subVector.length);
593 if (subVector.length == d) {
594
595
596 vector[d] = 0;
597 mdfft(mdcm, forward, d, vector);
598 } else {
599 for (int i = 0; i < dimensionSize[subVector.length]; i++) {
600 vector[subVector.length] = i;
601
602 mdfft(mdcm, forward, d, vector);
603 }
604 }
605 }
606 return;
607 }
608
609
610
611
612
613
614
615
616 private static class MultiDimensionalComplexMatrix
617 implements Cloneable {
618
619
620 protected int[] dimensionSize;
621
622
623 protected Object multiDimensionalComplexArray;
624
625
626
627
628 public MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) {
629
630 this.multiDimensionalComplexArray = multiDimensionalComplexArray;
631
632
633 int numOfDimensions = 0;
634 for (Object lastDimension = multiDimensionalComplexArray;
635 lastDimension instanceof Object[];) {
636 final Object[] array = (Object[]) lastDimension;
637 numOfDimensions++;
638 lastDimension = array[0];
639 }
640
641
642 dimensionSize = new int[numOfDimensions];
643
644
645 numOfDimensions = 0;
646 for (Object lastDimension = multiDimensionalComplexArray;
647 lastDimension instanceof Object[];) {
648 final Object[] array = (Object[]) lastDimension;
649 dimensionSize[numOfDimensions++] = array.length;
650 lastDimension = array[0];
651 }
652
653 }
654
655
656
657
658
659
660
661 public Complex get(int... vector)
662 throws IllegalArgumentException {
663 if (vector == null) {
664 if (dimensionSize.length > 0) {
665 throw MathRuntimeException.createIllegalArgumentException(
666 "some dimensions don't match: {0} != {1}",
667 0, dimensionSize.length);
668 }
669 return null;
670 }
671 if (vector.length != dimensionSize.length) {
672 throw MathRuntimeException.createIllegalArgumentException(
673 "some dimensions don't match: {0} != {1}",
674 vector.length, dimensionSize.length);
675 }
676
677 Object lastDimension = multiDimensionalComplexArray;
678
679 for (int i = 0; i < dimensionSize.length; i++) {
680 lastDimension = ((Object[]) lastDimension)[vector[i]];
681 }
682 return (Complex) lastDimension;
683 }
684
685
686
687
688
689
690
691
692 public Complex set(Complex magnitude, int... vector)
693 throws IllegalArgumentException {
694 if (vector == null) {
695 if (dimensionSize.length > 0) {
696 throw MathRuntimeException.createIllegalArgumentException(
697 "some dimensions don't match: {0} != {1}",
698 0, dimensionSize.length);
699 }
700 return null;
701 }
702 if (vector.length != dimensionSize.length) {
703 throw MathRuntimeException.createIllegalArgumentException(
704 "some dimensions don't match: {0} != {1}",
705 vector.length,dimensionSize.length);
706 }
707
708 Object[] lastDimension = (Object[]) multiDimensionalComplexArray;
709 for (int i = 0; i < dimensionSize.length - 1; i++) {
710 lastDimension = (Object[]) lastDimension[vector[i]];
711 }
712
713 Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]];
714 lastDimension[vector[dimensionSize.length - 1]] = magnitude;
715
716 return lastValue;
717 }
718
719
720
721
722
723 public int[] getDimensionSizes() {
724 return dimensionSize.clone();
725 }
726
727
728
729
730
731 public Object getArray() {
732 return multiDimensionalComplexArray;
733 }
734
735
736 @Override
737 public Object clone() {
738 MultiDimensionalComplexMatrix mdcm =
739 new MultiDimensionalComplexMatrix(Array.newInstance(
740 Complex.class, dimensionSize));
741 clone(mdcm);
742 return mdcm;
743 }
744
745
746
747
748
749 private void clone(MultiDimensionalComplexMatrix mdcm) {
750 int[] vector = new int[dimensionSize.length];
751 int size = 1;
752 for (int i = 0; i < dimensionSize.length; i++) {
753 size *= dimensionSize[i];
754 }
755 int[][] vectorList = new int[size][dimensionSize.length];
756 for (int[] nextVector: vectorList) {
757 System.arraycopy(vector, 0, nextVector, 0,
758 dimensionSize.length);
759 for (int i = 0; i < dimensionSize.length; i++) {
760 vector[i]++;
761 if (vector[i] < dimensionSize[i]) {
762 break;
763 } else {
764 vector[i] = 0;
765 }
766 }
767 }
768
769 for (int[] nextVector: vectorList) {
770 mdcm.set(get(nextVector), nextVector);
771 }
772 }
773 }
774
775
776
777
778
779 private static class RootsOfUnity implements Serializable {
780
781
782 private static final long serialVersionUID = 6404784357747329667L;
783
784
785 private int omegaCount;
786
787
788 private double[] omegaReal;
789
790
791 private double[] omegaImaginaryForward;
792
793
794 private double[] omegaImaginaryInverse;
795
796
797 private boolean isForward;
798
799
800
801
802 public RootsOfUnity() {
803
804 omegaCount = 0;
805 omegaReal = null;
806 omegaImaginaryForward = null;
807 omegaImaginaryInverse = null;
808 isForward = true;
809
810 }
811
812
813
814
815
816
817 public synchronized boolean isForward() throws IllegalStateException {
818
819 if (omegaCount == 0) {
820 throw MathRuntimeException.createIllegalStateException(
821 "roots of unity have not been computed yet");
822 }
823 return isForward;
824
825 }
826
827
828
829
830
831
832
833
834
835
836 public synchronized void computeOmega(int n) throws IllegalArgumentException {
837
838 if (n == 0) {
839 throw MathRuntimeException.createIllegalArgumentException(
840 "cannot compute 0-th root of unity, indefinite result");
841 }
842
843 isForward = (n > 0);
844
845
846 final int absN = Math.abs(n);
847
848 if (absN == omegaCount) {
849 return;
850 }
851
852
853 final double t = 2.0 * Math.PI / absN;
854 final double cosT = Math.cos(t);
855 final double sinT = Math.sin(t);
856 omegaReal = new double[absN];
857 omegaImaginaryForward = new double[absN];
858 omegaImaginaryInverse = new double[absN];
859 omegaReal[0] = 1.0;
860 omegaImaginaryForward[0] = 0.0;
861 omegaImaginaryInverse[0] = 0.0;
862 for (int i = 1; i < absN; i++) {
863 omegaReal[i] =
864 omegaReal[i-1] * cosT + omegaImaginaryForward[i-1] * sinT;
865 omegaImaginaryForward[i] =
866 omegaImaginaryForward[i-1] * cosT - omegaReal[i-1] * sinT;
867 omegaImaginaryInverse[i] = -omegaImaginaryForward[i];
868 }
869 omegaCount = absN;
870
871 }
872
873
874
875
876
877
878
879
880 public synchronized double getOmegaReal(int k)
881 throws IllegalStateException, IllegalArgumentException {
882
883 if (omegaCount == 0) {
884 throw MathRuntimeException.createIllegalStateException(
885 "roots of unity have not been computed yet");
886 }
887 if ((k < 0) || (k >= omegaCount)) {
888 throw MathRuntimeException.createIllegalArgumentException(
889 "out of range root of unity index {0} (must be in [{1};{2}])",
890 k, 0, omegaCount - 1);
891 }
892
893 return omegaReal[k];
894
895 }
896
897
898
899
900
901
902
903
904 public synchronized double getOmegaImaginary(int k)
905 throws IllegalStateException, IllegalArgumentException {
906
907 if (omegaCount == 0) {
908 throw MathRuntimeException.createIllegalStateException(
909 "roots of unity have not been computed yet");
910 }
911 if ((k < 0) || (k >= omegaCount)) {
912 throw MathRuntimeException.createIllegalArgumentException(
913 "out of range root of unity index {0} (must be in [{1};{2}])",
914 k, 0, omegaCount - 1);
915 }
916
917 return (isForward) ?
918 omegaImaginaryForward[k] : omegaImaginaryInverse[k];
919
920 }
921
922 }
923
924 }