View Javadoc

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.FunctionEvaluationException;
23  import org.apache.commons.math.MathRuntimeException;
24  import org.apache.commons.math.analysis.UnivariateRealFunction;
25  import org.apache.commons.math.complex.Complex;
26  
27  /**
28   * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html">
29   * Fast Fourier Transform</a> for transformation of one-dimensional data sets.
30   * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897,
31   * chapter 6.
32   * <p>
33   * There are several conventions for the definition of FFT and inverse FFT,
34   * mainly on different coefficient and exponent. Here the equations are listed
35   * in the comments of the corresponding methods.</p>
36   * <p>
37   * We require the length of data set to be power of 2, this greatly simplifies
38   * and speeds up the code. Users can pad the data with zeros to meet this
39   * requirement. There are other flavors of FFT, for reference, see S. Winograd,
40   * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation,
41   * 32 (1978), 175 - 199.</p>
42   *
43   * @version $Revision: 790243 $ $Date: 2009-07-01 12:03:28 -0400 (Wed, 01 Jul 2009) $
44   * @since 1.2
45   */
46  public class FastFourierTransformer implements Serializable {
47  
48      /** Serializable version identifier. */
49      static final long serialVersionUID = 5138259215438106000L;
50  
51      /** roots of unity */
52      private RootsOfUnity roots = new RootsOfUnity();
53  
54      /**
55       * Construct a default transformer.
56       */
57      public FastFourierTransformer() {
58          super();
59      }
60  
61      /**
62       * Transform the given real data set.
63       * <p>
64       * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
65       * </p>
66       * 
67       * @param f the real data array to be transformed
68       * @return the complex transformed array
69       * @throws IllegalArgumentException if any parameters are invalid
70       */
71      public Complex[] transform(double f[])
72          throws IllegalArgumentException {
73          return fft(f, false);
74      }
75  
76      /**
77       * Transform the given real function, sampled on the given interval.
78       * <p>
79       * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
80       * </p>
81       * 
82       * @param f the function to be sampled and transformed
83       * @param min the lower bound for the interval
84       * @param max the upper bound for the interval
85       * @param n the number of sample points
86       * @return the complex transformed array
87       * @throws FunctionEvaluationException if function cannot be evaluated
88       * at some point
89       * @throws IllegalArgumentException if any parameters are invalid
90       */
91      public Complex[] transform(UnivariateRealFunction f,
92                                 double min, double max, int n)
93          throws FunctionEvaluationException, IllegalArgumentException {
94          double data[] = sample(f, min, max, n);
95          return fft(data, false);
96      }
97  
98      /**
99       * 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&times;log(N)) complexity with
538      * N=n<sub>1</sub>&times;n<sub>2</sub>&times;n<sub>3</sub>&times;...&times;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 &pi; 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 }