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