View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.transform;
18  
19  import java.io.Serializable;
20  import java.lang.reflect.Array;
21  
22  import org.apache.commons.math3.analysis.FunctionUtils;
23  import org.apache.commons.math3.analysis.UnivariateFunction;
24  import org.apache.commons.math3.complex.Complex;
25  import org.apache.commons.math3.exception.DimensionMismatchException;
26  import org.apache.commons.math3.exception.MathIllegalArgumentException;
27  import org.apache.commons.math3.exception.MathIllegalStateException;
28  import org.apache.commons.math3.exception.util.LocalizedFormats;
29  import org.apache.commons.math3.util.ArithmeticUtils;
30  import org.apache.commons.math3.util.FastMath;
31  import org.apache.commons.math3.util.MathArrays;
32  
33  /**
34   * Implements the Fast Fourier Transform for transformation of one-dimensional
35   * real or complex data sets. For reference, see <em>Applied Numerical Linear
36   * Algebra</em>, ISBN 0898713897, chapter 6.
37   * <p>
38   * There are several variants of the discrete Fourier transform, with various
39   * normalization conventions, which are specified by the parameter
40   * {@link DftNormalization}.
41   * <p>
42   * The current implementation of the discrete Fourier transform as a fast
43   * Fourier transform requires the length of the data set to be a power of 2.
44   * This greatly simplifies and speeds up the code. Users can pad the data with
45   * zeros to meet this requirement. There are other flavors of FFT, for
46   * reference, see S. Winograd,
47   * <i>On computing the discrete Fourier transform</i>, Mathematics of
48   * Computation, 32 (1978), 175 - 199.
49   *
50   * @see DftNormalization
51   * @version $Id: FastFourierTransformer.java 1385310 2012-09-16 16:32:10Z tn $
52   * @since 1.2
53   */
54  public class FastFourierTransformer implements Serializable {
55  
56      /** Serializable version identifier. */
57      static final long serialVersionUID = 20120210L;
58  
59      /**
60       * {@code W_SUB_N_R[i]} is the real part of
61       * {@code exp(- 2 * i * pi / n)}:
62       * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}.
63       */
64      private static final double[] W_SUB_N_R =
65              {  0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1
66              , 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1
67              , 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1
68              , 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1
69              , 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1
70              , 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1
71              , 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1
72              , 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0
73              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
74              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
75              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
76              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
77              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
78              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
79              , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
80              , 0x1.0p0, 0x1.0p0, 0x1.0p0 };
81  
82      /**
83       * {@code W_SUB_N_I[i]} is the imaginary part of
84       * {@code exp(- 2 * i * pi / n)}:
85       * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}.
86       */
87      private static final double[] W_SUB_N_I =
88              {  0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1
89              , -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5
90              , -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9
91              , -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13
92              , -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17
93              , -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21
94              , -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25
95              , -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29
96              , -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33
97              , -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37
98              , -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41
99              , -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 }