diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs new file mode 100644 index 00000000..6b4efb25 --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs @@ -0,0 +1,298 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// https://numerics.mathdotnet.com +// +// Copyright (c) 2009-2018 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Runtime.CompilerServices; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Threading; + +using Complex = System.Numerics.Complex; + +namespace MathNet.Numerics.Providers.FourierTransform +{ + internal partial class ManagedFourierTransformProvider + { + /// + /// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1 + /// will cause k*k in the Bluestein sequence to overflow (GH-286). + /// + const int BluesteinSequenceLengthThreshold = 46341; + + /// + /// Generate the bluestein sequence for the provided problem size. + /// + /// Number of samples. + /// Bluestein sequence exp(I*Pi*k^2/N) + private static Complex32[] BluesteinSequence32(int n) + { + double s = Constants.Pi / n; + var sequence = new Complex32[n]; + + // TODO: benchmark whether the second variation is significantly + // faster than the former one. If not just use the former one always. + if (n > BluesteinSequenceLengthThreshold) + { + for (int k = 0; k < sequence.Length; k++) + { + double t = (s * k) * k; + sequence[k] = new Complex32((float)Math.Cos(t), (float)Math.Sin(t)); + } + } + else + { + for (int k = 0; k < sequence.Length; k++) + { + double t = s * (k * k); + sequence[k] = new Complex32((float)Math.Cos(t), (float)Math.Sin(t)); + } + } + + return sequence; + } + + /// + /// Generate the bluestein sequence for the provided problem size. + /// + /// Number of samples. + /// Bluestein sequence exp(I*Pi*k^2/N) + private static Complex[] BluesteinSequence(int n) + { + double s = Constants.Pi / n; + var sequence = new Complex[n]; + + // TODO: benchmark whether the second variation is significantly + // faster than the former one. If not just use the former one always. + if (n > BluesteinSequenceLengthThreshold) + { + for (int k = 0; k < sequence.Length; k++) + { + double t = (s * k) * k; + sequence[k] = new Complex(Math.Cos(t), Math.Sin(t)); + } + } + else + { + for (int k = 0; k < sequence.Length; k++) + { + double t = s * (k * k); + sequence[k] = new Complex(Math.Cos(t), Math.Sin(t)); + } + } + + return sequence; + } + + /// + /// Convolution with the bluestein sequence (Parallel Version). + /// + /// Sample Vector. + private static void BluesteinConvolutionParallel(Complex32[] samples) + { + int n = samples.Length; + Complex32[] sequence = BluesteinSequence32(n); + + // Padding to power of two >= 2N–1 so we can apply Radix-2 FFT. + int m = ((n << 1) - 1).CeilingToPowerOfTwo(); + var b = new Complex32[m]; + var a = new Complex32[m]; + + CommonParallel.Invoke( + () => + { + // Build and transform padded sequence b_k = exp(I*Pi*k^2/N) + for (int i = 0; i < n; i++) + { + b[i] = sequence[i]; + } + + for (int i = m - n + 1; i < b.Length; i++) + { + b[i] = sequence[m - i]; + } + + Radix2(b, -1); + }, + () => + { + // Build and transform padded sequence a_k = x_k * exp(-I*Pi*k^2/N) + for (int i = 0; i < samples.Length; i++) + { + a[i] = sequence[i].Conjugate() * samples[i]; + } + + Radix2(a, -1); + }); + + for (int i = 0; i < a.Length; i++) + { + a[i] *= b[i]; + } + + Radix2Parallel(a, 1); + + var nbinv = 1.0f / m; + for (int i = 0; i < samples.Length; i++) + { + samples[i] = nbinv * sequence[i].Conjugate() * a[i]; + } + } + + /// + /// Convolution with the bluestein sequence (Parallel Version). + /// + /// Sample Vector. + private static void BluesteinConvolutionParallel(Complex[] samples) + { + int n = samples.Length; + Complex[] sequence = BluesteinSequence(n); + + // Padding to power of two >= 2N–1 so we can apply Radix-2 FFT. + int m = ((n << 1) - 1).CeilingToPowerOfTwo(); + var b = new Complex[m]; + var a = new Complex[m]; + + CommonParallel.Invoke( + () => + { + // Build and transform padded sequence b_k = exp(I*Pi*k^2/N) + for (int i = 0; i < n; i++) + { + b[i] = sequence[i]; + } + + for (int i = m - n + 1; i < b.Length; i++) + { + b[i] = sequence[m - i]; + } + + Radix2(b, -1); + }, + () => + { + // Build and transform padded sequence a_k = x_k * exp(-I*Pi*k^2/N) + for (int i = 0; i < samples.Length; i++) + { + a[i] = sequence[i].Conjugate() * samples[i]; + } + + Radix2(a, -1); + }); + + for (int i = 0; i < a.Length; i++) + { + a[i] *= b[i]; + } + + Radix2Parallel(a, 1); + + var nbinv = 1.0 / m; + for (int i = 0; i < samples.Length; i++) + { + samples[i] = nbinv * sequence[i].Conjugate() * a[i]; + } + } + + /// + /// Swap the real and imaginary parts of each sample. + /// + /// Sample Vector. + private static void SwapRealImaginary(Complex32[] samples) + { + for (int i = 0; i < samples.Length; i++) + { + samples[i] = new Complex32(samples[i].Imaginary, samples[i].Real); + } + } + + /// + /// Swap the real and imaginary parts of each sample. + /// + /// Sample Vector. + private static void SwapRealImaginary(Complex[] samples) + { + for (int i = 0; i < samples.Length; i++) + { + samples[i] = new Complex(samples[i].Imaginary, samples[i].Real); + } + } + + /// + /// Bluestein generic FFT for arbitrary sized sample vectors. + /// + /// Time-space sample vector. + /// Fourier series exponent sign. + private static void Bluestein(Complex32[] samples, int exponentSign) + { + int n = samples.Length; + if (n.IsPowerOfTwo()) + { + Radix2Parallel(samples, exponentSign); + return; + } + + if (exponentSign == 1) + { + SwapRealImaginary(samples); + } + + BluesteinConvolutionParallel(samples); + + if (exponentSign == 1) + { + SwapRealImaginary(samples); + } + } + + /// + /// Bluestein generic FFT for arbitrary sized sample vectors. + /// + /// Time-space sample vector. + /// Fourier series exponent sign. + private static void Bluestein(Complex[] samples, int exponentSign) + { + int n = samples.Length; + if (n.IsPowerOfTwo()) + { + Radix2Parallel(samples, exponentSign); + return; + } + + if (exponentSign == 1) + { + SwapRealImaginary(samples); + } + + BluesteinConvolutionParallel(samples); + + if (exponentSign == 1) + { + SwapRealImaginary(samples); + } + } + } +} diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs new file mode 100644 index 00000000..f8bb0d0e --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs @@ -0,0 +1,222 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// https://numerics.mathdotnet.com +// +// Copyright (c) 2009-2018 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Runtime.CompilerServices; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Threading; + +using Complex = System.Numerics.Complex; + +namespace MathNet.Numerics.Providers.FourierTransform +{ + internal partial class ManagedFourierTransformProvider + { + /// + /// Radix-2 Reorder Helper Method + /// + /// Sample type + /// Sample vector + private static void Radix2Reorder(T[] samples) + { + var j = 0; + for (var i = 0; i < samples.Length - 1; i++) + { + if (i < j) + { + var temp = samples[i]; + samples[i] = samples[j]; + samples[j] = temp; + } + + var m = samples.Length; + + do + { + m >>= 1; + j ^= m; + } + while ((j & m) == 0); + } + } + + /// + /// Radix-2 Step Helper Method + /// + /// Sample vector. + /// Fourier series exponent sign. + /// Level Group Size. + /// Index inside of the level. +#if !NET40 + [MethodImpl(MethodImplOptions.AggressiveInlining)] +#endif + private static void Radix2Step(Complex32[] samples, int exponentSign, int levelSize, int k) + { + // Twiddle Factor + var exponent = (exponentSign * k) * Constants.Pi / levelSize; + var w = new Complex32((float)Math.Cos(exponent), (float)Math.Sin(exponent)); + + var step = levelSize << 1; + for (var i = k; i < samples.Length; i += step) + { + var ai = samples[i]; + var t = w * samples[i + levelSize]; + samples[i] = ai + t; + samples[i + levelSize] = ai - t; + } + } + + /// + /// Radix-2 Step Helper Method + /// + /// Sample vector. + /// Fourier series exponent sign. + /// Level Group Size. + /// Index inside of the level. +#if !NET40 + [MethodImpl(MethodImplOptions.AggressiveInlining)] +#endif + private static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k) + { + // Twiddle Factor + var exponent = (exponentSign * k) * Constants.Pi / levelSize; + var w = new Complex(Math.Cos(exponent), Math.Sin(exponent)); + + var step = levelSize << 1; + for (var i = k; i < samples.Length; i += step) + { + var ai = samples[i]; + var t = w * samples[i + levelSize]; + samples[i] = ai + t; + samples[i + levelSize] = ai - t; + } + } + + /// + /// Radix-2 generic FFT for power-of-two sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier series exponent sign. + /// + private static void Radix2(Complex32[] samples, int exponentSign) + { + if (!samples.Length.IsPowerOfTwo()) + { + throw new ArgumentException(Resources.ArgumentPowerOfTwo); + } + + Radix2Reorder(samples); + for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + { + for (var k = 0; k < levelSize; k++) + { + Radix2Step(samples, exponentSign, levelSize, k); + } + } + } + + /// + /// Radix-2 generic FFT for power-of-two sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier series exponent sign. + /// + private static void Radix2(Complex[] samples, int exponentSign) + { + if (!samples.Length.IsPowerOfTwo()) + { + throw new ArgumentException(Resources.ArgumentPowerOfTwo); + } + + Radix2Reorder(samples); + for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + { + for (var k = 0; k < levelSize; k++) + { + Radix2Step(samples, exponentSign, levelSize, k); + } + } + } + + /// + /// 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. + /// + private static void Radix2Parallel(Complex32[] samples, int exponentSign) + { + if (!samples.Length.IsPowerOfTwo()) + { + throw new ArgumentException(Resources.ArgumentPowerOfTwo); + } + + Radix2Reorder(samples); + for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + { + var size = levelSize; + + CommonParallel.For(0, size, 64, (u, v) => + { + for (int i = u; i < v; i++) + { + Radix2Step(samples, exponentSign, size, i); + } + }); + } + } + + /// + /// 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. + /// + private static void Radix2Parallel(Complex[] samples, int exponentSign) + { + if (!samples.Length.IsPowerOfTwo()) + { + throw new ArgumentException(Resources.ArgumentPowerOfTwo); + } + + Radix2Reorder(samples); + for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + { + var size = levelSize; + + CommonParallel.For(0, size, 64, (u, v) => + { + for (int i = u; i < v; i++) + { + Radix2Step(samples, exponentSign, size, i); + } + }); + } + } + } +} diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 3ac171b3..ba6e48f9 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // https://numerics.mathdotnet.com // -// Copyright (c) 2009-2016 Math.NET +// Copyright (c) 2009-2018 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -27,12 +27,13 @@ // using System; -using System.Numerics; using MathNet.Numerics.IntegralTransforms; +using Complex = System.Numerics.Complex; + namespace MathNet.Numerics.Providers.FourierTransform { - internal class ManagedFourierTransformProvider : IFourierTransformProvider + internal partial class ManagedFourierTransformProvider : IFourierTransformProvider { /// /// Try to find out whether the provider is available, at least in principle. @@ -65,67 +66,97 @@ namespace MathNet.Numerics.Providers.FourierTransform public void Forward(Complex32[] samples, FourierTransformScaling scaling) { + Bluestein(samples, -1); + switch (scaling) { case FourierTransformScaling.SymmetricScaling: - Fourier.BluesteinForward(samples, FourierOptions.Default); + { + ForwardScaleByOptions(FourierOptions.Default, samples); break; + } case FourierTransformScaling.ForwardScaling: - // Only backward scaling can be expressed with options, hence the double-inverse - Fourier.BluesteinInverse(samples, FourierOptions.AsymmetricScaling | FourierOptions.InverseExponent); + { + InverseScaleByOptions(FourierOptions.AsymmetricScaling | FourierOptions.InverseExponent, samples); break; + } default: - Fourier.BluesteinForward(samples, FourierOptions.NoScaling); + { + ForwardScaleByOptions(FourierOptions.NoScaling, samples); break; + } } } public void Forward(Complex[] samples, FourierTransformScaling scaling) { + Bluestein(samples, -1); + switch (scaling) { case FourierTransformScaling.SymmetricScaling: - Fourier.BluesteinForward(samples, FourierOptions.Default); + { + ForwardScaleByOptions(FourierOptions.Default, samples); break; + } case FourierTransformScaling.ForwardScaling: - // Only backward scaling can be expressed with options, hence the double-inverse - Fourier.BluesteinInverse(samples, FourierOptions.AsymmetricScaling | FourierOptions.InverseExponent); + { + InverseScaleByOptions(FourierOptions.AsymmetricScaling | FourierOptions.InverseExponent, samples); break; + } default: - Fourier.BluesteinForward(samples, FourierOptions.NoScaling); + { + ForwardScaleByOptions(FourierOptions.NoScaling, samples); break; + } } } public void Backward(Complex32[] spectrum, FourierTransformScaling scaling) { + Bluestein(spectrum, 1); + switch (scaling) { case FourierTransformScaling.SymmetricScaling: - Fourier.BluesteinInverse(spectrum, FourierOptions.Default); + { + InverseScaleByOptions(FourierOptions.Default, spectrum); break; + } case FourierTransformScaling.BackwardScaling: - Fourier.BluesteinInverse(spectrum, FourierOptions.AsymmetricScaling); + { + InverseScaleByOptions(FourierOptions.AsymmetricScaling, spectrum); break; + } default: - Fourier.BluesteinInverse(spectrum, FourierOptions.NoScaling); + { + InverseScaleByOptions(FourierOptions.NoScaling, spectrum); break; + } } } public void Backward(Complex[] spectrum, FourierTransformScaling scaling) { + Bluestein(spectrum, 1); + switch (scaling) { case FourierTransformScaling.SymmetricScaling: - Fourier.BluesteinInverse(spectrum, FourierOptions.Default); + { + InverseScaleByOptions(FourierOptions.Default, spectrum); break; + } case FourierTransformScaling.BackwardScaling: - Fourier.BluesteinInverse(spectrum, FourierOptions.AsymmetricScaling); + { + InverseScaleByOptions(FourierOptions.AsymmetricScaling, spectrum); break; + } default: - Fourier.BluesteinInverse(spectrum, FourierOptions.NoScaling); + { + InverseScaleByOptions(FourierOptions.NoScaling, spectrum); break; + } } } @@ -270,5 +301,93 @@ namespace MathNet.Numerics.Providers.FourierTransform { throw new NotSupportedException(); } + + /// + /// Rescale FFT-the resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + private static void ForwardScaleByOptions(FourierOptions options, Complex32[] samples) + { + if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling || + (options & FourierOptions.AsymmetricScaling) == FourierOptions.AsymmetricScaling) + { + return; + } + + var scalingFactor = (float)Math.Sqrt(1.0 / samples.Length); + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + + /// + /// Rescale FFT-the resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + private static void ForwardScaleByOptions(FourierOptions options, Complex[] samples) + { + if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling || + (options & FourierOptions.AsymmetricScaling) == FourierOptions.AsymmetricScaling) + { + return; + } + + var scalingFactor = Math.Sqrt(1.0 / samples.Length); + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + + /// + /// Rescale the iFFT-resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + private static void InverseScaleByOptions(FourierOptions options, Complex32[] samples) + { + if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling) + { + return; + } + + var scalingFactor = (float)1.0 / samples.Length; + if ((options & FourierOptions.AsymmetricScaling) != FourierOptions.AsymmetricScaling) + { + scalingFactor = (float)Math.Sqrt(scalingFactor); + } + + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + + /// + /// Rescale the iFFT-resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + private static void InverseScaleByOptions(FourierOptions options, Complex[] samples) + { + if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling) + { + return; + } + + var scalingFactor = 1.0 / samples.Length; + if ((options & FourierOptions.AsymmetricScaling) != FourierOptions.AsymmetricScaling) + { + scalingFactor = Math.Sqrt(scalingFactor); + } + + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } } }