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