| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math4.transform; |
| |
| import java.util.Random; |
| |
| import org.apache.commons.numbers.complex.Complex; |
| import org.apache.commons.math4.analysis.UnivariateFunction; |
| import org.apache.commons.math4.analysis.function.Sin; |
| import org.apache.commons.math4.analysis.function.Sinc; |
| import org.apache.commons.math4.exception.MathIllegalArgumentException; |
| import org.apache.commons.math4.exception.NotStrictlyPositiveException; |
| import org.apache.commons.math4.exception.NumberIsTooLargeException; |
| import org.apache.commons.math4.transform.DftNormalization; |
| import org.apache.commons.math4.transform.FastFourierTransformer; |
| import org.apache.commons.math4.transform.TransformType; |
| import org.apache.commons.math4.transform.TransformUtils; |
| import org.apache.commons.math4.util.FastMath; |
| import org.junit.Assert; |
| import org.junit.Test; |
| |
| /** |
| * Test case for fast Fourier transformer. |
| * <p> |
| * FFT algorithm is exact, the small tolerance number is used only |
| * to account for round-off errors. |
| * |
| */ |
| public final class FastFourierTransformerTest { |
| /** The common seed of all random number generators used in this test. */ |
| private final static long SEED = 20110111L; |
| |
| /* |
| * Precondition checks. |
| */ |
| |
| @Test |
| public void testTransformComplexSizeNotAPowerOfTwo() { |
| final int n = 127; |
| final Complex[] x = createComplexData(n); |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(norm[i]); |
| try { |
| fft.transform(x, type[j]); |
| Assert.fail(norm[i] + ", " + type[j] + |
| ": MathIllegalArgumentException was expected"); |
| } catch (MathIllegalArgumentException e) { |
| // Expected behaviour |
| } |
| } |
| } |
| } |
| |
| @Test |
| public void testTransformRealSizeNotAPowerOfTwo() { |
| final int n = 127; |
| final double[] x = createRealData(n); |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(norm[i]); |
| try { |
| fft.transform(x, type[j]); |
| Assert.fail(norm[i] + ", " + type[j] + |
| ": MathIllegalArgumentException was expected"); |
| } catch (MathIllegalArgumentException e) { |
| // Expected behaviour |
| } |
| } |
| } |
| } |
| |
| @Test |
| public void testTransformFunctionSizeNotAPowerOfTwo() { |
| final int n = 127; |
| final UnivariateFunction f = new Sin(); |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(norm[i]); |
| try { |
| fft.transform(f, 0.0, Math.PI, n, type[j]); |
| Assert.fail(norm[i] + ", " + type[j] + |
| ": MathIllegalArgumentException was expected"); |
| } catch (MathIllegalArgumentException e) { |
| // Expected behaviour |
| } |
| } |
| } |
| } |
| |
| @Test |
| public void testTransformFunctionNotStrictlyPositiveNumberOfSamples() { |
| final int n = -128; |
| final UnivariateFunction f = new Sin(); |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(norm[i]); |
| try { |
| fft.transform(f, 0.0, Math.PI, n, type[j]); |
| fft.transform(f, 0.0, Math.PI, n, type[j]); |
| Assert.fail(norm[i] + ", " + type[j] + |
| ": NotStrictlyPositiveException was expected"); |
| } catch (NotStrictlyPositiveException e) { |
| // Expected behaviour |
| } |
| } |
| } |
| } |
| |
| @Test |
| public void testTransformFunctionInvalidBounds() { |
| final int n = 128; |
| final UnivariateFunction f = new Sin(); |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(norm[i]); |
| try { |
| fft.transform(f, Math.PI, 0.0, n, type[j]); |
| Assert.fail(norm[i] + ", " + type[j] + |
| ": NumberIsTooLargeException was expected"); |
| } catch (NumberIsTooLargeException e) { |
| // Expected behaviour |
| } |
| } |
| } |
| } |
| |
| /* |
| * Utility methods for checking (successful) transforms. |
| */ |
| |
| private static Complex[] createComplexData(final int n) { |
| final Random random = new Random(SEED); |
| final Complex[] data = new Complex[n]; |
| for (int i = 0; i < n; i++) { |
| final double re = 2.0 * random.nextDouble() - 1.0; |
| final double im = 2.0 * random.nextDouble() - 1.0; |
| data[i] = Complex.ofCartesian(re, im); |
| } |
| return data; |
| } |
| |
| private static double[] createRealData(final int n) { |
| final Random random = new Random(SEED); |
| final double[] data = new double[n]; |
| for (int i = 0; i < n; i++) { |
| data[i] = 2.0 * random.nextDouble() - 1.0; |
| } |
| return data; |
| } |
| |
| /** Naive implementation of DFT, for reference. */ |
| private static Complex[] dft(final Complex[] x, final int sgn) { |
| final int n = x.length; |
| final double[] cos = new double[n]; |
| final double[] sin = new double[n]; |
| final Complex[] y = new Complex[n]; |
| for (int i = 0; i < n; i++) { |
| final double arg = 2.0 * FastMath.PI * i / n; |
| cos[i] = FastMath.cos(arg); |
| sin[i] = FastMath.sin(arg); |
| } |
| for (int i = 0; i < n; i++) { |
| double yr = 0.0; |
| double yi = 0.0; |
| for (int j = 0; j < n; j++) { |
| final int index = (i * j) % n; |
| final double c = cos[index]; |
| final double s = sin[index]; |
| final double xr = x[j].getReal(); |
| final double xi = x[j].getImaginary(); |
| yr += c * xr - sgn * s * xi; |
| yi += sgn * s * xr + c * xi; |
| } |
| y[i] = Complex.ofCartesian(yr, yi); |
| } |
| return y; |
| } |
| |
| private static void doTestTransformComplex(final int n, final double tol, |
| final DftNormalization normalization, |
| final TransformType type) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(normalization); |
| final Complex[] x = createComplexData(n); |
| final Complex[] expected; |
| final double s; |
| if (type==TransformType.FORWARD) { |
| expected = dft(x, -1); |
| if (normalization == DftNormalization.STANDARD){ |
| s = 1.0; |
| } else { |
| s = 1.0 / FastMath.sqrt(n); |
| } |
| } else { |
| expected = dft(x, 1); |
| if (normalization == DftNormalization.STANDARD) { |
| s = 1.0 / n; |
| } else { |
| s = 1.0 / FastMath.sqrt(n); |
| } |
| } |
| final Complex[] actual = fft.transform(x, type); |
| for (int i = 0; i < n; i++) { |
| final String msg; |
| msg = String.format("%s, %s, %d, %d", normalization, type, n, i); |
| final double re = s * expected[i].getReal(); |
| Assert.assertEquals(msg, re, actual[i].getReal(), |
| tol * FastMath.abs(re)); |
| final double im = s * expected[i].getImaginary(); |
| Assert.assertEquals(msg, im, actual[i].getImaginary(), tol * |
| FastMath.abs(re)); |
| } |
| } |
| |
| private static void doTestTransformReal(final int n, final double tol, |
| final DftNormalization normalization, |
| final TransformType type) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(normalization); |
| final double[] x = createRealData(n); |
| final Complex[] xc = new Complex[n]; |
| for (int i = 0; i < n; i++) { |
| xc[i] = Complex.ofCartesian(x[i], 0.0); |
| } |
| final Complex[] expected; |
| final double s; |
| if (type == TransformType.FORWARD) { |
| expected = dft(xc, -1); |
| if (normalization == DftNormalization.STANDARD) { |
| s = 1.0; |
| } else { |
| s = 1.0 / FastMath.sqrt(n); |
| } |
| } else { |
| expected = dft(xc, 1); |
| if (normalization == DftNormalization.STANDARD) { |
| s = 1.0 / n; |
| } else { |
| s = 1.0 / FastMath.sqrt(n); |
| } |
| } |
| final Complex[] actual = fft.transform(x, type); |
| for (int i = 0; i < n; i++) { |
| final String msg; |
| msg = String.format("%s, %s, %d, %d", normalization, type, n, i); |
| final double re = s * expected[i].getReal(); |
| Assert.assertEquals(msg, re, actual[i].getReal(), |
| tol * FastMath.abs(re)); |
| final double im = s * expected[i].getImaginary(); |
| Assert.assertEquals(msg, im, actual[i].getImaginary(), tol * |
| FastMath.abs(re)); |
| } |
| } |
| |
| private static void doTestTransformFunction(final UnivariateFunction f, |
| final double min, final double max, int n, final double tol, |
| final DftNormalization normalization, |
| final TransformType type) { |
| final FastFourierTransformer fft; |
| fft = new FastFourierTransformer(normalization); |
| final Complex[] x = new Complex[n]; |
| for (int i = 0; i < n; i++) { |
| final double t = min + i * (max - min) / n; |
| x[i] = Complex.ofCartesian(f.value(t), 0); |
| } |
| final Complex[] expected; |
| final double s; |
| if (type == TransformType.FORWARD) { |
| expected = dft(x, -1); |
| if (normalization == DftNormalization.STANDARD) { |
| s = 1.0; |
| } else { |
| s = 1.0 / FastMath.sqrt(n); |
| } |
| } else { |
| expected = dft(x, 1); |
| if (normalization == DftNormalization.STANDARD) { |
| s = 1.0 / n; |
| } else { |
| s = 1.0 / FastMath.sqrt(n); |
| } |
| } |
| final Complex[] actual = fft.transform(f, min, max, n, type); |
| for (int i = 0; i < n; i++) { |
| final String msg = String.format("%d, %d", n, i); |
| final double re = s * expected[i].getReal(); |
| Assert.assertEquals(msg, re, actual[i].getReal(), |
| tol * FastMath.abs(re)); |
| final double im = s * expected[i].getImaginary(); |
| Assert.assertEquals(msg, im, actual[i].getImaginary(), tol * |
| FastMath.abs(re)); |
| } |
| } |
| |
| /* |
| * Tests of standard transform (when data is valid). |
| */ |
| |
| @Test |
| public void testTransformComplex() { |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| doTestTransformComplex(2, 1.0E-15, norm[i], type[j]); |
| doTestTransformComplex(4, 1.0E-14, norm[i], type[j]); |
| doTestTransformComplex(8, 1.0E-14, norm[i], type[j]); |
| doTestTransformComplex(16, 1.0E-13, norm[i], type[j]); |
| doTestTransformComplex(32, 1.0E-13, norm[i], type[j]); |
| doTestTransformComplex(64, 1.0E-12, norm[i], type[j]); |
| doTestTransformComplex(128, 1.0E-12, norm[i], type[j]); |
| } |
| } |
| } |
| |
| @Test |
| public void testStandardTransformReal() { |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| doTestTransformReal(2, 1.0E-15, norm[i], type[j]); |
| doTestTransformReal(4, 1.0E-14, norm[i], type[j]); |
| doTestTransformReal(8, 1.0E-14, norm[i], type[j]); |
| doTestTransformReal(16, 1.0E-13, norm[i], type[j]); |
| doTestTransformReal(32, 1.0E-13, norm[i], type[j]); |
| doTestTransformReal(64, 1.0E-13, norm[i], type[j]); |
| doTestTransformReal(128, 1.0E-11, norm[i], type[j]); |
| } |
| } |
| } |
| |
| @Test |
| public void testStandardTransformFunction() { |
| final UnivariateFunction f = new Sinc(); |
| final double min = -FastMath.PI; |
| final double max = FastMath.PI; |
| final DftNormalization[] norm; |
| norm = DftNormalization.values(); |
| final TransformType[] type; |
| type = TransformType.values(); |
| for (int i = 0; i < norm.length; i++) { |
| for (int j = 0; j < type.length; j++) { |
| doTestTransformFunction(f, min, max, 2, 1.0E-15, norm[i], type[j]); |
| doTestTransformFunction(f, min, max, 4, 1.0E-14, norm[i], type[j]); |
| doTestTransformFunction(f, min, max, 8, 1.0E-14, norm[i], type[j]); |
| doTestTransformFunction(f, min, max, 16, 1.0E-13, norm[i], type[j]); |
| doTestTransformFunction(f, min, max, 32, 1.0E-13, norm[i], type[j]); |
| doTestTransformFunction(f, min, max, 64, 1.0E-12, norm[i], type[j]); |
| doTestTransformFunction(f, min, max, 128, 1.0E-11, norm[i], type[j]); |
| } |
| } |
| } |
| |
| /* |
| * Additional tests for 1D data. |
| */ |
| |
| /** |
| * Test of transformer for the ad hoc data taken from Mathematica. |
| */ |
| @Test |
| public void testAdHocData() { |
| FastFourierTransformer transformer; |
| transformer = new FastFourierTransformer(DftNormalization.STANDARD); |
| Complex result[]; double tolerance = 1E-12; |
| |
| double x[] = {1.3, 2.4, 1.7, 4.1, 2.9, 1.7, 5.1, 2.7}; |
| Complex y[] = { |
| Complex.ofCartesian(21.9, 0.0), |
| Complex.ofCartesian(-2.09497474683058, 1.91507575950825), |
| Complex.ofCartesian(-2.6, 2.7), |
| Complex.ofCartesian(-1.10502525316942, -4.88492424049175), |
| Complex.ofCartesian(0.1, 0.0), |
| Complex.ofCartesian(-1.10502525316942, 4.88492424049175), |
| Complex.ofCartesian(-2.6, -2.7), |
| Complex.ofCartesian(-2.09497474683058, -1.91507575950825)}; |
| |
| result = transformer.transform(x, TransformType.FORWARD); |
| for (int i = 0; i < result.length; i++) { |
| Assert.assertEquals(y[i].getReal(), result[i].getReal(), tolerance); |
| Assert.assertEquals(y[i].getImaginary(), result[i].getImaginary(), tolerance); |
| } |
| |
| result = transformer.transform(y, TransformType.INVERSE); |
| for (int i = 0; i < result.length; i++) { |
| Assert.assertEquals(x[i], result[i].getReal(), tolerance); |
| Assert.assertEquals(0.0, result[i].getImaginary(), tolerance); |
| } |
| |
| double x2[] = {10.4, 21.6, 40.8, 13.6, 23.2, 32.8, 13.6, 19.2}; |
| TransformUtils.scaleArray(x2, 1.0 / FastMath.sqrt(x2.length)); |
| Complex y2[] = y; |
| |
| transformer = new FastFourierTransformer(DftNormalization.UNITARY); |
| result = transformer.transform(y2, TransformType.FORWARD); |
| for (int i = 0; i < result.length; i++) { |
| Assert.assertEquals(x2[i], result[i].getReal(), tolerance); |
| Assert.assertEquals(0.0, result[i].getImaginary(), tolerance); |
| } |
| |
| result = transformer.transform(x2, TransformType.INVERSE); |
| for (int i = 0; i < result.length; i++) { |
| Assert.assertEquals(y2[i].getReal(), result[i].getReal(), tolerance); |
| Assert.assertEquals(y2[i].getImaginary(), result[i].getImaginary(), tolerance); |
| } |
| } |
| |
| /** |
| * Test of transformer for the sine function. |
| */ |
| @Test |
| public void testSinFunction() { |
| UnivariateFunction f = new Sin(); |
| FastFourierTransformer transformer; |
| transformer = new FastFourierTransformer(DftNormalization.STANDARD); |
| Complex result[]; int N = 1 << 8; |
| double min, max, tolerance = 1E-12; |
| |
| min = 0.0; max = 2.0 * FastMath.PI; |
| result = transformer.transform(f, min, max, N, TransformType.FORWARD); |
| Assert.assertEquals(0.0, result[1].getReal(), tolerance); |
| Assert.assertEquals(-(N >> 1), result[1].getImaginary(), tolerance); |
| Assert.assertEquals(0.0, result[N-1].getReal(), tolerance); |
| Assert.assertEquals(N >> 1, result[N-1].getImaginary(), tolerance); |
| for (int i = 0; i < N-1; i += (i == 0 ? 2 : 1)) { |
| Assert.assertEquals(0.0, result[i].getReal(), tolerance); |
| Assert.assertEquals(0.0, result[i].getImaginary(), tolerance); |
| } |
| |
| min = -FastMath.PI; max = FastMath.PI; |
| result = transformer.transform(f, min, max, N, TransformType.INVERSE); |
| Assert.assertEquals(0.0, result[1].getReal(), tolerance); |
| Assert.assertEquals(-0.5, result[1].getImaginary(), tolerance); |
| Assert.assertEquals(0.0, result[N-1].getReal(), tolerance); |
| Assert.assertEquals(0.5, result[N-1].getImaginary(), tolerance); |
| for (int i = 0; i < N-1; i += (i == 0 ? 2 : 1)) { |
| Assert.assertEquals(0.0, result[i].getReal(), tolerance); |
| Assert.assertEquals(0.0, result[i].getImaginary(), tolerance); |
| } |
| } |
| } |