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