001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.transform;
018    
019    import java.io.Serializable;
020    import java.lang.reflect.Array;
021    
022    import org.apache.commons.math.FunctionEvaluationException;
023    import org.apache.commons.math.MathRuntimeException;
024    import org.apache.commons.math.analysis.UnivariateRealFunction;
025    import org.apache.commons.math.complex.Complex;
026    
027    /**
028     * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html">
029     * Fast Fourier Transform</a> for transformation of one-dimensional data sets.
030     * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897,
031     * chapter 6.
032     * <p>
033     * There are several conventions for the definition of FFT and inverse FFT,
034     * mainly on different coefficient and exponent. Here the equations are listed
035     * in the comments of the corresponding methods.</p>
036     * <p>
037     * We require the length of data set to be power of 2, this greatly simplifies
038     * and speeds up the code. Users can pad the data with zeros to meet this
039     * requirement. There are other flavors of FFT, for reference, see S. Winograd,
040     * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation,
041     * 32 (1978), 175 - 199.</p>
042     *
043     * @version $Revision: 790243 $ $Date: 2009-07-01 12:03:28 -0400 (Wed, 01 Jul 2009) $
044     * @since 1.2
045     */
046    public class FastFourierTransformer implements Serializable {
047    
048        /** Serializable version identifier. */
049        static final long serialVersionUID = 5138259215438106000L;
050    
051        /** roots of unity */
052        private RootsOfUnity roots = new RootsOfUnity();
053    
054        /**
055         * Construct a default transformer.
056         */
057        public FastFourierTransformer() {
058            super();
059        }
060    
061        /**
062         * Transform the given real data set.
063         * <p>
064         * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
065         * </p>
066         * 
067         * @param f the real data array to be transformed
068         * @return the complex transformed array
069         * @throws IllegalArgumentException if any parameters are invalid
070         */
071        public Complex[] transform(double f[])
072            throws IllegalArgumentException {
073            return fft(f, false);
074        }
075    
076        /**
077         * Transform the given real function, sampled on the given interval.
078         * <p>
079         * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
080         * </p>
081         * 
082         * @param f the function to be sampled and transformed
083         * @param min the lower bound for the interval
084         * @param max the upper bound for the interval
085         * @param n the number of sample points
086         * @return the complex transformed array
087         * @throws FunctionEvaluationException if function cannot be evaluated
088         * at some point
089         * @throws IllegalArgumentException if any parameters are invalid
090         */
091        public Complex[] transform(UnivariateRealFunction f,
092                                   double min, double max, int n)
093            throws FunctionEvaluationException, IllegalArgumentException {
094            double data[] = sample(f, min, max, n);
095            return fft(data, false);
096        }
097    
098        /**
099         * Transform the given complex data set.
100         * <p>
101         * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $
102         * </p>
103         * 
104         * @param f the complex data array to be transformed
105         * @return the complex transformed array
106         * @throws IllegalArgumentException if any parameters are invalid
107         */
108        public Complex[] transform(Complex f[])
109            throws IllegalArgumentException {
110            roots.computeOmega(f.length);
111            return fft(f);
112        }
113    
114        /**
115         * Transform the given real data set.
116         * <p>
117         * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
118         * </p>
119         * 
120         * @param f the real data array to be transformed
121         * @return the complex transformed array
122         * @throws IllegalArgumentException if any parameters are invalid
123         */
124        public Complex[] transform2(double f[])
125            throws IllegalArgumentException {
126    
127            double scaling_coefficient = 1.0 / Math.sqrt(f.length);
128            return scaleArray(fft(f, false), scaling_coefficient);
129        }
130    
131        /**
132         * Transform the given real function, sampled on the given interval.
133         * <p>
134         * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
135         * </p>
136         * 
137         * @param f the function to be sampled and transformed
138         * @param min the lower bound for the interval
139         * @param max the upper bound for the interval
140         * @param n the number of sample points
141         * @return the complex transformed array
142         * @throws FunctionEvaluationException if function cannot be evaluated
143         * at some point
144         * @throws IllegalArgumentException if any parameters are invalid
145         */
146        public Complex[] transform2(UnivariateRealFunction f,
147                                    double min, double max, int n)
148            throws FunctionEvaluationException, IllegalArgumentException {
149    
150            double data[] = sample(f, min, max, n);
151            double scaling_coefficient = 1.0 / Math.sqrt(n);
152            return scaleArray(fft(data, false), scaling_coefficient);
153        }
154    
155        /**
156         * Transform the given complex data set.
157         * <p>
158         * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$
159         * </p>
160         * 
161         * @param f the complex data array to be transformed
162         * @return the complex transformed array
163         * @throws IllegalArgumentException if any parameters are invalid
164         */
165        public Complex[] transform2(Complex f[])
166            throws IllegalArgumentException {
167    
168            roots.computeOmega(f.length);
169            double scaling_coefficient = 1.0 / Math.sqrt(f.length);
170            return scaleArray(fft(f), scaling_coefficient);
171        }
172    
173        /**
174         * Inversely transform the given real data set.
175         * <p>
176         * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
177         * </p>
178         * 
179         * @param f the real data array to be inversely transformed
180         * @return the complex inversely transformed array
181         * @throws IllegalArgumentException if any parameters are invalid
182         */
183        public Complex[] inversetransform(double f[])
184            throws IllegalArgumentException {
185    
186            double scaling_coefficient = 1.0 / f.length;
187            return scaleArray(fft(f, true), scaling_coefficient);
188        }
189    
190        /**
191         * Inversely transform the given real function, sampled on the given interval.
192         * <p>
193         * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
194         * </p>
195         * 
196         * @param f the function to be sampled and inversely transformed
197         * @param min the lower bound for the interval
198         * @param max the upper bound for the interval
199         * @param n the number of sample points
200         * @return the complex inversely transformed array
201         * @throws FunctionEvaluationException if function cannot be evaluated
202         * at some point
203         * @throws IllegalArgumentException if any parameters are invalid
204         */
205        public Complex[] inversetransform(UnivariateRealFunction f,
206                                          double min, double max, int n)
207            throws FunctionEvaluationException, IllegalArgumentException {
208    
209            double data[] = sample(f, min, max, n);
210            double scaling_coefficient = 1.0 / n;
211            return scaleArray(fft(data, true), scaling_coefficient);
212        }
213    
214        /**
215         * Inversely transform the given complex data set.
216         * <p>
217         * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $
218         * </p>
219         * 
220         * @param f the complex data array to be inversely transformed
221         * @return the complex inversely transformed array
222         * @throws IllegalArgumentException if any parameters are invalid
223         */
224        public Complex[] inversetransform(Complex f[])
225            throws IllegalArgumentException {
226    
227            roots.computeOmega(-f.length);    // pass negative argument
228            double scaling_coefficient = 1.0 / f.length;
229            return scaleArray(fft(f), scaling_coefficient);
230        }
231    
232        /**
233         * Inversely transform the given real data set.
234         * <p>
235         * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
236         * </p>
237         * 
238         * @param f the real data array to be inversely transformed
239         * @return the complex inversely transformed array
240         * @throws IllegalArgumentException if any parameters are invalid
241         */
242        public Complex[] inversetransform2(double f[])
243            throws IllegalArgumentException {
244    
245            double scaling_coefficient = 1.0 / Math.sqrt(f.length);
246            return scaleArray(fft(f, true), scaling_coefficient);
247        }
248    
249        /**
250         * Inversely transform the given real function, sampled on the given interval.
251         * <p>
252         * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
253         * </p>
254         * 
255         * @param f the function to be sampled and inversely transformed
256         * @param min the lower bound for the interval
257         * @param max the upper bound for the interval
258         * @param n the number of sample points
259         * @return the complex inversely transformed array
260         * @throws FunctionEvaluationException if function cannot be evaluated
261         * at some point
262         * @throws IllegalArgumentException if any parameters are invalid
263         */
264        public Complex[] inversetransform2(UnivariateRealFunction f,
265                                           double min, double max, int n)
266            throws FunctionEvaluationException, IllegalArgumentException {
267    
268            double data[] = sample(f, min, max, n);
269            double scaling_coefficient = 1.0 / Math.sqrt(n);
270            return scaleArray(fft(data, true), scaling_coefficient);
271        }
272    
273        /**
274         * Inversely transform the given complex data set.
275         * <p>
276         * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$
277         * </p>
278         * 
279         * @param f the complex data array to be inversely transformed
280         * @return the complex inversely transformed array
281         * @throws IllegalArgumentException if any parameters are invalid
282         */
283        public Complex[] inversetransform2(Complex f[])
284            throws IllegalArgumentException {
285    
286            roots.computeOmega(-f.length);    // pass negative argument
287            double scaling_coefficient = 1.0 / Math.sqrt(f.length);
288            return scaleArray(fft(f), scaling_coefficient);
289        }
290    
291        /**
292         * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
293         *
294         * @param f the real data array to be transformed
295         * @param isInverse the indicator of forward or inverse transform
296         * @return the complex transformed array
297         * @throws IllegalArgumentException if any parameters are invalid
298         */
299        protected Complex[] fft(double f[], boolean isInverse)
300            throws IllegalArgumentException {
301    
302            verifyDataSet(f);
303            Complex F[] = new Complex[f.length];
304            if (f.length == 1) {
305                F[0] = new Complex(f[0], 0.0);
306                return F;
307            }
308    
309            // Rather than the naive real to complex conversion, pack 2N
310            // real numbers into N complex numbers for better performance.
311            int N = f.length >> 1;
312            Complex c[] = new Complex[N];
313            for (int i = 0; i < N; i++) {
314                c[i] = new Complex(f[2*i], f[2*i+1]);
315            }
316            roots.computeOmega(isInverse ? -N : N);
317            Complex z[] = fft(c);
318    
319            // reconstruct the FFT result for the original array
320            roots.computeOmega(isInverse ? -2*N : 2*N);
321            F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);
322            F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0);
323            for (int i = 1; i < N; i++) {
324                Complex A = z[N-i].conjugate();
325                Complex B = z[i].add(A);
326                Complex C = z[i].subtract(A);
327                //Complex D = roots.getOmega(i).multiply(Complex.I);
328                Complex D = new Complex(-roots.getOmegaImaginary(i),
329                                        roots.getOmegaReal(i));
330                F[i] = B.subtract(C.multiply(D));
331                F[2*N-i] = F[i].conjugate();
332            }
333    
334            return scaleArray(F, 0.5);
335        }
336    
337        /**
338         * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).
339         *
340         * @param data the complex data array to be transformed
341         * @return the complex transformed array
342         * @throws IllegalArgumentException if any parameters are invalid
343         */
344        protected Complex[] fft(Complex data[])
345            throws IllegalArgumentException {
346    
347            int i, j, k, m, N = data.length;
348            Complex A, B, C, D, E, F, z, f[] = new Complex[N];
349    
350            // initial simple cases
351            verifyDataSet(data);
352            if (N == 1) {
353                f[0] = data[0];
354                return f;
355            }
356            if (N == 2) {
357                f[0] = data[0].add(data[1]);
358                f[1] = data[0].subtract(data[1]);
359                return f;
360            }
361    
362            // permute original data array in bit-reversal order
363            j = 0;
364            for (i = 0; i < N; i++) {
365                f[i] = data[j];
366                k = N >> 1;
367                while (j >= k && k > 0) {
368                    j -= k; k >>= 1;
369                }
370                j += k;
371            }
372    
373            // the bottom base-4 round
374            for (i = 0; i < N; i += 4) {
375                A = f[i].add(f[i+1]);
376                B = f[i+2].add(f[i+3]);
377                C = f[i].subtract(f[i+1]);
378                D = f[i+2].subtract(f[i+3]);
379                E = C.add(D.multiply(Complex.I));
380                F = C.subtract(D.multiply(Complex.I));
381                f[i] = A.add(B);
382                f[i+2] = A.subtract(B);
383                // omegaCount indicates forward or inverse transform
384                f[i+1] = roots.isForward() ? F : E;
385                f[i+3] = roots.isForward() ? E : F;
386            }
387    
388            // iterations from bottom to top take O(N*logN) time
389            for (i = 4; i < N; i <<= 1) {
390                m = N / (i<<1);
391                for (j = 0; j < N; j += i<<1) {
392                    for (k = 0; k < i; k++) {
393                        //z = f[i+j+k].multiply(roots.getOmega(k*m));
394                        final int k_times_m = k*m;
395                        final double omega_k_times_m_real = roots.getOmegaReal(k_times_m);
396                        final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m);
397                        //z = f[i+j+k].multiply(omega[k*m]);
398                        z = new Complex(
399                            f[i+j+k].getReal() * omega_k_times_m_real -
400                            f[i+j+k].getImaginary() * omega_k_times_m_imaginary,
401                            f[i+j+k].getReal() * omega_k_times_m_imaginary +
402                            f[i+j+k].getImaginary() * omega_k_times_m_real);
403                      
404                        f[i+j+k] = f[j+k].subtract(z);
405                        f[j+k] = f[j+k].add(z);
406                    }
407                }
408            }
409            return f;
410        }
411    
412        /**
413         * Sample the given univariate real function on the given interval.
414         * <p>
415         * The interval is divided equally into N sections and sample points
416         * are taken from min to max-(max-min)/N. Usually f(x) is periodic
417         * such that f(min) = f(max) (note max is not sampled), but we don't
418         * require that.</p>
419         *
420         * @param f the function to be sampled
421         * @param min the lower bound for the interval
422         * @param max the upper bound for the interval
423         * @param n the number of sample points
424         * @return the samples array
425         * @throws FunctionEvaluationException if function cannot be evaluated
426         * at some point
427         * @throws IllegalArgumentException if any parameters are invalid
428         */
429        public static double[] sample(UnivariateRealFunction f,
430                                      double min, double max, int n)
431            throws FunctionEvaluationException, IllegalArgumentException {
432    
433            if (n <= 0) {
434                throw MathRuntimeException.createIllegalArgumentException(
435                        "number of sample is not positive: {0}",
436                        n);
437            }
438            verifyInterval(min, max);
439    
440            double s[] = new double[n];
441            double h = (max - min) / n;
442            for (int i = 0; i < n; i++) {
443                s[i] = f.value(min + i * h);
444            }
445            return s;
446        }
447    
448        /**
449         * Multiply every component in the given real array by the
450         * given real number. The change is made in place.
451         *
452         * @param f the real array to be scaled
453         * @param d the real scaling coefficient
454         * @return a reference to the scaled array
455         */
456        public static double[] scaleArray(double f[], double d) {
457            for (int i = 0; i < f.length; i++) {
458                f[i] *= d;
459            }
460            return f;
461        }
462    
463        /**
464         * Multiply every component in the given complex array by the
465         * given real number. The change is made in place.
466         *
467         * @param f the complex array to be scaled
468         * @param d the real scaling coefficient
469         * @return a reference to the scaled array
470         */
471        public static Complex[] scaleArray(Complex f[], double d) {
472            for (int i = 0; i < f.length; i++) {
473                f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary());
474            }
475            return f;
476        }
477    
478        /**
479         * Returns true if the argument is power of 2.
480         * 
481         * @param n the number to test
482         * @return true if the argument is power of 2
483         */
484        public static boolean isPowerOf2(long n) {
485            return (n > 0) && ((n & (n - 1)) == 0);
486        }
487    
488        /**
489         * Verifies that the data set has length of power of 2.
490         * 
491         * @param d the data array
492         * @throws IllegalArgumentException if array length is not power of 2
493         */
494        public static void verifyDataSet(double d[]) throws IllegalArgumentException {
495            if (!isPowerOf2(d.length)) {
496                throw MathRuntimeException.createIllegalArgumentException(
497                        "{0} is not a power of 2, consider padding for fix",
498                        d.length);
499            }       
500        }
501    
502        /**
503         * Verifies that the data set has length of power of 2.
504         * 
505         * @param o the data array
506         * @throws IllegalArgumentException if array length is not power of 2
507         */
508        public static void verifyDataSet(Object o[]) throws IllegalArgumentException {
509            if (!isPowerOf2(o.length)) {
510                throw MathRuntimeException.createIllegalArgumentException(
511                        "{0} is not a power of 2, consider padding for fix",
512                        o.length);
513            }       
514        }
515    
516        /**
517         * Verifies that the endpoints specify an interval.
518         * 
519         * @param lower lower endpoint
520         * @param upper upper endpoint
521         * @throws IllegalArgumentException if not interval
522         */
523        public static void verifyInterval(double lower, double upper)
524            throws IllegalArgumentException {
525    
526            if (lower >= upper) {
527                throw MathRuntimeException.createIllegalArgumentException(
528                        "endpoints do not specify an interval: [{0}, {1}]",
529                        lower, upper);
530            }       
531        }
532        
533        /**
534         * Performs a multi-dimensional Fourier transform on a given array.
535         * Use {@link #inversetransform2(Complex[])} and
536         * {@link #transform2(Complex[])} in a row-column implementation
537         * in any number of dimensions with O(N&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    }