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×log(N)) complexity with
426 * N = n<sub>1</sub> × n<sub>2</sub> ×n<sub>3</sub> × ...
427 * × 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 }