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