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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.transform;

  24. import java.io.Serializable;

  26. import org.hipparchus.analysis.FunctionUtils;
  27. import org.hipparchus.analysis.UnivariateFunction;
  28. import org.hipparchus.complex.Complex;
  29. import org.hipparchus.exception.MathIllegalArgumentException;
  30. import org.hipparchus.exception.MathRuntimeException;
  31. import org.hipparchus.util.ArithmeticUtils;
  32. import org.hipparchus.util.FastMath;
  33. import org.hipparchus.util.MathArrays;
  34. import org.hipparchus.util.MathUtils;

  36. /**
  37.  * Implements the Fast Fourier Transform for transformation of one-dimensional
  38.  * real or complex data sets. For reference, see <em>Applied Numerical Linear
  39.  * Algebra</em>, ISBN 0898713897, chapter 6.
  40.  * <p>
  41.  * There are several variants of the discrete Fourier transform, with various
  42.  * normalization conventions, which are specified by the parameter
  43.  * {@link DftNormalization}.
  44.  * <p>
  45.  * The current implementation of the discrete Fourier transform as a fast
  46.  * Fourier transform requires the length of the data set to be a power of 2.
  47.  * This greatly simplifies and speeds up the code. Users can pad the data with
  48.  * zeros to meet this requirement. There are other flavors of FFT, for
  49.  * reference, see S. Winograd,
  50.  * <i>On computing the discrete Fourier transform</i>, Mathematics of
  51.  * Computation, 32 (1978), 175 - 199.
  52.  *
  53.  * @see DftNormalization
  54.  */
  55. public class FastFourierTransformer implements Serializable {

  57.     /** Serializable version identifier. */
  58.     private static final long serialVersionUID = 20120210L;

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

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

  106.     /** The type of DFT to be performed. */
  107.     private final DftNormalization normalization;

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

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

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

  143.                 temp = b[i];
  144.                 b[i] = b[j];
  145.                 b[j] = temp;
  146.             }

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

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

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

  172.         switch (normalization) {
  173.             case STANDARD:
  174.                 if (type == TransformType.INVERSE) {
  175.                     final double scaleFactor = 1.0 / n;
  176.                     for (int i = 0; i < n; i++) {
  177.                         dataR[i] *= scaleFactor;
  178.                         dataI[i] *= scaleFactor;
  179.                     }
  180.                 }
  181.                 break;
  182.             case UNITARY:
  183.                 final double scaleFactor = 1.0 / FastMath.sqrt(n);
  184.                 for (int i = 0; i < n; i++) {
  185.                     dataR[i] *= scaleFactor;
  186.                     dataI[i] *= scaleFactor;
  187.                 }
  188.                 break;
  189.             default:
  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.                 throw MathRuntimeException.createInternalError();
  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 MathIllegalArgumentException 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.         MathUtils.checkDimension(dataRI.length, 2);
  218.         final double[] dataR = dataRI[0];
  219.         final double[] dataI = dataRI[1];
  220.         MathArrays.checkEqualLength(dataR, dataI);

  222.         final int n = dataR.length;
  223.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  224.             throw new MathIllegalArgumentException(LocalizedFFTFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
  225.                     n);
  226.         }

  228.         if (n == 1) {
  229.             return;
  230.         } else if (n == 2) {
  231.             final double srcR0 = dataR[0];
  232.             final double srcI0 = dataI[0];
  233.             final double srcR1 = dataR[1];
  234.             final double srcI1 = dataI[1];

  236.             // X_0 = x_0 + x_1
  237.             dataR[0] = srcR0 + srcR1;
  238.             dataI[0] = srcI0 + srcI1;
  239.             // X_1 = x_0 - x_1
  240.             dataR[1] = srcR0 - srcR1;
  241.             dataI[1] = srcI0 - srcI1;

  243.             normalizeTransformedData(dataRI, normalization, type);
  244.             return;
  245.         }

  247.         bitReversalShuffle2(dataR, dataI);

  249.         // Do 4-term DFT.
  250.         if (type == TransformType.INVERSE) {
  251.             for (int i0 = 0; i0 < n; i0 += 4) {
  252.                 final int i1 = i0 + 1;
  253.                 final int i2 = i0 + 2;
  254.                 final int i3 = i0 + 3;

  256.                 final double srcR0 = dataR[i0];
  257.                 final double srcI0 = dataI[i0];
  258.                 final double srcR1 = dataR[i2];
  259.                 final double srcI1 = dataI[i2];
  260.                 final double srcR2 = dataR[i1];
  261.                 final double srcI2 = dataI[i1];
  262.                 final double srcR3 = dataR[i3];
  263.                 final double srcI3 = dataI[i3];

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

  285.                 final double srcR0 = dataR[i0];
  286.                 final double srcI0 = dataI[i0];
  287.                 final double srcR1 = dataR[i2];
  288.                 final double srcI1 = dataI[i2];
  289.                 final double srcR2 = dataR[i1];
  290.                 final double srcI2 = dataI[i1];
  291.                 final double srcR3 = dataR[i3];
  292.                 final double srcI3 = dataI[i3];

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

  310.         int lastN0 = 4;
  311.         int lastLogN0 = 2;
  312.         while (lastN0 < n) {
  313.             int n0 = lastN0 << 1;
  314.             int logN0 = lastLogN0 + 1;
  315.             double wSubN0R = W_SUB_N_R[logN0];
  316.             double wSubN0I = W_SUB_N_I[logN0];
  317.             if (type == TransformType.INVERSE) {
  318.                 wSubN0I = -wSubN0I;
  319.             }

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

  325.                 double wSubN0ToRR = 1;
  326.                 double wSubN0ToRI = 0;

  328.                 for (int r = 0; r < lastN0; r++) {
  329.                     double grR = dataR[destEvenStartIndex + r];
  330.                     double grI = dataI[destEvenStartIndex + r];
  331.                     double hrR = dataR[destOddStartIndex + r];
  332.                     double hrI = dataI[destOddStartIndex + r];

  334.                     // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr
  335.                     dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI * hrI;
  336.                     dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI * hrR;
  337.                     // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr
  338.                     dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI * hrI);
  339.                     dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI * hrR);

  341.                     // WsubN0ToR *= WsubN0R
  342.                     double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I;
  343.                     double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R;
  344.                     wSubN0ToRR = nextWsubN0ToRR;
  345.                     wSubN0ToRI = nextWsubN0ToRI;
  346.                 }
  347.             }

  349.             lastN0 = n0;
  350.             lastLogN0 = logN0;
  351.         }

  353.         normalizeTransformedData(dataRI, normalization, type);
  354.     }

  356.     /**
  357.      * Returns the (forward, inverse) transform of the specified real data set.
  358.      *
  359.      * @param f the real data array to be transformed
  360.      * @param type the type of transform (forward, inverse) to be performed
  361.      * @return the complex transformed array
  362.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  363.      */
  364.     public Complex[] transform(final double[] f, final TransformType type) {
  365.         final double[][] dataRI = { f.clone(), new double[f.length] };
  366.         transformInPlace(dataRI, normalization, type);
  367.         return TransformUtils.createComplexArray(dataRI);
  368.     }

  370.     /**
  371.      * Returns the (forward, inverse) transform of the specified real function,
  372.      * sampled on the specified interval.
  373.      *
  374.      * @param f the function to be sampled and transformed
  375.      * @param min the (inclusive) lower bound for the interval
  376.      * @param max the (exclusive) upper bound for the interval
  377.      * @param n the number of sample points
  378.      * @param type the type of transform (forward, inverse) to be performed
  379.      * @return the complex transformed array
  380.      * @throws org.hipparchus.exception.MathIllegalArgumentException
  381.      *   if the lower bound is greater than, or equal to the upper bound
  382.      * @throws org.hipparchus.exception.MathIllegalArgumentException
  383.      *   if the number of sample points {@code n} is negative
  384.      * @throws MathIllegalArgumentException if the number of sample points
  385.      *   {@code n} is not a power of two
  386.      */
  387.     public Complex[] transform(final UnivariateFunction f,
  388.                                final double min, final double max, final int n,
  389.                                final TransformType type) {

  391.         final double[] data = FunctionUtils.sample(f, min, max, n);
  392.         return transform(data, type);
  393.     }

  395.     /**
  396.      * Returns the (forward, inverse) transform of the specified complex data set.
  397.      *
  398.      * @param f the complex data array to be transformed
  399.      * @param type the type of transform (forward, inverse) to be performed
  400.      * @return the complex transformed array
  401.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  402.      */
  403.     public Complex[] transform(final Complex[] f, final TransformType type) {
  404.         final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);

  406.         transformInPlace(dataRI, normalization, type);

  408.         return TransformUtils.createComplexArray(dataRI);
  409.     }

  411. }
Created with JaCoCo 0.8.12.202403310830