    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    }