001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.transform;
018
019 import java.io.Serializable;
020 import java.lang.reflect.Array;
021
022 import org.apache.commons.math.FunctionEvaluationException;
023 import org.apache.commons.math.MathRuntimeException;
024 import org.apache.commons.math.analysis.UnivariateRealFunction;
025 import org.apache.commons.math.complex.Complex;
026
027 /**
028 * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html">
029 * Fast Fourier Transform</a> for transformation of one-dimensional data sets.
030 * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897,
031 * chapter 6.
032 * <p>
033 * There are several conventions for the definition of FFT and inverse FFT,
034 * mainly on different coefficient and exponent. Here the equations are listed
035 * in the comments of the corresponding methods.</p>
036 * <p>
037 * We require the length of data set to be power of 2, this greatly simplifies
038 * and speeds up the code. Users can pad the data with zeros to meet this
039 * requirement. There are other flavors of FFT, for reference, see S. Winograd,
040 * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation,
041 * 32 (1978), 175 - 199.</p>
042 *
043 * @version $Revision: 790243 $ $Date: 2009-07-01 12:03:28 -0400 (Wed, 01 Jul 2009) $
044 * @since 1.2
045 */
046 public class FastFourierTransformer implements Serializable {
047
048 /** Serializable version identifier. */
049 static final long serialVersionUID = 5138259215438106000L;
050
051 /** roots of unity */
052 private RootsOfUnity roots = new RootsOfUnity();
053
054 /**
055 * Construct a default transformer.
056 */
057 public FastFourierTransformer() {
058 super();
059 }
060
061 /**
062 * Transform the given real data set.
063 * <p>
064 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
065 * </p>
066 *
067 * @param f the real data array to be transformed
068 * @return the complex transformed array
069 * @throws IllegalArgumentException if any parameters are invalid
070 */
071 public Complex[] transform(double f[])
072 throws IllegalArgumentException {
073 return fft(f, false);
074 }
075
076 /**
077 * Transform the given real function, sampled on the given interval.
078 * <p>
079 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
080 * </p>
081 *
082 * @param f the function to be sampled and transformed
083 * @param min the lower bound for the interval
084 * @param max the upper bound for the interval
085 * @param n the number of sample points
086 * @return the complex transformed array
087 * @throws FunctionEvaluationException if function cannot be evaluated
088 * at some point
089 * @throws IllegalArgumentException if any parameters are invalid
090 */
091 public Complex[] transform(UnivariateRealFunction f,
092 double min, double max, int n)
093 throws FunctionEvaluationException, IllegalArgumentException {
094 double data[] = sample(f, min, max, n);
095 return fft(data, false);
096 }
097
098 /**
099 * Transform the given complex data set.
100 * <p>
101 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
102 * </p>
103 *
104 * @param f the complex data array to be transformed
105 * @return the complex transformed array
106 * @throws IllegalArgumentException if any parameters are invalid
107 */
108 public Complex[] transform(Complex f[])
109 throws IllegalArgumentException {
110 roots.computeOmega(f.length);
111 return fft(f);
112 }
113
114 /**
115 * Transform the given real data set.
116 * <p>
117 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
118 * </p>
119 *
120 * @param f the real data array to be transformed
121 * @return the complex transformed array
122 * @throws IllegalArgumentException if any parameters are invalid
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 * Transform the given real function, sampled on the given interval.
133 * <p>
134 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
135 * </p>
136 *
137 * @param f the function to be sampled and transformed
138 * @param min the lower bound for the interval
139 * @param max the upper bound for the interval
140 * @param n the number of sample points
141 * @return the complex transformed array
142 * @throws FunctionEvaluationException if function cannot be evaluated
143 * at some point
144 * @throws IllegalArgumentException if any parameters are invalid
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 * Transform the given complex data set.
157 * <p>
158 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
159 * </p>
160 *
161 * @param f the complex data array to be transformed
162 * @return the complex transformed array
163 * @throws IllegalArgumentException if any parameters are invalid
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 * Inversely transform the given real data set.
175 * <p>
176 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
177 * </p>
178 *
179 * @param f the real data array to be inversely transformed
180 * @return the complex inversely transformed array
181 * @throws IllegalArgumentException if any parameters are invalid
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 * Inversely transform the given real function, sampled on the given interval.
192 * <p>
193 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
194 * </p>
195 *
196 * @param f the function to be sampled and inversely transformed
197 * @param min the lower bound for the interval
198 * @param max the upper bound for the interval
199 * @param n the number of sample points
200 * @return the complex inversely transformed array
201 * @throws FunctionEvaluationException if function cannot be evaluated
202 * at some point
203 * @throws IllegalArgumentException if any parameters are invalid
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 * Inversely transform the given complex data set.
216 * <p>
217 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
218 * </p>
219 *
220 * @param f the complex data array to be inversely transformed
221 * @return the complex inversely transformed array
222 * @throws IllegalArgumentException if any parameters are invalid
223 */
224 public Complex[] inversetransform(Complex f[])
225 throws IllegalArgumentException {
226
227 roots.computeOmega(-f.length); // pass negative argument
228 double scaling_coefficient = 1.0 / f.length;
229 return scaleArray(fft(f), scaling_coefficient);
230 }
231
232 /**
233 * Inversely transform the given real data set.
234 * <p>
235 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
236 * </p>
237 *
238 * @param f the real data array to be inversely transformed
239 * @return the complex inversely transformed array
240 * @throws IllegalArgumentException if any parameters are invalid
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 * Inversely transform the given real function, sampled on the given interval.
251 * <p>
252 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
253 * </p>
254 *
255 * @param f the function to be sampled and inversely transformed
256 * @param min the lower bound for the interval
257 * @param max the upper bound for the interval
258 * @param n the number of sample points
259 * @return the complex inversely transformed array
260 * @throws FunctionEvaluationException if function cannot be evaluated
261 * at some point
262 * @throws IllegalArgumentException if any parameters are invalid
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 * Inversely transform the given complex data set.
275 * <p>
276 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
277 * </p>
278 *
279 * @param f the complex data array to be inversely transformed
280 * @return the complex inversely transformed array
281 * @throws IllegalArgumentException if any parameters are invalid
282 */
283 public Complex[] inversetransform2(Complex f[])
284 throws IllegalArgumentException {
285
286 roots.computeOmega(-f.length); // pass negative argument
287 double scaling_coefficient = 1.0 / Math.sqrt(f.length);
288 return scaleArray(fft(f), scaling_coefficient);
289 }
290
291 /**
292 * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
293 *
294 * @param f the real data array to be transformed
295 * @param isInverse the indicator of forward or inverse transform
296 * @return the complex transformed array
297 * @throws IllegalArgumentException if any parameters are invalid
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 // Rather than the naive real to complex conversion, pack 2N
310 // real numbers into N complex numbers for better performance.
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 // reconstruct the FFT result for the original array
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 //Complex D = roots.getOmega(i).multiply(Complex.I);
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 * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
339 *
340 * @param data the complex data array to be transformed
341 * @return the complex transformed array
342 * @throws IllegalArgumentException if any parameters are invalid
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 // initial simple cases
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 // permute original data array in bit-reversal order
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 // the bottom base-4 round
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 // omegaCount indicates forward or inverse transform
384 f[i+1] = roots.isForward() ? F : E;
385 f[i+3] = roots.isForward() ? E : F;
386 }
387
388 // iterations from bottom to top take O(N*logN) time
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 //z = f[i+j+k].multiply(roots.getOmega(k*m));
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 //z = f[i+j+k].multiply(omega[k*m]);
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 * Sample the given univariate real function on the given interval.
414 * <p>
415 * The interval is divided equally into N sections and sample points
416 * are taken from min to max-(max-min)/N. Usually f(x) is periodic
417 * such that f(min) = f(max) (note max is not sampled), but we don't
418 * require that.</p>
419 *
420 * @param f the function to be sampled
421 * @param min the lower bound for the interval
422 * @param max the upper bound for the interval
423 * @param n the number of sample points
424 * @return the samples array
425 * @throws FunctionEvaluationException if function cannot be evaluated
426 * at some point
427 * @throws IllegalArgumentException if any parameters are invalid
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 * Multiply every component in the given real array by the
450 * given real number. The change is made in place.
451 *
452 * @param f the real array to be scaled
453 * @param d the real scaling coefficient
454 * @return a reference to the scaled array
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 * Multiply every component in the given complex array by the
465 * given real number. The change is made in place.
466 *
467 * @param f the complex array to be scaled
468 * @param d the real scaling coefficient
469 * @return a reference to the scaled array
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 * Returns true if the argument is power of 2.
480 *
481 * @param n the number to test
482 * @return true if the argument is power of 2
483 */
484 public static boolean isPowerOf2(long n) {
485 return (n > 0) && ((n & (n - 1)) == 0);
486 }
487
488 /**
489 * Verifies that the data set has length of power of 2.
490 *
491 * @param d the data array
492 * @throws IllegalArgumentException if array length is not power of 2
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 * Verifies that the data set has length of power of 2.
504 *
505 * @param o the data array
506 * @throws IllegalArgumentException if array length is not power of 2
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 * Verifies that the endpoints specify an interval.
518 *
519 * @param lower lower endpoint
520 * @param upper upper endpoint
521 * @throws IllegalArgumentException if not interval
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 * Performs a multi-dimensional Fourier transform on a given array.
535 * Use {@link #inversetransform2(Complex[])} and
536 * {@link #transform2(Complex[])} in a row-column implementation
537 * in any number of dimensions with O(N×log(N)) complexity with
538 * N=n<sub>1</sub>×n<sub>2</sub>×n<sub>3</sub>×...×n<sub>d</sub>,
539 * n<sub>x</sub>=number of elements in dimension x,
540 * and d=total number of dimensions.
541 *
542 * @param mdca Multi-Dimensional Complex Array id est Complex[][][][]
543 * @param forward inverseTransform2 is preformed if this is false
544 * @return transform of mdca as a Multi-Dimensional Complex Array id est Complex[][][][]
545 * @throws IllegalArgumentException if any dimension is not a power of two
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 //cycle through each dimension
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 * Performs one dimension of a multi-dimensional Fourier transform.
561 *
562 * @param mdcm input matrix
563 * @param forward inverseTransform2 is preformed if this is false
564 * @param d index of the dimension to process
565 * @param subVector recursion subvector
566 * @throws IllegalArgumentException if any dimension is not a power of two
567 */
568 private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward,
569 int d, int[] subVector)
570 throws IllegalArgumentException {
571 int[] dimensionSize = mdcm.getDimensionSizes();
572 //if done
573 if (subVector.length == dimensionSize.length) {
574 Complex[] temp = new Complex[dimensionSize[d]];
575 for (int i = 0; i < dimensionSize[d]; i++) {
576 //fft along dimension d
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 //value is not important once the recursion is done.
595 //then an fft will be applied along the dimension d.
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 //further split along the next dimension
602 mdfft(mdcm, forward, d, vector);
603 }
604 }
605 }
606 return;
607 }
608
609 /**
610 * Complex matrix implementation.
611 * Not designed for synchronized access
612 * may eventually be replaced by jsr-83 of the java community process
613 * http://jcp.org/en/jsr/detail?id=83
614 * may require additional exception throws for other basic requirements.
615 */
616 private static class MultiDimensionalComplexMatrix
617 implements Cloneable {
618
619 /** Size in all dimensions. */
620 protected int[] dimensionSize;
621
622 /** Storage array. */
623 protected Object multiDimensionalComplexArray;
624
625 /** Simple constructor.
626 * @param multiDimensionalComplexArray array containing the matrix elements
627 */
628 public MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) {
629
630 this.multiDimensionalComplexArray = multiDimensionalComplexArray;
631
632 // count dimensions
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 // allocate array with exact count
642 dimensionSize = new int[numOfDimensions];
643
644 // fill array
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 * Get a matrix element.
657 * @param vector indices of the element
658 * @return matrix element
659 * @exception IllegalArgumentException if dimensions do not match
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 * Set a matrix element.
687 * @param magnitude magnitude of the element
688 * @param vector indices of the element
689 * @return the previous value
690 * @exception IllegalArgumentException if dimensions do not match
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 * Get the size in all dimensions.
721 * @return size in all dimensions
722 */
723 public int[] getDimensionSizes() {
724 return dimensionSize.clone();
725 }
726
727 /**
728 * Get the underlying storage array
729 * @return underlying storage array
730 */
731 public Object getArray() {
732 return multiDimensionalComplexArray;
733 }
734
735 /** {@inheritDoc} */
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 * Copy contents of current array into mdcm.
747 * @param mdcm array where to copy data
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 /** Computes the n<sup>th</sup> roots of unity.
777 * A cache of already computed values is maintained.
778 */
779 private static class RootsOfUnity implements Serializable {
780
781 /** Serializable version id. */
782 private static final long serialVersionUID = 6404784357747329667L;
783
784 /** Number of roots of unity. */
785 private int omegaCount;
786
787 /** Real part of the roots. */
788 private double[] omegaReal;
789
790 /** Imaginary part of the roots for forward transform. */
791 private double[] omegaImaginaryForward;
792
793 /** Imaginary part of the roots for reverse transform. */
794 private double[] omegaImaginaryInverse;
795
796 /** Forward/reverse indicator. */
797 private boolean isForward;
798
799 /**
800 * Build an engine for computing then <sup>th</sup> roots of unity
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 * Check if computation has been done for forward or reverse transform.
814 * @return true if computation has been done for forward transform
815 * @throws IllegalStateException if no roots of unity have been computed yet
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 /** Computes the n<sup>th</sup> roots of unity.
828 * <p>The computed omega[] = { 1, w, w<sup>2</sup>, ... w<sup>(n-1)</sup> } where
829 * w = exp(-2 π i / n), i = &sqrt;(-1).</p>
830 * <p>Note that n is positive for
831 * forward transform and negative for inverse transform.</p>
832 * @param n number of roots of unity to compute,
833 * positive for forward transform, negative for inverse transform
834 * @throws IllegalArgumentException if n = 0
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 // avoid repetitive calculations
846 final int absN = Math.abs(n);
847
848 if (absN == omegaCount) {
849 return;
850 }
851
852 // calculate everything from scratch, for both forward and inverse versions
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 * Get the real part of the k<sup>th</sup> n<sup>th</sup> root of unity
875 * @param k index of the n<sup>th</sup> root of unity
876 * @return real part of the k<sup>th</sup> n<sup>th</sup> root of unity
877 * @throws IllegalStateException if no roots of unity have been computed yet
878 * @throws IllegalArgumentException if k is out of range
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 * Get the imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity
899 * @param k index of the n<sup>th</sup> root of unity
900 * @return imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity
901 * @throws IllegalStateException if no roots of unity have been computed yet
902 * @throws IllegalArgumentException if k is out of range
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 }