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 java.lang.reflect.Array;
21
22 import org.apache.commons.math.MathRuntimeException;
23 import org.apache.commons.math.analysis.UnivariateRealFunction;
24 import org.apache.commons.math.complex.Complex;
25 import org.apache.commons.math.exception.util.LocalizedFormats;
26 import org.apache.commons.math.util.FastMath;
27
28 /**
29 * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html">
30 * Fast Fourier Transform</a> for transformation of one-dimensional data sets.
31 * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897,
32 * chapter 6.
33 * <p>
34 * There are several conventions for the definition of FFT and inverse FFT,
35 * mainly on different coefficient and exponent. Here the equations are listed
36 * in the comments of the corresponding methods.</p>
37 * <p>
38 * We require the length of data set to be power of 2, this greatly simplifies
39 * and speeds up the code. Users can pad the data with zeros to meet this
40 * requirement. There are other flavors of FFT, for reference, see S. Winograd,
41 * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation,
42 * 32 (1978), 175 - 199.</p>
43 *
44 * @version $Id: FastFourierTransformer.java 1165808 2011-09-06 19:59:47Z luc $
45 * @since 1.2
46 */
47 public class FastFourierTransformer implements Serializable {
48
49 /** Serializable version identifier. */
50 static final long serialVersionUID = 5138259215438106000L;
51
52 /** roots of unity */
53 private RootsOfUnity roots = new RootsOfUnity();
54
55 /**
56 * Construct a default transformer.
57 */
58 public FastFourierTransformer() {
59 super();
60 }
61
62 /**
63 * Transform the given real data set.
64 * <p>
65 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
66 * </p>
67 *
68 * @param f the real data array to be transformed
69 * @return the complex transformed array
70 * @throws IllegalArgumentException if any parameters are invalid
71 */
72 public Complex[] transform(double f[])
73 throws IllegalArgumentException {
74 return fft(f, false);
75 }
76
77 /**
78 * Transform the given real function, sampled on the given interval.
79 * <p>
80 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
81 * </p>
82 *
83 * @param f the function to be sampled and transformed
84 * @param min the lower bound for the interval
85 * @param max the upper bound for the interval
86 * @param n the number of sample points
87 * @return the complex transformed array
88 * @throws IllegalArgumentException if any parameters are invalid
89 */
90 public Complex[] transform(UnivariateRealFunction f,
91 double min, double max, int n)
92 throws IllegalArgumentException {
93 double data[] = sample(f, min, max, n);
94 return fft(data, false);
95 }
96
97 /**
98 * Transform the given complex data set.
99 * <p>
100 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
101 * </p>
102 *
103 * @param f the complex data array to be transformed
104 * @return the complex transformed array
105 * @throws IllegalArgumentException if any parameters are invalid
106 */
107 public Complex[] transform(Complex f[])
108 throws IllegalArgumentException {
109 roots.computeOmega(f.length);
110 return fft(f);
111 }
112
113 /**
114 * Transform the given real data set.
115 * <p>
116 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
117 * </p>
118 *
119 * @param f the real data array to be transformed
120 * @return the complex transformed array
121 * @throws IllegalArgumentException if any parameters are invalid
122 */
123 public Complex[] transform2(double f[])
124 throws IllegalArgumentException {
125
126 double scaling_coefficient = 1.0 / FastMath.sqrt(f.length);
127 return scaleArray(fft(f, false), scaling_coefficient);
128 }
129
130 /**
131 * Transform the given real function, sampled on the given interval.
132 * <p>
133 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
134 * </p>
135 *
136 * @param f the function to be sampled and transformed
137 * @param min the lower bound for the interval
138 * @param max the upper bound for the interval
139 * @param n the number of sample points
140 * @return the complex transformed array
141 * @throws IllegalArgumentException if any parameters are invalid
142 */
143 public Complex[] transform2(UnivariateRealFunction f,
144 double min, double max, int n)
145 throws IllegalArgumentException {
146
147 double data[] = sample(f, min, max, n);
148 double scaling_coefficient = 1.0 / FastMath.sqrt(n);
149 return scaleArray(fft(data, false), scaling_coefficient);
150 }
151
152 /**
153 * Transform the given complex data set.
154 * <p>
155 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
156 * </p>
157 *
158 * @param f the complex data array to be transformed
159 * @return the complex transformed array
160 * @throws IllegalArgumentException if any parameters are invalid
161 */
162 public Complex[] transform2(Complex f[])
163 throws IllegalArgumentException {
164
165 roots.computeOmega(f.length);
166 double scaling_coefficient = 1.0 / FastMath.sqrt(f.length);
167 return scaleArray(fft(f), scaling_coefficient);
168 }
169
170 /**
171 * Inversely transform the given real data set.
172 * <p>
173 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
174 * </p>
175 *
176 * @param f the real data array to be inversely transformed
177 * @return the complex inversely transformed array
178 * @throws IllegalArgumentException if any parameters are invalid
179 */
180 public Complex[] inversetransform(double f[])
181 throws IllegalArgumentException {
182
183 double scaling_coefficient = 1.0 / f.length;
184 return scaleArray(fft(f, true), scaling_coefficient);
185 }
186
187 /**
188 * Inversely transform the given real function, sampled on the given interval.
189 * <p>
190 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
191 * </p>
192 *
193 * @param f the function to be sampled and inversely transformed
194 * @param min the lower bound for the interval
195 * @param max the upper bound for the interval
196 * @param n the number of sample points
197 * @return the complex inversely transformed array
198 * @throws IllegalArgumentException if any parameters are invalid
199 */
200 public Complex[] inversetransform(UnivariateRealFunction f,
201 double min, double max, int n)
202 throws IllegalArgumentException {
203
204 double data[] = sample(f, min, max, n);
205 double scaling_coefficient = 1.0 / n;
206 return scaleArray(fft(data, true), scaling_coefficient);
207 }
208
209 /**
210 * Inversely transform the given complex data set.
211 * <p>
212 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
213 * </p>
214 *
215 * @param f the complex data array to be inversely transformed
216 * @return the complex inversely transformed array
217 * @throws IllegalArgumentException if any parameters are invalid
218 */
219 public Complex[] inversetransform(Complex f[])
220 throws IllegalArgumentException {
221
222 roots.computeOmega(-f.length); // pass negative argument
223 double scaling_coefficient = 1.0 / f.length;
224 return scaleArray(fft(f), scaling_coefficient);
225 }
226
227 /**
228 * Inversely transform the given real data set.
229 * <p>
230 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
231 * </p>
232 *
233 * @param f the real data array to be inversely transformed
234 * @return the complex inversely transformed array
235 * @throws IllegalArgumentException if any parameters are invalid
236 */
237 public Complex[] inversetransform2(double f[])
238 throws IllegalArgumentException {
239
240 double scaling_coefficient = 1.0 / FastMath.sqrt(f.length);
241 return scaleArray(fft(f, true), scaling_coefficient);
242 }
243
244 /**
245 * Inversely transform the given real function, sampled on the given interval.
246 * <p>
247 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
248 * </p>
249 *
250 * @param f the function to be sampled and inversely transformed
251 * @param min the lower bound for the interval
252 * @param max the upper bound for the interval
253 * @param n the number of sample points
254 * @return the complex inversely transformed array
255 * @throws IllegalArgumentException if any parameters are invalid
256 */
257 public Complex[] inversetransform2(UnivariateRealFunction f,
258 double min, double max, int n)
259 throws IllegalArgumentException {
260
261 double data[] = sample(f, min, max, n);
262 double scaling_coefficient = 1.0 / FastMath.sqrt(n);
263 return scaleArray(fft(data, true), scaling_coefficient);
264 }
265
266 /**
267 * Inversely transform the given complex data set.
268 * <p>
269 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
270 * </p>
271 *
272 * @param f the complex data array to be inversely transformed
273 * @return the complex inversely transformed array
274 * @throws IllegalArgumentException if any parameters are invalid
275 */
276 public Complex[] inversetransform2(Complex f[])
277 throws IllegalArgumentException {
278
279 roots.computeOmega(-f.length); // pass negative argument
280 double scaling_coefficient = 1.0 / FastMath.sqrt(f.length);
281 return scaleArray(fft(f), scaling_coefficient);
282 }
283
284 /**
285 * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
286 *
287 * @param f the real data array to be transformed
288 * @param isInverse the indicator of forward or inverse transform
289 * @return the complex transformed array
290 * @throws IllegalArgumentException if any parameters are invalid
291 */
292 protected Complex[] fft(double f[], boolean isInverse)
293 throws IllegalArgumentException {
294
295 verifyDataSet(f);
296 Complex F[] = new Complex[f.length];
297 if (f.length == 1) {
298 F[0] = new Complex(f[0], 0.0);
299 return F;
300 }
301
302 // Rather than the naive real to complex conversion, pack 2N
303 // real numbers into N complex numbers for better performance.
304 int N = f.length >> 1;
305 Complex c[] = new Complex[N];
306 for (int i = 0; i < N; i++) {
307 c[i] = new Complex(f[2*i], f[2*i+1]);
308 }
309 roots.computeOmega(isInverse ? -N : N);
310 Complex z[] = fft(c);
311
312 // reconstruct the FFT result for the original array
313 roots.computeOmega(isInverse ? -2*N : 2*N);
314 F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);
315 F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0);
316 for (int i = 1; i < N; i++) {
317 Complex A = z[N-i].conjugate();
318 Complex B = z[i].add(A);
319 Complex C = z[i].subtract(A);
320 //Complex D = roots.getOmega(i).multiply(Complex.I);
321 Complex D = new Complex(-roots.getOmegaImaginary(i),
322 roots.getOmegaReal(i));
323 F[i] = B.subtract(C.multiply(D));
324 F[2*N-i] = F[i].conjugate();
325 }
326
327 return scaleArray(F, 0.5);
328 }
329
330 /**
331 * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
332 *
333 * @param data the complex data array to be transformed
334 * @return the complex transformed array
335 * @throws IllegalArgumentException if any parameters are invalid
336 */
337 protected Complex[] fft(Complex data[])
338 throws IllegalArgumentException {
339
340 final int n = data.length;
341 final Complex f[] = new Complex[n];
342
343 // initial simple cases
344 verifyDataSet(data);
345 if (n == 1) {
346 f[0] = data[0];
347 return f;
348 }
349 if (n == 2) {
350 f[0] = data[0].add(data[1]);
351 f[1] = data[0].subtract(data[1]);
352 return f;
353 }
354
355 // permute original data array in bit-reversal order
356 int ii = 0;
357 for (int i = 0; i < n; i++) {
358 f[i] = data[ii];
359 int k = n >> 1;
360 while (ii >= k && k > 0) {
361 ii -= k; k >>= 1;
362 }
363 ii += k;
364 }
365
366 // the bottom base-4 round
367 for (int i = 0; i < n; i += 4) {
368 final Complex a = f[i].add(f[i+1]);
369 final Complex b = f[i+2].add(f[i+3]);
370 final Complex c = f[i].subtract(f[i+1]);
371 final Complex d = f[i+2].subtract(f[i+3]);
372 final Complex e1 = c.add(d.multiply(Complex.I));
373 final Complex e2 = c.subtract(d.multiply(Complex.I));
374 f[i] = a.add(b);
375 f[i+2] = a.subtract(b);
376 // omegaCount indicates forward or inverse transform
377 f[i+1] = roots.isForward() ? e2 : e1;
378 f[i+3] = roots.isForward() ? e1 : e2;
379 }
380
381 // iterations from bottom to top take O(N*logN) time
382 for (int i = 4; i < n; i <<= 1) {
383 final int m = n / (i<<1);
384 for (int j = 0; j < n; j += i<<1) {
385 for (int k = 0; k < i; k++) {
386 //z = f[i+j+k].multiply(roots.getOmega(k*m));
387 final int k_times_m = k*m;
388 final double omega_k_times_m_real = roots.getOmegaReal(k_times_m);
389 final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m);
390 //z = f[i+j+k].multiply(omega[k*m]);
391 final Complex z = new Complex(
392 f[i+j+k].getReal() * omega_k_times_m_real -
393 f[i+j+k].getImaginary() * omega_k_times_m_imaginary,
394 f[i+j+k].getReal() * omega_k_times_m_imaginary +
395 f[i+j+k].getImaginary() * omega_k_times_m_real);
396
397 f[i+j+k] = f[j+k].subtract(z);
398 f[j+k] = f[j+k].add(z);
399 }
400 }
401 }
402 return f;
403 }
404
405 /**
406 * Sample the given univariate real function on the given interval.
407 * <p>
408 * The interval is divided equally into N sections and sample points
409 * are taken from min to max-(max-min)/N. Usually f(x) is periodic
410 * such that f(min) = f(max) (note max is not sampled), but we don't
411 * require that.</p>
412 *
413 * @param f the function to be sampled
414 * @param min the lower bound for the interval
415 * @param max the upper bound for the interval
416 * @param n the number of sample points
417 * @return the samples array
418 * @throws IllegalArgumentException if any parameters are invalid
419 */
420 public static double[] sample(UnivariateRealFunction f, double min, double max, int n)
421 throws IllegalArgumentException {
422
423 if (n <= 0) {
424 throw MathRuntimeException.createIllegalArgumentException(
425 LocalizedFormats.NOT_POSITIVE_NUMBER_OF_SAMPLES,
426 n);
427 }
428 verifyInterval(min, max);
429
430 double s[] = new double[n];
431 double h = (max - min) / n;
432 for (int i = 0; i < n; i++) {
433 s[i] = f.value(min + i * h);
434 }
435 return s;
436 }
437
438 /**
439 * Multiply every component in the given real array by the
440 * given real number. The change is made in place.
441 *
442 * @param f the real array to be scaled
443 * @param d the real scaling coefficient
444 * @return a reference to the scaled array
445 */
446 public static double[] scaleArray(double f[], double d) {
447 for (int i = 0; i < f.length; i++) {
448 f[i] *= d;
449 }
450 return f;
451 }
452
453 /**
454 * Multiply every component in the given complex array by the
455 * given real number. The change is made in place.
456 *
457 * @param f the complex array to be scaled
458 * @param d the real scaling coefficient
459 * @return a reference to the scaled array
460 */
461 public static Complex[] scaleArray(Complex f[], double d) {
462 for (int i = 0; i < f.length; i++) {
463 f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary());
464 }
465 return f;
466 }
467
468 /**
469 * Returns true if the argument is power of 2.
470 *
471 * @param n the number to test
472 * @return true if the argument is power of 2
473 */
474 public static boolean isPowerOf2(long n) {
475 return (n > 0) && ((n & (n - 1)) == 0);
476 }
477
478 /**
479 * Verifies that the data set has length of power of 2.
480 *
481 * @param d the data array
482 * @throws IllegalArgumentException if array length is not power of 2
483 */
484 public static void verifyDataSet(double d[]) throws IllegalArgumentException {
485 if (!isPowerOf2(d.length)) {
486 throw MathRuntimeException.createIllegalArgumentException(
487 LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, d.length);
488 }
489 }
490
491 /**
492 * Verifies that the data set has length of power of 2.
493 *
494 * @param o the data array
495 * @throws IllegalArgumentException if array length is not power of 2
496 */
497 public static void verifyDataSet(Object o[]) throws IllegalArgumentException {
498 if (!isPowerOf2(o.length)) {
499 throw MathRuntimeException.createIllegalArgumentException(
500 LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, o.length);
501 }
502 }
503
504 /**
505 * Verifies that the endpoints specify an interval.
506 *
507 * @param lower lower endpoint
508 * @param upper upper endpoint
509 * @throws IllegalArgumentException if not interval
510 */
511 public static void verifyInterval(double lower, double upper)
512 throws IllegalArgumentException {
513
514 if (lower >= upper) {
515 throw MathRuntimeException.createIllegalArgumentException(
516 LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
517 lower, upper);
518 }
519 }
520
521 /**
522 * Performs a multi-dimensional Fourier transform on a given array.
523 * Use {@link #inversetransform2(Complex[])} and
524 * {@link #transform2(Complex[])} in a row-column implementation
525 * in any number of dimensions with O(N×log(N)) complexity with
526 * N=n<sub>1</sub>×n<sub>2</sub>×n<sub>3</sub>×...×n<sub>d</sub>,
527 * n<sub>x</sub>=number of elements in dimension x,
528 * and d=total number of dimensions.
529 *
530 * @param mdca Multi-Dimensional Complex Array id est Complex[][][][]
531 * @param forward inverseTransform2 is preformed if this is false
532 * @return transform of mdca as a Multi-Dimensional Complex Array id est Complex[][][][]
533 * @throws IllegalArgumentException if any dimension is not a power of two
534 */
535 public Object mdfft(Object mdca, boolean forward)
536 throws IllegalArgumentException {
537 MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)
538 new MultiDimensionalComplexMatrix(mdca).clone();
539 int[] dimensionSize = mdcm.getDimensionSizes();
540 //cycle through each dimension
541 for (int i = 0; i < dimensionSize.length; i++) {
542 mdfft(mdcm, forward, i, new int[0]);
543 }
544 return mdcm.getArray();
545 }
546
547 /**
548 * Performs one dimension of a multi-dimensional Fourier transform.
549 *
550 * @param mdcm input matrix
551 * @param forward inverseTransform2 is preformed if this is false
552 * @param d index of the dimension to process
553 * @param subVector recursion subvector
554 * @throws IllegalArgumentException if any dimension is not a power of two
555 */
556 private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward,
557 int d, int[] subVector)
558 throws IllegalArgumentException {
559 int[] dimensionSize = mdcm.getDimensionSizes();
560 //if done
561 if (subVector.length == dimensionSize.length) {
562 Complex[] temp = new Complex[dimensionSize[d]];
563 for (int i = 0; i < dimensionSize[d]; i++) {
564 //fft along dimension d
565 subVector[d] = i;
566 temp[i] = mdcm.get(subVector);
567 }
568
569 if (forward) {
570 temp = transform2(temp);
571 } else {
572 temp = inversetransform2(temp);
573 }
574
575 for (int i = 0; i < dimensionSize[d]; i++) {
576 subVector[d] = i;
577 mdcm.set(temp[i], subVector);
578 }
579 } else {
580 int[] vector = new int[subVector.length + 1];
581 System.arraycopy(subVector, 0, vector, 0, subVector.length);
582 if (subVector.length == d) {
583 //value is not important once the recursion is done.
584 //then an fft will be applied along the dimension d.
585 vector[d] = 0;
586 mdfft(mdcm, forward, d, vector);
587 } else {
588 for (int i = 0; i < dimensionSize[subVector.length]; i++) {
589 vector[subVector.length] = i;
590 //further split along the next dimension
591 mdfft(mdcm, forward, d, vector);
592 }
593 }
594 }
595 return;
596 }
597
598 /**
599 * Complex matrix implementation.
600 * Not designed for synchronized access
601 * may eventually be replaced by jsr-83 of the java community process
602 * http://jcp.org/en/jsr/detail?id=83
603 * may require additional exception throws for other basic requirements.
604 */
605 private static class MultiDimensionalComplexMatrix
606 implements Cloneable {
607
608 /** Size in all dimensions. */
609 protected int[] dimensionSize;
610
611 /** Storage array. */
612 protected Object multiDimensionalComplexArray;
613
614 /** Simple constructor.
615 * @param multiDimensionalComplexArray array containing the matrix elements
616 */
617 public MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) {
618
619 this.multiDimensionalComplexArray = multiDimensionalComplexArray;
620
621 // count dimensions
622 int numOfDimensions = 0;
623 for (Object lastDimension = multiDimensionalComplexArray;
624 lastDimension instanceof Object[];) {
625 final Object[] array = (Object[]) lastDimension;
626 numOfDimensions++;
627 lastDimension = array[0];
628 }
629
630 // allocate array with exact count
631 dimensionSize = new int[numOfDimensions];
632
633 // fill array
634 numOfDimensions = 0;
635 for (Object lastDimension = multiDimensionalComplexArray;
636 lastDimension instanceof Object[];) {
637 final Object[] array = (Object[]) lastDimension;
638 dimensionSize[numOfDimensions++] = array.length;
639 lastDimension = array[0];
640 }
641
642 }
643
644 /**
645 * Get a matrix element.
646 * @param vector indices of the element
647 * @return matrix element
648 * @exception IllegalArgumentException if dimensions do not match
649 */
650 public Complex get(int... vector)
651 throws IllegalArgumentException {
652 if (vector == null) {
653 if (dimensionSize.length > 0) {
654 throw MathRuntimeException.createIllegalArgumentException(
655 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 0, dimensionSize.length);
656 }
657 return null;
658 }
659 if (vector.length != dimensionSize.length) {
660 throw MathRuntimeException.createIllegalArgumentException(
661 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, vector.length, dimensionSize.length);
662 }
663
664 Object lastDimension = multiDimensionalComplexArray;
665
666 for (int i = 0; i < dimensionSize.length; i++) {
667 lastDimension = ((Object[]) lastDimension)[vector[i]];
668 }
669 return (Complex) lastDimension;
670 }
671
672 /**
673 * Set a matrix element.
674 * @param magnitude magnitude of the element
675 * @param vector indices of the element
676 * @return the previous value
677 * @exception IllegalArgumentException if dimensions do not match
678 */
679 public Complex set(Complex magnitude, int... vector)
680 throws IllegalArgumentException {
681 if (vector == null) {
682 if (dimensionSize.length > 0) {
683 throw MathRuntimeException.createIllegalArgumentException(
684 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 0, dimensionSize.length);
685 }
686 return null;
687 }
688 if (vector.length != dimensionSize.length) {
689 throw MathRuntimeException.createIllegalArgumentException(
690 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, vector.length,dimensionSize.length);
691 }
692
693 Object[] lastDimension = (Object[]) multiDimensionalComplexArray;
694 for (int i = 0; i < dimensionSize.length - 1; i++) {
695 lastDimension = (Object[]) lastDimension[vector[i]];
696 }
697
698 Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]];
699 lastDimension[vector[dimensionSize.length - 1]] = magnitude;
700
701 return lastValue;
702 }
703
704 /**
705 * Get the size in all dimensions.
706 * @return size in all dimensions
707 */
708 public int[] getDimensionSizes() {
709 return dimensionSize.clone();
710 }
711
712 /**
713 * Get the underlying storage array
714 * @return underlying storage array
715 */
716 public Object getArray() {
717 return multiDimensionalComplexArray;
718 }
719
720 /** {@inheritDoc} */
721 @Override
722 public Object clone() {
723 MultiDimensionalComplexMatrix mdcm =
724 new MultiDimensionalComplexMatrix(Array.newInstance(
725 Complex.class, dimensionSize));
726 clone(mdcm);
727 return mdcm;
728 }
729
730 /**
731 * Copy contents of current array into mdcm.
732 * @param mdcm array where to copy data
733 */
734 private void clone(MultiDimensionalComplexMatrix mdcm) {
735 int[] vector = new int[dimensionSize.length];
736 int size = 1;
737 for (int i = 0; i < dimensionSize.length; i++) {
738 size *= dimensionSize[i];
739 }
740 int[][] vectorList = new int[size][dimensionSize.length];
741 for (int[] nextVector: vectorList) {
742 System.arraycopy(vector, 0, nextVector, 0,
743 dimensionSize.length);
744 for (int i = 0; i < dimensionSize.length; i++) {
745 vector[i]++;
746 if (vector[i] < dimensionSize[i]) {
747 break;
748 } else {
749 vector[i] = 0;
750 }
751 }
752 }
753
754 for (int[] nextVector: vectorList) {
755 mdcm.set(get(nextVector), nextVector);
756 }
757 }
758 }
759
760
761 /** Computes the n<sup>th</sup> roots of unity.
762 * A cache of already computed values is maintained.
763 */
764 private static class RootsOfUnity implements Serializable {
765
766 /** Serializable version id. */
767 private static final long serialVersionUID = 6404784357747329667L;
768
769 /** Number of roots of unity. */
770 private int omegaCount;
771
772 /** Real part of the roots. */
773 private double[] omegaReal;
774
775 /** Imaginary part of the roots for forward transform. */
776 private double[] omegaImaginaryForward;
777
778 /** Imaginary part of the roots for reverse transform. */
779 private double[] omegaImaginaryInverse;
780
781 /** Forward/reverse indicator. */
782 private boolean isForward;
783
784 /**
785 * Build an engine for computing then <sup>th</sup> roots of unity
786 */
787 public RootsOfUnity() {
788
789 omegaCount = 0;
790 omegaReal = null;
791 omegaImaginaryForward = null;
792 omegaImaginaryInverse = null;
793 isForward = true;
794
795 }
796
797 /**
798 * Check if computation has been done for forward or reverse transform.
799 * @return true if computation has been done for forward transform
800 * @throws IllegalStateException if no roots of unity have been computed yet
801 */
802 public synchronized boolean isForward() throws IllegalStateException {
803
804 if (omegaCount == 0) {
805 throw MathRuntimeException.createIllegalStateException(LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
806 }
807 return isForward;
808
809 }
810
811 /** Computes the n<sup>th</sup> roots of unity.
812 * <p>The computed omega[] = { 1, w, w<sup>2</sup>, ... w<sup>(n-1)</sup> } where
813 * w = exp(-2 π i / n), i = &sqrt;(-1).</p>
814 * <p>Note that n is positive for
815 * forward transform and negative for inverse transform.</p>
816 * @param n number of roots of unity to compute,
817 * positive for forward transform, negative for inverse transform
818 * @throws IllegalArgumentException if n = 0
819 */
820 public synchronized void computeOmega(int n) throws IllegalArgumentException {
821
822 if (n == 0) {
823 throw MathRuntimeException.createIllegalArgumentException(
824 LocalizedFormats.CANNOT_COMPUTE_0TH_ROOT_OF_UNITY);
825 }
826
827 isForward = n > 0;
828
829 // avoid repetitive calculations
830 final int absN = FastMath.abs(n);
831
832 if (absN == omegaCount) {
833 return;
834 }
835
836 // calculate everything from scratch, for both forward and inverse versions
837 final double t = 2.0 * FastMath.PI / absN;
838 final double cosT = FastMath.cos(t);
839 final double sinT = FastMath.sin(t);
840 omegaReal = new double[absN];
841 omegaImaginaryForward = new double[absN];
842 omegaImaginaryInverse = new double[absN];
843 omegaReal[0] = 1.0;
844 omegaImaginaryForward[0] = 0.0;
845 omegaImaginaryInverse[0] = 0.0;
846 for (int i = 1; i < absN; i++) {
847 omegaReal[i] =
848 omegaReal[i-1] * cosT + omegaImaginaryForward[i-1] * sinT;
849 omegaImaginaryForward[i] =
850 omegaImaginaryForward[i-1] * cosT - omegaReal[i-1] * sinT;
851 omegaImaginaryInverse[i] = -omegaImaginaryForward[i];
852 }
853 omegaCount = absN;
854
855 }
856
857 /**
858 * Get the real part of the k<sup>th</sup> n<sup>th</sup> root of unity
859 * @param k index of the n<sup>th</sup> root of unity
860 * @return real part of the k<sup>th</sup> n<sup>th</sup> root of unity
861 * @throws IllegalStateException if no roots of unity have been computed yet
862 * @throws IllegalArgumentException if k is out of range
863 */
864 public synchronized double getOmegaReal(int k)
865 throws IllegalStateException, IllegalArgumentException {
866
867 if (omegaCount == 0) {
868 throw MathRuntimeException.createIllegalStateException(LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
869 }
870 if ((k < 0) || (k >= omegaCount)) {
871 throw MathRuntimeException.createIllegalArgumentException(
872 LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX, k, 0, omegaCount - 1);
873 }
874
875 return omegaReal[k];
876
877 }
878
879 /**
880 * Get the imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity
881 * @param k index of the n<sup>th</sup> root of unity
882 * @return imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity
883 * @throws IllegalStateException if no roots of unity have been computed yet
884 * @throws IllegalArgumentException if k is out of range
885 */
886 public synchronized double getOmegaImaginary(int k)
887 throws IllegalStateException, IllegalArgumentException {
888
889 if (omegaCount == 0) {
890 throw MathRuntimeException.createIllegalStateException(LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET);
891 }
892 if ((k < 0) || (k >= omegaCount)) {
893 throw MathRuntimeException.createIllegalArgumentException(
894 LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX, k, 0, omegaCount - 1);
895 }
896
897 return isForward ? omegaImaginaryForward[k] : omegaImaginaryInverse[k];
898
899 }
900
901 }
902
903 }