diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs index 6b4efb25..44da0aca 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs @@ -28,7 +28,6 @@ using System; using System.Runtime.CompilerServices; -using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; @@ -38,9 +37,9 @@ 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). - /// + /// 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; /// @@ -135,7 +134,7 @@ namespace MathNet.Numerics.Providers.FourierTransform b[i] = sequence[m - i]; } - Radix2(b, -1); + Radix2Forward(b); }, () => { @@ -145,7 +144,7 @@ namespace MathNet.Numerics.Providers.FourierTransform a[i] = sequence[i].Conjugate() * samples[i]; } - Radix2(a, -1); + Radix2Forward(a); }); for (int i = 0; i < a.Length; i++) @@ -153,7 +152,7 @@ namespace MathNet.Numerics.Providers.FourierTransform a[i] *= b[i]; } - Radix2Parallel(a, 1); + Radix2InverseParallel(a); var nbinv = 1.0f / m; for (int i = 0; i < samples.Length; i++) @@ -190,7 +189,7 @@ namespace MathNet.Numerics.Providers.FourierTransform b[i] = sequence[m - i]; } - Radix2(b, -1); + Radix2Forward(b); }, () => { @@ -200,7 +199,7 @@ namespace MathNet.Numerics.Providers.FourierTransform a[i] = sequence[i].Conjugate() * samples[i]; } - Radix2(a, -1); + Radix2Forward(a); }); for (int i = 0; i < a.Length; i++) @@ -208,7 +207,7 @@ namespace MathNet.Numerics.Providers.FourierTransform a[i] *= b[i]; } - Radix2Parallel(a, 1); + Radix2InverseParallel(a); var nbinv = 1.0 / m; for (int i = 0; i < samples.Length; i++) @@ -244,55 +243,37 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Bluestein generic FFT for arbitrary sized sample vectors. /// - /// Time-space sample vector. - /// Fourier series exponent sign. - private static void Bluestein(Complex32[] samples, int exponentSign) + private static void BluesteinForward(Complex[] samples) { - 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) + private static void BluesteinInverse(Complex[] spectrum) { - int n = samples.Length; - if (n.IsPowerOfTwo()) - { - Radix2Parallel(samples, exponentSign); - return; - } - - if (exponentSign == 1) - { - SwapRealImaginary(samples); - } + SwapRealImaginary(spectrum); + BluesteinConvolutionParallel(spectrum); + SwapRealImaginary(spectrum); + } + /// + /// Bluestein generic FFT for arbitrary sized sample vectors. + /// + private static void BluesteinForward(Complex32[] samples) + { BluesteinConvolutionParallel(samples); + } - if (exponentSign == 1) - { - SwapRealImaginary(samples); - } + /// + /// Bluestein generic FFT for arbitrary sized sample vectors. + /// + private static void BluesteinInverse(Complex32[] spectrum) + { + SwapRealImaginary(spectrum); + BluesteinConvolutionParallel(spectrum); + SwapRealImaginary(spectrum); } } } diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs index f8bb0d0e..484c9a44 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs @@ -28,7 +28,6 @@ using System; using System.Runtime.CompilerServices; -using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; using Complex = System.Numerics.Complex; @@ -120,22 +119,29 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// 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) + private static void Radix2Forward(Complex32[] data) { - if (!samples.Length.IsPowerOfTwo()) + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { - throw new ArgumentException(Resources.ArgumentPowerOfTwo); + for (var k = 0; k < levelSize; k++) + { + Radix2Step(data, -1, levelSize, k); + } } + } - Radix2Reorder(samples); - for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + /// + /// Radix-2 generic FFT for power-of-two sized sample vectors. + /// + private static void Radix2Forward(Complex[] data) + { + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { for (var k = 0; k < levelSize; k++) { - Radix2Step(samples, exponentSign, levelSize, k); + Radix2Step(data, -1, levelSize, k); } } } @@ -143,22 +149,29 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// 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) + private static void Radix2Inverse(Complex32[] data) { - if (!samples.Length.IsPowerOfTwo()) + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { - throw new ArgumentException(Resources.ArgumentPowerOfTwo); + for (var k = 0; k < levelSize; k++) + { + Radix2Step(data, 1, levelSize, k); + } } + } - Radix2Reorder(samples); - for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + /// + /// Radix-2 generic FFT for power-of-two sized sample vectors. + /// + private static void Radix2Inverse(Complex[] data) + { + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { for (var k = 0; k < levelSize; k++) { - Radix2Step(samples, exponentSign, levelSize, k); + Radix2Step(data, 1, levelSize, k); } } } @@ -166,18 +179,30 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// 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) + private static void Radix2ForwardParallel(Complex32[] data) { - if (!samples.Length.IsPowerOfTwo()) + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { - throw new ArgumentException(Resources.ArgumentPowerOfTwo); + var size = levelSize; + + CommonParallel.For(0, size, 64, (u, v) => + { + for (int i = u; i < v; i++) + { + Radix2Step(data, -1, size, i); + } + }); } + } - Radix2Reorder(samples); - for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + /// + /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). + /// + private static void Radix2ForwardParallel(Complex[] data) + { + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { var size = levelSize; @@ -185,7 +210,7 @@ namespace MathNet.Numerics.Providers.FourierTransform { for (int i = u; i < v; i++) { - Radix2Step(samples, exponentSign, size, i); + Radix2Step(data, -1, size, i); } }); } @@ -194,18 +219,30 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// 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) + private static void Radix2InverseParallel(Complex32[] data) { - if (!samples.Length.IsPowerOfTwo()) + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { - throw new ArgumentException(Resources.ArgumentPowerOfTwo); + var size = levelSize; + + CommonParallel.For(0, size, 64, (u, v) => + { + for (int i = u; i < v; i++) + { + Radix2Step(data, 1, size, i); + } + }); } + } - Radix2Reorder(samples); - for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) + /// + /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). + /// + private static void Radix2InverseParallel(Complex[] data) + { + Radix2Reorder(data); + for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) { var size = levelSize; @@ -213,7 +250,7 @@ namespace MathNet.Numerics.Providers.FourierTransform { for (int i = u; i < v; i++) { - Radix2Step(samples, exponentSign, size, i); + Radix2Step(data, 1, size, i); } }); } diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 4b484bbd..68a4725c 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -65,7 +65,21 @@ namespace MathNet.Numerics.Providers.FourierTransform public void Forward(Complex32[] samples, FourierTransformScaling scaling) { - Bluestein(samples, -1); + if (samples.Length.IsPowerOfTwo()) + { + if (samples.Length >= 1024) + { + Radix2ForwardParallel(samples); + } + else + { + Radix2Forward(samples); + } + } + else + { + BluesteinForward(samples); + } switch (scaling) { @@ -84,7 +98,21 @@ namespace MathNet.Numerics.Providers.FourierTransform public void Forward(Complex[] samples, FourierTransformScaling scaling) { - Bluestein(samples, -1); + if (samples.Length.IsPowerOfTwo()) + { + if (samples.Length >= 1024) + { + Radix2ForwardParallel(samples); + } + else + { + Radix2Forward(samples); + } + } + else + { + BluesteinForward(samples); + } switch (scaling) { @@ -103,7 +131,21 @@ namespace MathNet.Numerics.Providers.FourierTransform public void Backward(Complex32[] spectrum, FourierTransformScaling scaling) { - Bluestein(spectrum, 1); + if (spectrum.Length.IsPowerOfTwo()) + { + if (spectrum.Length >= 1024) + { + Radix2InverseParallel(spectrum); + } + else + { + Radix2Inverse(spectrum); + } + } + else + { + BluesteinInverse(spectrum); + } switch (scaling) { @@ -122,7 +164,21 @@ namespace MathNet.Numerics.Providers.FourierTransform public void Backward(Complex[] spectrum, FourierTransformScaling scaling) { - Bluestein(spectrum, 1); + if (spectrum.Length.IsPowerOfTwo()) + { + if (spectrum.Length >= 1024) + { + Radix2InverseParallel(spectrum); + } + else + { + Radix2Inverse(spectrum); + } + } + else + { + BluesteinInverse(spectrum); + } switch (scaling) {