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.math3.transform;
018    
019    import java.io.Serializable;
020    import java.lang.reflect.Array;
021    
022    import org.apache.commons.math3.analysis.FunctionUtils;
023    import org.apache.commons.math3.analysis.UnivariateFunction;
024    import org.apache.commons.math3.complex.Complex;
025    import org.apache.commons.math3.exception.DimensionMismatchException;
026    import org.apache.commons.math3.exception.MathIllegalArgumentException;
027    import org.apache.commons.math3.exception.MathIllegalStateException;
028    import org.apache.commons.math3.exception.util.LocalizedFormats;
029    import org.apache.commons.math3.util.ArithmeticUtils;
030    import org.apache.commons.math3.util.FastMath;
031    import org.apache.commons.math3.util.MathArrays;
032    
033    /**
034     * Implements the Fast Fourier Transform for transformation of one-dimensional
035     * real or complex data sets. For reference, see <em>Applied Numerical Linear
036     * Algebra</em>, ISBN 0898713897, chapter 6.
037     * <p>
038     * There are several variants of the discrete Fourier transform, with various
039     * normalization conventions, which are specified by the parameter
040     * {@link DftNormalization}.
041     * <p>
042     * The current implementation of the discrete Fourier transform as a fast
043     * Fourier transform requires the length of the data set to be a power of 2.
044     * This greatly simplifies and speeds up the code. Users can pad the data with
045     * zeros to meet this requirement. There are other flavors of FFT, for
046     * reference, see S. Winograd,
047     * <i>On computing the discrete Fourier transform</i>, Mathematics of
048     * Computation, 32 (1978), 175 - 199.
049     *
050     * @see DftNormalization
051     * @version $Id: FastFourierTransformer.java 1385310 2012-09-16 16:32:10Z tn $
052     * @since 1.2
053     */
054    public class FastFourierTransformer implements Serializable {
055    
056        /** Serializable version identifier. */
057        static final long serialVersionUID = 20120210L;
058    
059        /**
060         * {@code W_SUB_N_R[i]} is the real part of
061         * {@code exp(- 2 * i * pi / n)}:
062         * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}.
063         */
064        private static final double[] W_SUB_N_R =
065                {  0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1
066                , 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1
067                , 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1
068                , 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1
069                , 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1
070                , 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1
071                , 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1
072                , 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0
073                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
074                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
075                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
076                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
077                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
078                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
079                , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
080                , 0x1.0p0, 0x1.0p0, 0x1.0p0 };
081    
082        /**
083         * {@code W_SUB_N_I[i]} is the imaginary part of
084         * {@code exp(- 2 * i * pi / n)}:
085         * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}.
086         */
087        private static final double[] W_SUB_N_I =
088                {  0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1
089                , -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5
090                , -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9
091                , -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13
092                , -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17
093                , -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21
094                , -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25
095                , -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29
096                , -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33
097                , -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37
098                , -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41
099                , -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45
100                , -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49
101                , -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53
102                , -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57
103                , -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 };
104    
105        /** The type of DFT to be performed. */
106        private final DftNormalization normalization;
107    
108        /**
109         * Creates a new instance of this class, with various normalization
110         * conventions.
111         *
112         * @param normalization the type of normalization to be applied to the
113         * transformed data
114         */
115        public FastFourierTransformer(final DftNormalization normalization) {
116            this.normalization = normalization;
117        }
118    
119        /**
120         * Performs identical index bit reversal shuffles on two arrays of identical
121         * size. Each element in the array is swapped with another element based on
122         * the bit-reversal of the index. For example, in an array with length 16,
123         * item at binary index 0011 (decimal 3) would be swapped with the item at
124         * binary index 1100 (decimal 12).
125         *
126         * @param a the first array to be shuffled
127         * @param b the second array to be shuffled
128         */
129        private static void bitReversalShuffle2(double[] a, double[] b) {
130            final int n = a.length;
131            assert b.length == n;
132            final int halfOfN = n >> 1;
133    
134            int j = 0;
135            for (int i = 0; i < n; i++) {
136                if (i < j) {
137                    // swap indices i & j
138                    double temp = a[i];
139                    a[i] = a[j];
140                    a[j] = temp;
141    
142                    temp = b[i];
143                    b[i] = b[j];
144                    b[j] = temp;
145                }
146    
147                int k = halfOfN;
148                while (k <= j && k > 0) {
149                    j -= k;
150                    k >>= 1;
151                }
152                j += k;
153            }
154        }
155    
156        /**
157         * Applies the proper normalization to the specified transformed data.
158         *
159         * @param dataRI the unscaled transformed data
160         * @param normalization the normalization to be applied
161         * @param type the type of transform (forward, inverse) which resulted in the specified data
162         */
163        private static void normalizeTransformedData(final double[][] dataRI,
164            final DftNormalization normalization, final TransformType type) {
165    
166            final double[] dataR = dataRI[0];
167            final double[] dataI = dataRI[1];
168            final int n = dataR.length;
169            assert dataI.length == n;
170    
171            switch (normalization) {
172                case STANDARD:
173                    if (type == TransformType.INVERSE) {
174                        final double scaleFactor = 1.0 / ((double) n);
175                        for (int i = 0; i < n; i++) {
176                            dataR[i] *= scaleFactor;
177                            dataI[i] *= scaleFactor;
178                        }
179                    }
180                    break;
181                case UNITARY:
182                    final double scaleFactor = 1.0 / FastMath.sqrt(n);
183                    for (int i = 0; i < n; i++) {
184                        dataR[i] *= scaleFactor;
185                        dataI[i] *= scaleFactor;
186                    }
187                    break;
188                default:
189                    /*
190                     * This should never occur in normal conditions. However this
191                     * clause has been added as a safeguard if other types of
192                     * normalizations are ever implemented, and the corresponding
193                     * test is forgotten in the present switch.
194                     */
195                    throw new MathIllegalStateException();
196            }
197        }
198    
199        /**
200         * Computes the standard transform of the specified complex data. The
201         * computation is done in place. The input data is laid out as follows
202         * <ul>
203         *   <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point,</li>
204         *   <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data point.</li>
205         * </ul>
206         *
207         * @param dataRI the two dimensional array of real and imaginary parts of the data
208         * @param normalization the normalization to be applied to the transformed data
209         * @param type the type of transform (forward, inverse) to be performed
210         * @throws DimensionMismatchException if the number of rows of the specified
211         *   array is not two, or the array is not rectangular
212         * @throws MathIllegalArgumentException if the number of data points is not
213         *   a power of two
214         */
215        public static void transformInPlace(final double[][] dataRI,
216            final DftNormalization normalization, final TransformType type) {
217    
218            if (dataRI.length != 2) {
219                throw new DimensionMismatchException(dataRI.length, 2);
220            }
221            final double[] dataR = dataRI[0];
222            final double[] dataI = dataRI[1];
223            if (dataR.length != dataI.length) {
224                throw new DimensionMismatchException(dataI.length, dataR.length);
225            }
226    
227            final int n = dataR.length;
228            if (!ArithmeticUtils.isPowerOfTwo(n)) {
229                throw new MathIllegalArgumentException(
230                    LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
231                    Integer.valueOf(n));
232            }
233    
234            if (n == 1) {
235                return;
236            } else if (n == 2) {
237                final double srcR0 = dataR[0];
238                final double srcI0 = dataI[0];
239                final double srcR1 = dataR[1];
240                final double srcI1 = dataI[1];
241    
242                // X_0 = x_0 + x_1
243                dataR[0] = srcR0 + srcR1;
244                dataI[0] = srcI0 + srcI1;
245                // X_1 = x_0 - x_1
246                dataR[1] = srcR0 - srcR1;
247                dataI[1] = srcI0 - srcI1;
248    
249                normalizeTransformedData(dataRI, normalization, type);
250                return;
251            }
252    
253            bitReversalShuffle2(dataR, dataI);
254    
255            // Do 4-term DFT.
256            if (type == TransformType.INVERSE) {
257                for (int i0 = 0; i0 < n; i0 += 4) {
258                    final int i1 = i0 + 1;
259                    final int i2 = i0 + 2;
260                    final int i3 = i0 + 3;
261    
262                    final double srcR0 = dataR[i0];
263                    final double srcI0 = dataI[i0];
264                    final double srcR1 = dataR[i2];
265                    final double srcI1 = dataI[i2];
266                    final double srcR2 = dataR[i1];
267                    final double srcI2 = dataI[i1];
268                    final double srcR3 = dataR[i3];
269                    final double srcI3 = dataI[i3];
270    
271                    // 4-term DFT
272                    // X_0 = x_0 + x_1 + x_2 + x_3
273                    dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
274                    dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
275                    // X_1 = x_0 - x_2 + j * (x_3 - x_1)
276                    dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1);
277                    dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3);
278                    // X_2 = x_0 - x_1 + x_2 - x_3
279                    dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
280                    dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
281                    // X_3 = x_0 - x_2 + j * (x_1 - x_3)
282                    dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3);
283                    dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1);
284                }
285            } else {
286                for (int i0 = 0; i0 < n; i0 += 4) {
287                    final int i1 = i0 + 1;
288                    final int i2 = i0 + 2;
289                    final int i3 = i0 + 3;
290    
291                    final double srcR0 = dataR[i0];
292                    final double srcI0 = dataI[i0];
293                    final double srcR1 = dataR[i2];
294                    final double srcI1 = dataI[i2];
295                    final double srcR2 = dataR[i1];
296                    final double srcI2 = dataI[i1];
297                    final double srcR3 = dataR[i3];
298                    final double srcI3 = dataI[i3];
299    
300                    // 4-term DFT
301                    // X_0 = x_0 + x_1 + x_2 + x_3
302                    dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
303                    dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
304                    // X_1 = x_0 - x_2 + j * (x_3 - x_1)
305                    dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3);
306                    dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1);
307                    // X_2 = x_0 - x_1 + x_2 - x_3
308                    dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
309                    dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
310                    // X_3 = x_0 - x_2 + j * (x_1 - x_3)
311                    dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1);
312                    dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3);
313                }
314            }
315    
316            int lastN0 = 4;
317            int lastLogN0 = 2;
318            while (lastN0 < n) {
319                int n0 = lastN0 << 1;
320                int logN0 = lastLogN0 + 1;
321                double wSubN0R = W_SUB_N_R[logN0];
322                double wSubN0I = W_SUB_N_I[logN0];
323                if (type == TransformType.INVERSE) {
324                    wSubN0I = -wSubN0I;
325                }
326    
327                // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2).
328                for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) {
329                    int destOddStartIndex = destEvenStartIndex + lastN0;
330    
331                    double wSubN0ToRR = 1;
332                    double wSubN0ToRI = 0;
333    
334                    for (int r = 0; r < lastN0; r++) {
335                        double grR = dataR[destEvenStartIndex + r];
336                        double grI = dataI[destEvenStartIndex + r];
337                        double hrR = dataR[destOddStartIndex + r];
338                        double hrI = dataI[destOddStartIndex + r];
339    
340                        // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr
341                        dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI * hrI;
342                        dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI * hrR;
343                        // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr
344                        dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI * hrI);
345                        dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI * hrR);
346    
347                        // WsubN0ToR *= WsubN0R
348                        double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I;
349                        double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R;
350                        wSubN0ToRR = nextWsubN0ToRR;
351                        wSubN0ToRI = nextWsubN0ToRI;
352                    }
353                }
354    
355                lastN0 = n0;
356                lastLogN0 = logN0;
357            }
358    
359            normalizeTransformedData(dataRI, normalization, type);
360        }
361    
362        /**
363         * Returns the (forward, inverse) transform of the specified real data set.
364         *
365         * @param f the real data array to be transformed
366         * @param type the type of transform (forward, inverse) to be performed
367         * @return the complex transformed array
368         * @throws MathIllegalArgumentException if the length of the data array is not a power of two
369         */
370        public Complex[] transform(final double[] f, final TransformType type) {
371            final double[][] dataRI = new double[][] {
372                MathArrays.copyOf(f, f.length), new double[f.length]
373            };
374    
375            transformInPlace(dataRI, normalization, type);
376    
377            return TransformUtils.createComplexArray(dataRI);
378        }
379    
380        /**
381         * Returns the (forward, inverse) transform of the specified real function,
382         * sampled on the specified interval.
383         *
384         * @param f the function to be sampled and transformed
385         * @param min the (inclusive) lower bound for the interval
386         * @param max the (exclusive) upper bound for the interval
387         * @param n the number of sample points
388         * @param type the type of transform (forward, inverse) to be performed
389         * @return the complex transformed array
390         * @throws org.apache.commons.math3.exception.NumberIsTooLargeException
391         *   if the lower bound is greater than, or equal to the upper bound
392         * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
393         *   if the number of sample points {@code n} is negative
394         * @throws MathIllegalArgumentException if the number of sample points
395         *   {@code n} is not a power of two
396         */
397        public Complex[] transform(final UnivariateFunction f,
398                                   final double min, final double max, final int n,
399                                   final TransformType type) {
400    
401            final double[] data = FunctionUtils.sample(f, min, max, n);
402            return transform(data, type);
403        }
404    
405        /**
406         * Returns the (forward, inverse) transform of the specified complex data set.
407         *
408         * @param f the complex data array to be transformed
409         * @param type the type of transform (forward, inverse) to be performed
410         * @return the complex transformed array
411         * @throws MathIllegalArgumentException if the length of the data array is not a power of two
412         */
413        public Complex[] transform(final Complex[] f, final TransformType type) {
414            final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);
415    
416            transformInPlace(dataRI, normalization, type);
417    
418            return TransformUtils.createComplexArray(dataRI);
419        }
420    
421        /**
422         * Performs a multi-dimensional Fourier transform on a given array. Use
423         * {@link #transform(Complex[], TransformType)} in a row-column
424         * implementation in any number of dimensions with
425         * O(N&times;log(N)) complexity with
426         * N = n<sub>1</sub> &times; n<sub>2</sub> &times;n<sub>3</sub> &times; ...
427         * &times; n<sub>d</sub>, where n<sub>k</sub> is the number of elements in
428         * dimension k, and d is the total number of dimensions.
429         *
430         * @param mdca Multi-Dimensional Complex Array, i.e. {@code Complex[][][][]}
431         * @param type the type of transform (forward, inverse) to be performed
432         * @return transform of {@code mdca} as a Multi-Dimensional Complex Array, i.e. {@code Complex[][][][]}
433         * @throws IllegalArgumentException if any dimension is not a power of two
434         * @deprecated see MATH-736
435         */
436        @Deprecated
437        public Object mdfft(Object mdca, TransformType type) {
438            MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)
439                    new MultiDimensionalComplexMatrix(mdca).clone();
440            int[] dimensionSize = mdcm.getDimensionSizes();
441            //cycle through each dimension
442            for (int i = 0; i < dimensionSize.length; i++) {
443                mdfft(mdcm, type, i, new int[0]);
444            }
445            return mdcm.getArray();
446        }
447    
448        /**
449         * Performs one dimension of a multi-dimensional Fourier transform.
450         *
451         * @param mdcm input matrix
452         * @param type the type of transform (forward, inverse) to be performed
453         * @param d index of the dimension to process
454         * @param subVector recursion subvector
455         * @throws IllegalArgumentException if any dimension is not a power of two
456         * @deprecated see MATH-736
457         */
458        @Deprecated
459        private void mdfft(MultiDimensionalComplexMatrix mdcm,
460                TransformType type, int d, int[] subVector) {
461    
462            int[] dimensionSize = mdcm.getDimensionSizes();
463            //if done
464            if (subVector.length == dimensionSize.length) {
465                Complex[] temp = new Complex[dimensionSize[d]];
466                for (int i = 0; i < dimensionSize[d]; i++) {
467                    //fft along dimension d
468                    subVector[d] = i;
469                    temp[i] = mdcm.get(subVector);
470                }
471    
472                temp = transform(temp, type);
473    
474                for (int i = 0; i < dimensionSize[d]; i++) {
475                    subVector[d] = i;
476                    mdcm.set(temp[i], subVector);
477                }
478            } else {
479                int[] vector = new int[subVector.length + 1];
480                System.arraycopy(subVector, 0, vector, 0, subVector.length);
481                if (subVector.length == d) {
482                    //value is not important once the recursion is done.
483                    //then an fft will be applied along the dimension d.
484                    vector[d] = 0;
485                    mdfft(mdcm, type, d, vector);
486                } else {
487                    for (int i = 0; i < dimensionSize[subVector.length]; i++) {
488                        vector[subVector.length] = i;
489                        //further split along the next dimension
490                        mdfft(mdcm, type, d, vector);
491                    }
492                }
493            }
494        }
495    
496        /**
497         * Complex matrix implementation. Not designed for synchronized access may
498         * eventually be replaced by jsr-83 of the java community process
499         * http://jcp.org/en/jsr/detail?id=83
500         * may require additional exception throws for other basic requirements.
501         *
502         * @deprecated see MATH-736
503         */
504        @Deprecated
505        private static class MultiDimensionalComplexMatrix
506            implements Cloneable {
507    
508            /** Size in all dimensions. */
509            protected int[] dimensionSize;
510    
511            /** Storage array. */
512            protected Object multiDimensionalComplexArray;
513    
514            /**
515             * Simple constructor.
516             *
517             * @param multiDimensionalComplexArray array containing the matrix
518             * elements
519             */
520            public MultiDimensionalComplexMatrix(
521                    Object multiDimensionalComplexArray) {
522    
523                this.multiDimensionalComplexArray = multiDimensionalComplexArray;
524    
525                // count dimensions
526                int numOfDimensions = 0;
527                for (Object lastDimension = multiDimensionalComplexArray;
528                     lastDimension instanceof Object[];) {
529                    final Object[] array = (Object[]) lastDimension;
530                    numOfDimensions++;
531                    lastDimension = array[0];
532                }
533    
534                // allocate array with exact count
535                dimensionSize = new int[numOfDimensions];
536    
537                // fill array
538                numOfDimensions = 0;
539                for (Object lastDimension = multiDimensionalComplexArray;
540                     lastDimension instanceof Object[];) {
541                    final Object[] array = (Object[]) lastDimension;
542                    dimensionSize[numOfDimensions++] = array.length;
543                    lastDimension = array[0];
544                }
545    
546            }
547    
548            /**
549             * Get a matrix element.
550             *
551             * @param vector indices of the element
552             * @return matrix element
553             * @exception DimensionMismatchException if dimensions do not match
554             */
555            public Complex get(int... vector)
556                    throws DimensionMismatchException {
557    
558                if (vector == null) {
559                    if (dimensionSize.length > 0) {
560                        throw new DimensionMismatchException(
561                                0,
562                                dimensionSize.length);
563                    }
564                    return null;
565                }
566                if (vector.length != dimensionSize.length) {
567                    throw new DimensionMismatchException(
568                            vector.length,
569                            dimensionSize.length);
570                }
571    
572                Object lastDimension = multiDimensionalComplexArray;
573    
574                for (int i = 0; i < dimensionSize.length; i++) {
575                    lastDimension = ((Object[]) lastDimension)[vector[i]];
576                }
577                return (Complex) lastDimension;
578            }
579    
580            /**
581             * Set a matrix element.
582             *
583             * @param magnitude magnitude of the element
584             * @param vector indices of the element
585             * @return the previous value
586             * @exception DimensionMismatchException if dimensions do not match
587             */
588            public Complex set(Complex magnitude, int... vector)
589                    throws DimensionMismatchException {
590    
591                if (vector == null) {
592                    if (dimensionSize.length > 0) {
593                        throw new DimensionMismatchException(
594                                0,
595                                dimensionSize.length);
596                    }
597                    return null;
598                }
599                if (vector.length != dimensionSize.length) {
600                    throw new DimensionMismatchException(
601                            vector.length,
602                            dimensionSize.length);
603                }
604    
605                Object[] lastDimension = (Object[]) multiDimensionalComplexArray;
606                for (int i = 0; i < dimensionSize.length - 1; i++) {
607                    lastDimension = (Object[]) lastDimension[vector[i]];
608                }
609    
610                Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]];
611                lastDimension[vector[dimensionSize.length - 1]] = magnitude;
612    
613                return lastValue;
614            }
615    
616            /**
617             * Get the size in all dimensions.
618             *
619             * @return size in all dimensions
620             */
621            public int[] getDimensionSizes() {
622                return dimensionSize.clone();
623            }
624    
625            /**
626             * Get the underlying storage array.
627             *
628             * @return underlying storage array
629             */
630            public Object getArray() {
631                return multiDimensionalComplexArray;
632            }
633    
634            /** {@inheritDoc} */
635            @Override
636            public Object clone() {
637                MultiDimensionalComplexMatrix mdcm =
638                        new MultiDimensionalComplexMatrix(Array.newInstance(
639                        Complex.class, dimensionSize));
640                clone(mdcm);
641                return mdcm;
642            }
643    
644            /**
645             * Copy contents of current array into mdcm.
646             *
647             * @param mdcm array where to copy data
648             */
649            private void clone(MultiDimensionalComplexMatrix mdcm) {
650    
651                int[] vector = new int[dimensionSize.length];
652                int size = 1;
653                for (int i = 0; i < dimensionSize.length; i++) {
654                    size *= dimensionSize[i];
655                }
656                int[][] vectorList = new int[size][dimensionSize.length];
657                for (int[] nextVector : vectorList) {
658                    System.arraycopy(vector, 0, nextVector, 0,
659                                     dimensionSize.length);
660                    for (int i = 0; i < dimensionSize.length; i++) {
661                        vector[i]++;
662                        if (vector[i] < dimensionSize[i]) {
663                            break;
664                        } else {
665                            vector[i] = 0;
666                        }
667                    }
668                }
669    
670                for (int[] nextVector : vectorList) {
671                    mdcm.set(get(nextVector), nextVector);
672                }
673            }
674        }
675    }