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 */
017package org.apache.commons.math3.transform;
018
019import java.io.Serializable;
020import java.lang.reflect.Array;
021
022import org.apache.commons.math3.analysis.FunctionUtils;
023import org.apache.commons.math3.analysis.UnivariateFunction;
024import org.apache.commons.math3.complex.Complex;
025import org.apache.commons.math3.exception.DimensionMismatchException;
026import org.apache.commons.math3.exception.MathIllegalArgumentException;
027import org.apache.commons.math3.exception.MathIllegalStateException;
028import org.apache.commons.math3.exception.util.LocalizedFormats;
029import org.apache.commons.math3.util.ArithmeticUtils;
030import org.apache.commons.math3.util.FastMath;
031import 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 */
054public 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}