From f5f0354e253c6601d856a3eafa23ff665a8b00ff Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 4 Nov 2016 11:49:40 +0100 Subject: [PATCH] FFT: precompute twiddle factors --- .../IntegralTransforms/Fourier.Bluestein.cs | 16 ++--- .../IntegralTransforms/Fourier.Naive.cs | 6 +- .../IntegralTransforms/Fourier.RadixN.cs | 68 +++++++++++++++---- src/Numerics/IntegralTransforms/Fourier.cs | 16 ++--- .../IntegralTransformsTests/FourierTest.cs | 4 +- 5 files changed, 77 insertions(+), 33 deletions(-) diff --git a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs index 820c88e3..192240f7 100644 --- a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs +++ b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs @@ -107,7 +107,7 @@ namespace MathNet.Numerics.IntegralTransforms b[i] = sequence[m - i]; } - Radix2(b, -1); + Radix2(b, false); }, () => { @@ -117,7 +117,7 @@ namespace MathNet.Numerics.IntegralTransforms a[i] = sequence[i].Conjugate()*samples[i]; } - Radix2(a, -1); + Radix2(a, false); }); for (int i = 0; i < a.Length; i++) @@ -125,7 +125,7 @@ namespace MathNet.Numerics.IntegralTransforms a[i] *= b[i]; } - Radix2Parallel(a, 1); + Radix2Parallel(a, true); var nbinv = 1.0/m; for (int i = 0; i < samples.Length; i++) @@ -150,24 +150,24 @@ namespace MathNet.Numerics.IntegralTransforms /// Bluestein generic FFT for arbitrary sized sample vectors. /// /// Time-space sample vector. - /// Fourier series exponent sign. - internal static void Bluestein(Complex[] samples, int exponentSign) + /// Fourier series exponent sign: true for positive, false for negative. + internal static void Bluestein(Complex[] samples, bool positiveExponentSign) { int n = samples.Length; if (n.IsPowerOfTwo()) { - Radix2Parallel(samples, exponentSign); + Radix2Parallel(samples, positiveExponentSign); return; } - if (exponentSign == 1) + if (positiveExponentSign) { SwapRealImaginary(samples); } BluesteinConvolutionParallel(samples); - if (exponentSign == 1) + if (positiveExponentSign) { SwapRealImaginary(samples); } diff --git a/src/Numerics/IntegralTransforms/Fourier.Naive.cs b/src/Numerics/IntegralTransforms/Fourier.Naive.cs index 609abbde..0e768201 100644 --- a/src/Numerics/IntegralTransforms/Fourier.Naive.cs +++ b/src/Numerics/IntegralTransforms/Fourier.Naive.cs @@ -45,11 +45,11 @@ namespace MathNet.Numerics.IntegralTransforms /// Naive generic DFT, useful e.g. to verify faster algorithms. /// /// Time-space sample vector. - /// Fourier series exponent sign. + /// Fourier series exponent sign: true for positive, false for negative. /// Corresponding frequency-space vector. - internal static Complex[] Naive(Complex[] samples, int exponentSign) + internal static Complex[] Naive(Complex[] samples, bool positiveExponentSign) { - var w0 = exponentSign*Constants.Pi2/samples.Length; + var w0 = positiveExponentSign ? Constants.Pi2/samples.Length : -Constants.Pi2/samples.Length; var spectrum = new Complex[samples.Length]; CommonParallel.For(0, samples.Length, (u, v) => diff --git a/src/Numerics/IntegralTransforms/Fourier.RadixN.cs b/src/Numerics/IntegralTransforms/Fourier.RadixN.cs index 98564ee5..e075a589 100644 --- a/src/Numerics/IntegralTransforms/Fourier.RadixN.cs +++ b/src/Numerics/IntegralTransforms/Fourier.RadixN.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2014 Math.NET +// Copyright (c) 2009-2016 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -97,9 +97,9 @@ namespace MathNet.Numerics.IntegralTransforms /// Radix-2 generic FFT for power-of-two sized sample vectors. /// /// Sample vector, where the FFT is evaluated in place. - /// Fourier series exponent sign. + /// Fourier series exponent sign: true for positive, false for negative. /// - internal static void Radix2(Complex[] samples, int exponentSign) + internal static void Radix2(Complex[] samples, bool positiveExponentSign) { if (!samples.Length.IsPowerOfTwo()) { @@ -107,11 +107,21 @@ namespace MathNet.Numerics.IntegralTransforms } Radix2Reorder(samples); - for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + + for (int level = 0, levelSize = 1; levelSize < samples.Length; level++, levelSize *= 2) { - for (var k = 0; k < levelSize; k++) + Complex[] twiddles = TwiddleFactors(level); + for (int k = 0; k < levelSize; k++) { - Radix2Step(samples, exponentSign, levelSize, k); + Complex twiddle = positiveExponentSign ? twiddles[k] : twiddles[k].Conjugate(); + int step = levelSize << 1; + for (int i = k; i < samples.Length; i += step) + { + Complex ai = samples[i]; + Complex t = twiddle * samples[i + levelSize]; + samples[i] = ai + t; + samples[i + levelSize] = ai - t; + } } } } @@ -120,9 +130,9 @@ namespace MathNet.Numerics.IntegralTransforms /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// /// Sample vector, where the FFT is evaluated in place. - /// Fourier series exponent sign. + /// Fourier series exponent sign: true for positive, false for negative. /// - internal static void Radix2Parallel(Complex[] samples, int exponentSign) + internal static void Radix2Parallel(Complex[] samples, bool positiveExponentSign) { if (!samples.Length.IsPowerOfTwo()) { @@ -130,18 +140,52 @@ namespace MathNet.Numerics.IntegralTransforms } Radix2Reorder(samples); - for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + + for (int level = 0, levelSize = 1; levelSize < samples.Length; level++, levelSize *= 2) { - var size = levelSize; + Complex[] twiddles = TwiddleFactors(level); + var step = levelSize << 1; + + int size = levelSize; CommonParallel.For(0, size, 64, (u, v) => { - for (int i = u; i < v; i++) + for (int k = u; k < v; k++) { - Radix2Step(samples, exponentSign, size, i); + Complex twiddle = twiddles[k]; + double tr = twiddle.Real; + double ti = positiveExponentSign ? twiddle.Imaginary : -twiddle.Imaginary; + for (var i = k; i < samples.Length; i += step) + { + var ai = samples[i]; + var b = samples[i + size]; + double ttr = tr * b.Real - ti * b.Imaginary; + double tti = ti * b.Real + tr * b.Imaginary; + samples[i] = new Complex(ai.Real + ttr, ai.Imaginary + tti); + samples[i + size] = new Complex(ai.Real - ttr, ai.Imaginary - tti); + } } }); } } + + static readonly Complex[][] Twiddle = new Complex[64][]; + + static Complex[] TwiddleFactors(int level) + { + Complex[] tf = Twiddle[level]; + if (tf == null) + { + tf = new Complex[level.PowerOfTwo()]; + for (int i = 0; i < tf.Length; i++) + { + double exponent = i * Constants.Pi / tf.Length; + tf[i] = new Complex(Math.Cos(exponent), Math.Sin(exponent)); + } + Twiddle[level] = tf; + } + + return tf; + } } } diff --git a/src/Numerics/IntegralTransforms/Fourier.cs b/src/Numerics/IntegralTransforms/Fourier.cs index 2364f9ba..ada232dc 100644 --- a/src/Numerics/IntegralTransforms/Fourier.cs +++ b/src/Numerics/IntegralTransforms/Fourier.cs @@ -392,7 +392,7 @@ namespace MathNet.Numerics.IntegralTransforms /// Corresponding frequency-space vector. public static Complex[] NaiveForward(Complex[] samples, FourierOptions options = FourierOptions.Default) { - var frequencySpace = Naive(samples, SignByOptions(options)); + var frequencySpace = Naive(samples, PositiveSignByOptions(options)); ForwardScaleByOptions(options, frequencySpace); return frequencySpace; } @@ -405,7 +405,7 @@ namespace MathNet.Numerics.IntegralTransforms /// Corresponding time-space vector. public static Complex[] NaiveInverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default) { - var timeSpace = Naive(spectrum, -SignByOptions(options)); + var timeSpace = Naive(spectrum, !PositiveSignByOptions(options)); InverseScaleByOptions(options, timeSpace); return timeSpace; } @@ -418,7 +418,7 @@ namespace MathNet.Numerics.IntegralTransforms /// public static void Radix2Forward(Complex[] samples, FourierOptions options = FourierOptions.Default) { - Radix2Parallel(samples, SignByOptions(options)); + Radix2Parallel(samples, PositiveSignByOptions(options)); ForwardScaleByOptions(options, samples); } @@ -430,7 +430,7 @@ namespace MathNet.Numerics.IntegralTransforms /// public static void Radix2Inverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default) { - Radix2Parallel(spectrum, -SignByOptions(options)); + Radix2Parallel(spectrum, !PositiveSignByOptions(options)); InverseScaleByOptions(options, spectrum); } @@ -441,7 +441,7 @@ namespace MathNet.Numerics.IntegralTransforms /// Fourier Transform Convention Options. public static void BluesteinForward(Complex[] samples, FourierOptions options = FourierOptions.Default) { - Bluestein(samples, SignByOptions(options)); + Bluestein(samples, PositiveSignByOptions(options)); ForwardScaleByOptions(options, samples); } @@ -452,7 +452,7 @@ namespace MathNet.Numerics.IntegralTransforms /// Fourier Transform Convention Options. public static void BluesteinInverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default) { - Bluestein(spectrum, -SignByOptions(options)); + Bluestein(spectrum, !PositiveSignByOptions(options)); InverseScaleByOptions(options, spectrum); } @@ -462,9 +462,9 @@ namespace MathNet.Numerics.IntegralTransforms /// /// Fourier Transform Convention Options. /// Fourier series exponent sign. - static int SignByOptions(FourierOptions options) + static bool PositiveSignByOptions(FourierOptions options) { - return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent ? 1 : -1; + return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent; } /// diff --git a/src/UnitTests/IntegralTransformsTests/FourierTest.cs b/src/UnitTests/IntegralTransformsTests/FourierTest.cs index 17443794..3f210364 100644 --- a/src/UnitTests/IntegralTransformsTests/FourierTest.cs +++ b/src/UnitTests/IntegralTransformsTests/FourierTest.cs @@ -98,8 +98,8 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Forward(samples, FourierOptions.Default)); Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Inverse(samples, FourierOptions.Default)); - Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2(samples, -1)); - Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Parallel(samples, -1)); + Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2(samples, false)); + Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Parallel(samples, false)); } } }