FastFourierTransformer.java

  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;

  19. import java.io.Serializable;
  20. import java.lang.reflect.Array;

  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;

  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.  * @since 1.2
  52.  */
  53. public class FastFourierTransformer implements Serializable {

  55.     /** Serializable version identifier. */
  56.     static final long serialVersionUID = 20120210L;

  58.     /**
  59.      * {@code W_SUB_N_R[i]} is the real part of
  60.      * {@code exp(- 2 * i * pi / n)}:
  61.      * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}.
  62.      */
  63.     private static final double[] W_SUB_N_R =
  64.             {  0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1
  65.             , 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1
  66.             , 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1
  67.             , 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1
  68.             , 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1
  69.             , 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1
  70.             , 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1
  71.             , 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0
  72.             , 0x1.0p0, 0x1.0p0, 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 };

  81.     /**
  82.      * {@code W_SUB_N_I[i]} is the imaginary part of
  83.      * {@code exp(- 2 * i * pi / n)}:
  84.      * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}.
  85.      */
  86.     private static final double[] W_SUB_N_I =
  87.             {  0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1
  88.             , -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5
  89.             , -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9
  90.             , -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13
  91.             , -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17
  92.             , -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21
  93.             , -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25
  94.             , -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29
  95.             , -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33
  96.             , -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37
  97.             , -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41
  98.             , -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45
  99.             , -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49
  100.             , -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53
  101.             , -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57
  102.             , -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 };

  104.     /** The type of DFT to be performed. */
  105.     private final DftNormalization normalization;

  107.     /**
  108.      * Creates a new instance of this class, with various normalization
  109.      * conventions.
  110.      *
  111.      * @param normalization the type of normalization to be applied to the
  112.      * transformed data
  113.      */
  114.     public FastFourierTransformer(final DftNormalization normalization) {
  115.         this.normalization = normalization;
  116.     }

  118.     /**
  119.      * Performs identical index bit reversal shuffles on two arrays of identical
  120.      * size. Each element in the array is swapped with another element based on
  121.      * the bit-reversal of the index. For example, in an array with length 16,
  122.      * item at binary index 0011 (decimal 3) would be swapped with the item at
  123.      * binary index 1100 (decimal 12).
  124.      *
  125.      * @param a the first array to be shuffled
  126.      * @param b the second array to be shuffled
  127.      */
  128.     private static void bitReversalShuffle2(double[] a, double[] b) {
  129.         final int n = a.length;
  130.         assert b.length == n;
  131.         final int halfOfN = n >> 1;

  133.         int j = 0;
  134.         for (int i = 0; i < n; i++) {
  135.             if (i < j) {
  136.                 // swap indices i & j
  137.                 double temp = a[i];
  138.                 a[i] = a[j];
  139.                 a[j] = temp;

  141.                 temp = b[i];
  142.                 b[i] = b[j];
  143.                 b[j] = temp;
  144.             }

  146.             int k = halfOfN;
  147.             while (k <= j && k > 0) {
  148.                 j -= k;
  149.                 k >>= 1;
  150.             }
  151.             j += k;
  152.         }
  153.     }

  155.     /**
  156.      * Applies the proper normalization to the specified transformed data.
  157.      *
  158.      * @param dataRI the unscaled transformed data
  159.      * @param normalization the normalization to be applied
  160.      * @param type the type of transform (forward, inverse) which resulted in the specified data
  161.      */
  162.     private static void normalizeTransformedData(final double[][] dataRI,
  163.         final DftNormalization normalization, final TransformType type) {

  165.         final double[] dataR = dataRI[0];
  166.         final double[] dataI = dataRI[1];
  167.         final int n = dataR.length;
  168.         assert dataI.length == n;

  170.         switch (normalization) {
  171.             case STANDARD:
  172.                 if (type == TransformType.INVERSE) {
  173.                     final double scaleFactor = 1.0 / ((double) n);
  174.                     for (int i = 0; i < n; i++) {
  175.                         dataR[i] *= scaleFactor;
  176.                         dataI[i] *= scaleFactor;
  177.                     }
  178.                 }
  179.                 break;
  180.             case UNITARY:
  181.                 final double scaleFactor = 1.0 / FastMath.sqrt(n);
  182.                 for (int i = 0; i < n; i++) {
  183.                     dataR[i] *= scaleFactor;
  184.                     dataI[i] *= scaleFactor;
  185.                 }
  186.                 break;
  187.             default:
  188.                 /*
  189.                  * This should never occur in normal conditions. However this
  190.                  * clause has been added as a safeguard if other types of
  191.                  * normalizations are ever implemented, and the corresponding
  192.                  * test is forgotten in the present switch.
  193.                  */
  194.                 throw new MathIllegalStateException();
  195.         }
  196.     }

  198.     /**
  199.      * Computes the standard transform of the specified complex data. The
  200.      * computation is done in place. The input data is laid out as follows
  201.      * <ul>
  202.      *   <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point,</li>
  203.      *   <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data point.</li>
  204.      * </ul>
  205.      *
  206.      * @param dataRI the two dimensional array of real and imaginary parts of the data
  207.      * @param normalization the normalization to be applied to the transformed data
  208.      * @param type the type of transform (forward, inverse) to be performed
  209.      * @throws DimensionMismatchException if the number of rows of the specified
  210.      *   array is not two, or the array is not rectangular
  211.      * @throws MathIllegalArgumentException if the number of data points is not
  212.      *   a power of two
  213.      */
  214.     public static void transformInPlace(final double[][] dataRI,
  215.         final DftNormalization normalization, final TransformType type) {

  217.         if (dataRI.length != 2) {
  218.             throw new DimensionMismatchException(dataRI.length, 2);
  219.         }
  220.         final double[] dataR = dataRI[0];
  221.         final double[] dataI = dataRI[1];
  222.         if (dataR.length != dataI.length) {
  223.             throw new DimensionMismatchException(dataI.length, dataR.length);
  224.         }

  226.         final int n = dataR.length;
  227.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  228.             throw new MathIllegalArgumentException(
  229.                 LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
  230.                 Integer.valueOf(n));
  231.         }

  233.         if (n == 1) {
  234.             return;
  235.         } else if (n == 2) {
  236.             final double srcR0 = dataR[0];
  237.             final double srcI0 = dataI[0];
  238.             final double srcR1 = dataR[1];
  239.             final double srcI1 = dataI[1];

  241.             // X_0 = x_0 + x_1
  242.             dataR[0] = srcR0 + srcR1;
  243.             dataI[0] = srcI0 + srcI1;
  244.             // X_1 = x_0 - x_1
  245.             dataR[1] = srcR0 - srcR1;
  246.             dataI[1] = srcI0 - srcI1;

  248.             normalizeTransformedData(dataRI, normalization, type);
  249.             return;
  250.         }

  252.         bitReversalShuffle2(dataR, dataI);

  254.         // Do 4-term DFT.
  255.         if (type == TransformType.INVERSE) {
  256.             for (int i0 = 0; i0 < n; i0 += 4) {
  257.                 final int i1 = i0 + 1;
  258.                 final int i2 = i0 + 2;
  259.                 final int i3 = i0 + 3;

  261.                 final double srcR0 = dataR[i0];
  262.                 final double srcI0 = dataI[i0];
  263.                 final double srcR1 = dataR[i2];
  264.                 final double srcI1 = dataI[i2];
  265.                 final double srcR2 = dataR[i1];
  266.                 final double srcI2 = dataI[i1];
  267.                 final double srcR3 = dataR[i3];
  268.                 final double srcI3 = dataI[i3];

  270.                 // 4-term DFT
  271.                 // X_0 = x_0 + x_1 + x_2 + x_3
  272.                 dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
  273.                 dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
  274.                 // X_1 = x_0 - x_2 + j * (x_3 - x_1)
  275.                 dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1);
  276.                 dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3);
  277.                 // X_2 = x_0 - x_1 + x_2 - x_3
  278.                 dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
  279.                 dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
  280.                 // X_3 = x_0 - x_2 + j * (x_1 - x_3)
  281.                 dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3);
  282.                 dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1);
  283.             }
  284.         } else {
  285.             for (int i0 = 0; i0 < n; i0 += 4) {
  286.                 final int i1 = i0 + 1;
  287.                 final int i2 = i0 + 2;
  288.                 final int i3 = i0 + 3;

  290.                 final double srcR0 = dataR[i0];
  291.                 final double srcI0 = dataI[i0];
  292.                 final double srcR1 = dataR[i2];
  293.                 final double srcI1 = dataI[i2];
  294.                 final double srcR2 = dataR[i1];
  295.                 final double srcI2 = dataI[i1];
  296.                 final double srcR3 = dataR[i3];
  297.                 final double srcI3 = dataI[i3];

  299.                 // 4-term DFT
  300.                 // X_0 = x_0 + x_1 + x_2 + x_3
  301.                 dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
  302.                 dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
  303.                 // X_1 = x_0 - x_2 + j * (x_3 - x_1)
  304.                 dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3);
  305.                 dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1);
  306.                 // X_2 = x_0 - x_1 + x_2 - x_3
  307.                 dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
  308.                 dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
  309.                 // X_3 = x_0 - x_2 + j * (x_1 - x_3)
  310.                 dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1);
  311.                 dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3);
  312.             }
  313.         }

  315.         int lastN0 = 4;
  316.         int lastLogN0 = 2;
  317.         while (lastN0 < n) {
  318.             int n0 = lastN0 << 1;
  319.             int logN0 = lastLogN0 + 1;
  320.             double wSubN0R = W_SUB_N_R[logN0];
  321.             double wSubN0I = W_SUB_N_I[logN0];
  322.             if (type == TransformType.INVERSE) {
  323.                 wSubN0I = -wSubN0I;
  324.             }

  326.             // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2).
  327.             for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) {
  328.                 int destOddStartIndex = destEvenStartIndex + lastN0;

  330.                 double wSubN0ToRR = 1;
  331.                 double wSubN0ToRI = 0;

  333.                 for (int r = 0; r < lastN0; r++) {
  334.                     double grR = dataR[destEvenStartIndex + r];
  335.                     double grI = dataI[destEvenStartIndex + r];
  336.                     double hrR = dataR[destOddStartIndex + r];
  337.                     double hrI = dataI[destOddStartIndex + r];

  339.                     // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr
  340.                     dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI * hrI;
  341.                     dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI * hrR;
  342.                     // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr
  343.                     dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI * hrI);
  344.                     dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI * hrR);

  346.                     // WsubN0ToR *= WsubN0R
  347.                     double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I;
  348.                     double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R;
  349.                     wSubN0ToRR = nextWsubN0ToRR;
  350.                     wSubN0ToRI = nextWsubN0ToRI;
  351.                 }
  352.             }

  354.             lastN0 = n0;
  355.             lastLogN0 = logN0;
  356.         }

  358.         normalizeTransformedData(dataRI, normalization, type);
  359.     }

  361.     /**
  362.      * Returns the (forward, inverse) transform of the specified real data set.
  363.      *
  364.      * @param f the real data array to be transformed
  365.      * @param type the type of transform (forward, inverse) to be performed
  366.      * @return the complex transformed array
  367.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  368.      */
  369.     public Complex[] transform(final double[] f, final TransformType type) {
  370.         final double[][] dataRI = new double[][] {
  371.             MathArrays.copyOf(f, f.length), new double[f.length]
  372.         };

  374.         transformInPlace(dataRI, normalization, type);

  376.         return TransformUtils.createComplexArray(dataRI);
  377.     }

  379.     /**
  380.      * Returns the (forward, inverse) transform of the specified real function,
  381.      * sampled on the specified interval.
  382.      *
  383.      * @param f the function to be sampled and transformed
  384.      * @param min the (inclusive) lower bound for the interval
  385.      * @param max the (exclusive) upper bound for the interval
  386.      * @param n the number of sample points
  387.      * @param type the type of transform (forward, inverse) to be performed
  388.      * @return the complex transformed array
  389.      * @throws org.apache.commons.math3.exception.NumberIsTooLargeException
  390.      *   if the lower bound is greater than, or equal to the upper bound
  391.      * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
  392.      *   if the number of sample points {@code n} is negative
  393.      * @throws MathIllegalArgumentException if the number of sample points
  394.      *   {@code n} is not a power of two
  395.      */
  396.     public Complex[] transform(final UnivariateFunction f,
  397.                                final double min, final double max, final int n,
  398.                                final TransformType type) {

  400.         final double[] data = FunctionUtils.sample(f, min, max, n);
  401.         return transform(data, type);
  402.     }

  404.     /**
  405.      * Returns the (forward, inverse) transform of the specified complex data set.
  406.      *
  407.      * @param f the complex data array to be transformed
  408.      * @param type the type of transform (forward, inverse) to be performed
  409.      * @return the complex transformed array
  410.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  411.      */
  412.     public Complex[] transform(final Complex[] f, final TransformType type) {
  413.         final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);

  415.         transformInPlace(dataRI, normalization, type);

  417.         return TransformUtils.createComplexArray(dataRI);
  418.     }

  420.     /**
  421.      * Performs a multi-dimensional Fourier transform on a given array. Use
  422.      * {@link #transform(Complex[], TransformType)} in a row-column
  423.      * implementation in any number of dimensions with
  424.      * O(N&times;log(N)) complexity with
  425.      * N = n<sub>1</sub> &times; n<sub>2</sub> &times;n<sub>3</sub> &times; ...
  426.      * &times; n<sub>d</sub>, where n<sub>k</sub> is the number of elements in
  427.      * dimension k, and d is the total number of dimensions.
  428.      *
  429.      * @param mdca Multi-Dimensional Complex Array, i.e. {@code Complex[][][][]}
  430.      * @param type the type of transform (forward, inverse) to be performed
  431.      * @return transform of {@code mdca} as a Multi-Dimensional Complex Array, i.e. {@code Complex[][][][]}
  432.      * @throws IllegalArgumentException if any dimension is not a power of two
  433.      * @deprecated see MATH-736
  434.      */
  435.     @Deprecated
  436.     public Object mdfft(Object mdca, TransformType type) {
  437.         MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)
  438.                 new MultiDimensionalComplexMatrix(mdca).clone();
  439.         int[] dimensionSize = mdcm.getDimensionSizes();
  440.         //cycle through each dimension
  441.         for (int i = 0; i < dimensionSize.length; i++) {
  442.             mdfft(mdcm, type, i, new int[0]);
  443.         }
  444.         return mdcm.getArray();
  445.     }

  447.     /**
  448.      * Performs one dimension of a multi-dimensional Fourier transform.
  449.      *
  450.      * @param mdcm input matrix
  451.      * @param type the type of transform (forward, inverse) to be performed
  452.      * @param d index of the dimension to process
  453.      * @param subVector recursion subvector
  454.      * @throws IllegalArgumentException if any dimension is not a power of two
  455.      * @deprecated see MATH-736
  456.      */
  457.     @Deprecated
  458.     private void mdfft(MultiDimensionalComplexMatrix mdcm,
  459.             TransformType type, int d, int[] subVector) {

  461.         int[] dimensionSize = mdcm.getDimensionSizes();
  462.         //if done
  463.         if (subVector.length == dimensionSize.length) {
  464.             Complex[] temp = new Complex[dimensionSize[d]];
  465.             for (int i = 0; i < dimensionSize[d]; i++) {
  466.                 //fft along dimension d
  467.                 subVector[d] = i;
  468.                 temp[i] = mdcm.get(subVector);
  469.             }

  471.             temp = transform(temp, type);

  473.             for (int i = 0; i < dimensionSize[d]; i++) {
  474.                 subVector[d] = i;
  475.                 mdcm.set(temp[i], subVector);
  476.             }
  477.         } else {
  478.             int[] vector = new int[subVector.length + 1];
  479.             System.arraycopy(subVector, 0, vector, 0, subVector.length);
  480.             if (subVector.length == d) {
  481.                 //value is not important once the recursion is done.
  482.                 //then an fft will be applied along the dimension d.
  483.                 vector[d] = 0;
  484.                 mdfft(mdcm, type, d, vector);
  485.             } else {
  486.                 for (int i = 0; i < dimensionSize[subVector.length]; i++) {
  487.                     vector[subVector.length] = i;
  488.                     //further split along the next dimension
  489.                     mdfft(mdcm, type, d, vector);
  490.                 }
  491.             }
  492.         }
  493.     }

  495.     /**
  496.      * Complex matrix implementation. Not designed for synchronized access may
  497.      * eventually be replaced by jsr-83 of the java community process
  498.      * http://jcp.org/en/jsr/detail?id=83
  499.      * may require additional exception throws for other basic requirements.
  500.      *
  501.      * @deprecated see MATH-736
  502.      */
  503.     @Deprecated
  504.     private static class MultiDimensionalComplexMatrix
  505.         implements Cloneable {

  507.         /** Size in all dimensions. */
  508.         protected int[] dimensionSize;

  510.         /** Storage array. */
  511.         protected Object multiDimensionalComplexArray;

  513.         /**
  514.          * Simple constructor.
  515.          *
  516.          * @param multiDimensionalComplexArray array containing the matrix
  517.          * elements
  518.          */
  519.         public MultiDimensionalComplexMatrix(
  520.                 Object multiDimensionalComplexArray) {

  522.             this.multiDimensionalComplexArray = multiDimensionalComplexArray;

  524.             // count dimensions
  525.             int numOfDimensions = 0;
  526.             for (Object lastDimension = multiDimensionalComplexArray;
  527.                  lastDimension instanceof Object[];) {
  528.                 final Object[] array = (Object[]) lastDimension;
  529.                 numOfDimensions++;
  530.                 lastDimension = array[0];
  531.             }

  533.             // allocate array with exact count
  534.             dimensionSize = new int[numOfDimensions];

  536.             // fill array
  537.             numOfDimensions = 0;
  538.             for (Object lastDimension = multiDimensionalComplexArray;
  539.                  lastDimension instanceof Object[];) {
  540.                 final Object[] array = (Object[]) lastDimension;
  541.                 dimensionSize[numOfDimensions++] = array.length;
  542.                 lastDimension = array[0];
  543.             }

  545.         }

  547.         /**
  548.          * Get a matrix element.
  549.          *
  550.          * @param vector indices of the element
  551.          * @return matrix element
  552.          * @exception DimensionMismatchException if dimensions do not match
  553.          */
  554.         public Complex get(int... vector)
  555.                 throws DimensionMismatchException {

  557.             if (vector == null) {
  558.                 if (dimensionSize.length > 0) {
  559.                     throw new DimensionMismatchException(
  560.                             0,
  561.                             dimensionSize.length);
  562.                 }
  563.                 return null;
  564.             }
  565.             if (vector.length != dimensionSize.length) {
  566.                 throw new DimensionMismatchException(
  567.                         vector.length,
  568.                         dimensionSize.length);
  569.             }

  571.             Object lastDimension = multiDimensionalComplexArray;

  573.             for (int i = 0; i < dimensionSize.length; i++) {
  574.                 lastDimension = ((Object[]) lastDimension)[vector[i]];
  575.             }
  576.             return (Complex) lastDimension;
  577.         }

  579.         /**
  580.          * Set a matrix element.
  581.          *
  582.          * @param magnitude magnitude of the element
  583.          * @param vector indices of the element
  584.          * @return the previous value
  585.          * @exception DimensionMismatchException if dimensions do not match
  586.          */
  587.         public Complex set(Complex magnitude, int... vector)
  588.                 throws DimensionMismatchException {

  590.             if (vector == null) {
  591.                 if (dimensionSize.length > 0) {
  592.                     throw new DimensionMismatchException(
  593.                             0,
  594.                             dimensionSize.length);
  595.                 }
  596.                 return null;
  597.             }
  598.             if (vector.length != dimensionSize.length) {
  599.                 throw new DimensionMismatchException(
  600.                         vector.length,
  601.                         dimensionSize.length);
  602.             }

  604.             Object[] lastDimension = (Object[]) multiDimensionalComplexArray;
  605.             for (int i = 0; i < dimensionSize.length - 1; i++) {
  606.                 lastDimension = (Object[]) lastDimension[vector[i]];
  607.             }

  609.             Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]];
  610.             lastDimension[vector[dimensionSize.length - 1]] = magnitude;

  612.             return lastValue;
  613.         }

  615.         /**
  616.          * Get the size in all dimensions.
  617.          *
  618.          * @return size in all dimensions
  619.          */
  620.         public int[] getDimensionSizes() {
  621.             return dimensionSize.clone();
  622.         }

  624.         /**
  625.          * Get the underlying storage array.
  626.          *
  627.          * @return underlying storage array
  628.          */
  629.         public Object getArray() {
  630.             return multiDimensionalComplexArray;
  631.         }

  633.         /** {@inheritDoc} */
  634.         @Override
  635.         public Object clone() {
  636.             MultiDimensionalComplexMatrix mdcm =
  637.                     new MultiDimensionalComplexMatrix(Array.newInstance(
  638.                     Complex.class, dimensionSize));
  639.             clone(mdcm);
  640.             return mdcm;
  641.         }

  643.         /**
  644.          * Copy contents of current array into mdcm.
  645.          *
  646.          * @param mdcm array where to copy data
  647.          */
  648.         private void clone(MultiDimensionalComplexMatrix mdcm) {

  650.             int[] vector = new int[dimensionSize.length];
  651.             int size = 1;
  652.             for (int i = 0; i < dimensionSize.length; i++) {
  653.                 size *= dimensionSize[i];
  654.             }
  655.             int[][] vectorList = new int[size][dimensionSize.length];
  656.             for (int[] nextVector : vectorList) {
  657.                 System.arraycopy(vector, 0, nextVector, 0,
  658.                                  dimensionSize.length);
  659.                 for (int i = 0; i < dimensionSize.length; i++) {
  660.                     vector[i]++;
  661.                     if (vector[i] < dimensionSize[i]) {
  662.                         break;
  663.                     } else {
  664.                         vector[i] = 0;
  665.                     }
  666.                 }
  667.             }

  669.             for (int[] nextVector : vectorList) {
  670.                 mdcm.set(get(nextVector), nextVector);
  671.             }
  672.         }
  673.     }
  674. }
Created with JaCoCo 0.7.2.201409121644