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 */
017package org.apache.commons.math4.transform;
018
019import java.io.Serializable;
020import org.apache.commons.math4.analysis.FunctionUtils;
021import org.apache.commons.math4.analysis.UnivariateFunction;
022import org.apache.commons.math4.complex.Complex;
023import org.apache.commons.math4.exception.DimensionMismatchException;
024import org.apache.commons.math4.exception.MathIllegalArgumentException;
025import org.apache.commons.math4.exception.MathIllegalStateException;
026import org.apache.commons.math4.exception.util.LocalizedFormats;
027import org.apache.commons.math4.util.ArithmeticUtils;
028import org.apache.commons.math4.util.FastMath;
029import org.apache.commons.math4.util.MathArrays;
030
031/**
032 * Implements the Fast Fourier Transform for transformation of one-dimensional
033 * real or complex data sets. For reference, see <em>Applied Numerical Linear
034 * Algebra</em>, ISBN 0898713897, chapter 6.
035 * <p>
036 * There are several variants of the discrete Fourier transform, with various
037 * normalization conventions, which are specified by the parameter
038 * {@link DftNormalization}.
039 * <p>
040 * The current implementation of the discrete Fourier transform as a fast
041 * Fourier transform requires the length of the data set to be a power of 2.
042 * This greatly simplifies and speeds up the code. Users can pad the data with
043 * zeros to meet this requirement. There are other flavors of FFT, for
044 * reference, see S. Winograd,
045 * <i>On computing the discrete Fourier transform</i>, Mathematics of
046 * Computation, 32 (1978), 175 - 199.
047 *
048 * @see DftNormalization
049 * @since 1.2
050 */
051public class FastFourierTransformer implements Serializable {
052    /** Serializable version identifier. */
053    static final long serialVersionUID = 20120210L;
054
055    /**
056     * {@code W_SUB_N_R[i]} is the real part of
057     * {@code exp(- 2 * i * pi / n)}:
058     * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}.
059     */
060    private static final double[] W_SUB_N_R =
061            {  0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1
062            , 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1
063            , 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1
064            , 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1
065            , 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1
066            , 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1
067            , 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1
068            , 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0
069            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
070            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
071            , 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0
072            , 0x1.0p0, 0x1.0p0, 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 };
077
078    /**
079     * {@code W_SUB_N_I[i]} is the imaginary part of
080     * {@code exp(- 2 * i * pi / n)}:
081     * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}.
082     */
083    private static final double[] W_SUB_N_I =
084            {  0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1
085            , -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5
086            , -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9
087            , -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13
088            , -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17
089            , -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21
090            , -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25
091            , -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29
092            , -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33
093            , -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37
094            , -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41
095            , -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45
096            , -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49
097            , -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53
098            , -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57
099            , -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 };
100
101    /** The type of DFT to be performed. */
102    private final DftNormalization normalization;
103
104    /**
105     * Creates a new instance of this class, with various normalization
106     * conventions.
107     *
108     * @param normalization the type of normalization to be applied to the
109     * transformed data
110     */
111    public FastFourierTransformer(final DftNormalization normalization) {
112        this.normalization = normalization;
113    }
114
115    /**
116     * Performs identical index bit reversal shuffles on two arrays of identical
117     * size. Each element in the array is swapped with another element based on
118     * the bit-reversal of the index. For example, in an array with length 16,
119     * item at binary index 0011 (decimal 3) would be swapped with the item at
120     * binary index 1100 (decimal 12).
121     *
122     * @param a the first array to be shuffled
123     * @param b the second array to be shuffled
124     */
125    private static void bitReversalShuffle2(double[] a, double[] b) {
126        final int n = a.length;
127        assert b.length == n;
128        final int halfOfN = n >> 1;
129
130        int j = 0;
131        for (int i = 0; i < n; i++) {
132            if (i < j) {
133                // swap indices i & j
134                double temp = a[i];
135                a[i] = a[j];
136                a[j] = temp;
137
138                temp = b[i];
139                b[i] = b[j];
140                b[j] = temp;
141            }
142
143            int k = halfOfN;
144            while (k <= j && k > 0) {
145                j -= k;
146                k >>= 1;
147            }
148            j += k;
149        }
150    }
151
152    /**
153     * Applies the proper normalization to the specified transformed data.
154     *
155     * @param dataRI the unscaled transformed data
156     * @param normalization the normalization to be applied
157     * @param type the type of transform (forward, inverse) which resulted in the specified data
158     */
159    private static void normalizeTransformedData(final double[][] dataRI,
160        final DftNormalization normalization, final TransformType type) {
161
162        final double[] dataR = dataRI[0];
163        final double[] dataI = dataRI[1];
164        final int n = dataR.length;
165        assert dataI.length == n;
166
167        switch (normalization) {
168            case STANDARD:
169                if (type == TransformType.INVERSE) {
170                    final double scaleFactor = 1.0 / ((double) n);
171                    for (int i = 0; i < n; i++) {
172                        dataR[i] *= scaleFactor;
173                        dataI[i] *= scaleFactor;
174                    }
175                }
176                break;
177            case UNITARY:
178                final double scaleFactor = 1.0 / FastMath.sqrt(n);
179                for (int i = 0; i < n; i++) {
180                    dataR[i] *= scaleFactor;
181                    dataI[i] *= scaleFactor;
182                }
183                break;
184            default:
185                /*
186                 * This should never occur in normal conditions. However this
187                 * clause has been added as a safeguard if other types of
188                 * normalizations are ever implemented, and the corresponding
189                 * test is forgotten in the present switch.
190                 */
191                throw new MathIllegalStateException();
192        }
193    }
194
195    /**
196     * Computes the standard transform of the specified complex data. The
197     * computation is done in place. The input data is laid out as follows
198     * <ul>
199     *   <li>{@code dataRI[0][i]} is the real part of the {@code i}-th data point,</li>
200     *   <li>{@code dataRI[1][i]} is the imaginary part of the {@code i}-th data point.</li>
201     * </ul>
202     *
203     * @param dataRI the two dimensional array of real and imaginary parts of the data
204     * @param normalization the normalization to be applied to the transformed data
205     * @param type the type of transform (forward, inverse) to be performed
206     * @throws DimensionMismatchException if the number of rows of the specified
207     *   array is not two, or the array is not rectangular
208     * @throws MathIllegalArgumentException if the number of data points is not
209     *   a power of two
210     */
211    public static void transformInPlace(final double[][] dataRI,
212        final DftNormalization normalization, final TransformType type) {
213
214        if (dataRI.length != 2) {
215            throw new DimensionMismatchException(dataRI.length, 2);
216        }
217        final double[] dataR = dataRI[0];
218        final double[] dataI = dataRI[1];
219        if (dataR.length != dataI.length) {
220            throw new DimensionMismatchException(dataI.length, dataR.length);
221        }
222
223        final int n = dataR.length;
224        if (!ArithmeticUtils.isPowerOfTwo(n)) {
225            throw new MathIllegalArgumentException(
226                LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
227                Integer.valueOf(n));
228        }
229
230        if (n == 1) {
231            return;
232        } else if (n == 2) {
233            final double srcR0 = dataR[0];
234            final double srcI0 = dataI[0];
235            final double srcR1 = dataR[1];
236            final double srcI1 = dataI[1];
237
238            // X_0 = x_0 + x_1
239            dataR[0] = srcR0 + srcR1;
240            dataI[0] = srcI0 + srcI1;
241            // X_1 = x_0 - x_1
242            dataR[1] = srcR0 - srcR1;
243            dataI[1] = srcI0 - srcI1;
244
245            normalizeTransformedData(dataRI, normalization, type);
246            return;
247        }
248
249        bitReversalShuffle2(dataR, dataI);
250
251        // Do 4-term DFT.
252        if (type == TransformType.INVERSE) {
253            for (int i0 = 0; i0 < n; i0 += 4) {
254                final int i1 = i0 + 1;
255                final int i2 = i0 + 2;
256                final int i3 = i0 + 3;
257
258                final double srcR0 = dataR[i0];
259                final double srcI0 = dataI[i0];
260                final double srcR1 = dataR[i2];
261                final double srcI1 = dataI[i2];
262                final double srcR2 = dataR[i1];
263                final double srcI2 = dataI[i1];
264                final double srcR3 = dataR[i3];
265                final double srcI3 = dataI[i3];
266
267                // 4-term DFT
268                // X_0 = x_0 + x_1 + x_2 + x_3
269                dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
270                dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
271                // X_1 = x_0 - x_2 + j * (x_3 - x_1)
272                dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1);
273                dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3);
274                // X_2 = x_0 - x_1 + x_2 - x_3
275                dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
276                dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
277                // X_3 = x_0 - x_2 + j * (x_1 - x_3)
278                dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3);
279                dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1);
280            }
281        } else {
282            for (int i0 = 0; i0 < n; i0 += 4) {
283                final int i1 = i0 + 1;
284                final int i2 = i0 + 2;
285                final int i3 = i0 + 3;
286
287                final double srcR0 = dataR[i0];
288                final double srcI0 = dataI[i0];
289                final double srcR1 = dataR[i2];
290                final double srcI1 = dataI[i2];
291                final double srcR2 = dataR[i1];
292                final double srcI2 = dataI[i1];
293                final double srcR3 = dataR[i3];
294                final double srcI3 = dataI[i3];
295
296                // 4-term DFT
297                // X_0 = x_0 + x_1 + x_2 + x_3
298                dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3;
299                dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3;
300                // X_1 = x_0 - x_2 + j * (x_3 - x_1)
301                dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3);
302                dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1);
303                // X_2 = x_0 - x_1 + x_2 - x_3
304                dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3;
305                dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3;
306                // X_3 = x_0 - x_2 + j * (x_1 - x_3)
307                dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1);
308                dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3);
309            }
310        }
311
312        int lastN0 = 4;
313        int lastLogN0 = 2;
314        while (lastN0 < n) {
315            int n0 = lastN0 << 1;
316            int logN0 = lastLogN0 + 1;
317            double wSubN0R = W_SUB_N_R[logN0];
318            double wSubN0I = W_SUB_N_I[logN0];
319            if (type == TransformType.INVERSE) {
320                wSubN0I = -wSubN0I;
321            }
322
323            // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2).
324            for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) {
325                int destOddStartIndex = destEvenStartIndex + lastN0;
326
327                double wSubN0ToRR = 1;
328                double wSubN0ToRI = 0;
329
330                for (int r = 0; r < lastN0; r++) {
331                    double grR = dataR[destEvenStartIndex + r];
332                    double grI = dataI[destEvenStartIndex + r];
333                    double hrR = dataR[destOddStartIndex + r];
334                    double hrI = dataI[destOddStartIndex + r];
335
336                    // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr
337                    dataR[destEvenStartIndex + r] = grR + wSubN0ToRR * hrR - wSubN0ToRI * hrI;
338                    dataI[destEvenStartIndex + r] = grI + wSubN0ToRR * hrI + wSubN0ToRI * hrR;
339                    // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr
340                    dataR[destOddStartIndex + r] = grR - (wSubN0ToRR * hrR - wSubN0ToRI * hrI);
341                    dataI[destOddStartIndex + r] = grI - (wSubN0ToRR * hrI + wSubN0ToRI * hrR);
342
343                    // WsubN0ToR *= WsubN0R
344                    double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I;
345                    double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R;
346                    wSubN0ToRR = nextWsubN0ToRR;
347                    wSubN0ToRI = nextWsubN0ToRI;
348                }
349            }
350
351            lastN0 = n0;
352            lastLogN0 = logN0;
353        }
354
355        normalizeTransformedData(dataRI, normalization, type);
356    }
357
358    /**
359     * Returns the (forward, inverse) transform of the specified real data set.
360     *
361     * @param f the real data array to be transformed
362     * @param type the type of transform (forward, inverse) to be performed
363     * @return the complex transformed array
364     * @throws MathIllegalArgumentException if the length of the data array is not a power of two
365     */
366    public Complex[] transform(final double[] f, final TransformType type) {
367        final double[][] dataRI = new double[][] {
368            MathArrays.copyOf(f, f.length), new double[f.length]
369        };
370
371        transformInPlace(dataRI, normalization, type);
372
373        return TransformUtils.createComplexArray(dataRI);
374    }
375
376    /**
377     * Returns the (forward, inverse) transform of the specified real function,
378     * sampled on the specified interval.
379     *
380     * @param f the function to be sampled and transformed
381     * @param min the (inclusive) lower bound for the interval
382     * @param max the (exclusive) upper bound for the interval
383     * @param n the number of sample points
384     * @param type the type of transform (forward, inverse) to be performed
385     * @return the complex transformed array
386     * @throws org.apache.commons.math4.exception.NumberIsTooLargeException
387     *   if the lower bound is greater than, or equal to the upper bound
388     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException
389     *   if the number of sample points {@code n} is negative
390     * @throws MathIllegalArgumentException if the number of sample points
391     *   {@code n} is not a power of two
392     */
393    public Complex[] transform(final UnivariateFunction f,
394                               final double min, final double max, final int n,
395                               final TransformType type) {
396
397        final double[] data = FunctionUtils.sample(f, min, max, n);
398        return transform(data, type);
399    }
400
401    /**
402     * Returns the (forward, inverse) transform of the specified complex data set.
403     *
404     * @param f the complex data array to be transformed
405     * @param type the type of transform (forward, inverse) to be performed
406     * @return the complex transformed array
407     * @throws MathIllegalArgumentException if the length of the data array is not a power of two
408     */
409    public Complex[] transform(final Complex[] f, final TransformType type) {
410        final double[][] dataRI = TransformUtils.createRealImaginaryArray(f);
411
412        transformInPlace(dataRI, normalization, type);
413
414        return TransformUtils.createComplexArray(dataRI);
415    }
416}