FFT: migrate radix2 and bluestein algorithms into managed provider

src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs

@@ -0,0 +1,298 @@
+// <copyright file="ManagedFourierTransformProvider.Bluestein.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+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
+    {
+        /// <summary>
+     /// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1
+     /// will cause k*k in the Bluestein sequence to overflow (GH-286).
+     /// </summary>
+        const int BluesteinSequenceLengthThreshold = 46341;
+
+        /// <summary>
+        /// Generate the bluestein sequence for the provided problem size.
+        /// </summary>
+        /// <param name="n">Number of samples.</param>
+        /// <returns>Bluestein sequence exp(I*Pi*k^2/N)</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Generate the bluestein sequence for the provided problem size.
+        /// </summary>
+        /// <param name="n">Number of samples.</param>
+        /// <returns>Bluestein sequence exp(I*Pi*k^2/N)</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Convolution with the bluestein sequence (Parallel Version).
+        /// </summary>
+        /// <param name="samples">Sample Vector.</param>
+        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];
+            }
+        }
+
+        /// <summary>
+        /// Convolution with the bluestein sequence (Parallel Version).
+        /// </summary>
+        /// <param name="samples">Sample Vector.</param>
+        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];
+            }
+        }
+
+        /// <summary>
+        /// Swap the real and imaginary parts of each sample.
+        /// </summary>
+        /// <param name="samples">Sample Vector.</param>
+        private static void SwapRealImaginary(Complex32[] samples)
+        {
+            for (int i = 0; i < samples.Length; i++)
+            {
+                samples[i] = new Complex32(samples[i].Imaginary, samples[i].Real);
+            }
+        }
+
+        /// <summary>
+        /// Swap the real and imaginary parts of each sample.
+        /// </summary>
+        /// <param name="samples">Sample Vector.</param>
+        private static void SwapRealImaginary(Complex[] samples)
+        {
+            for (int i = 0; i < samples.Length; i++)
+            {
+                samples[i] = new Complex(samples[i].Imaginary, samples[i].Real);
+            }
+        }
+
+        /// <summary>
+        /// Bluestein generic FFT for arbitrary sized sample vectors.
+        /// </summary>
+        /// <param name="samples">Time-space sample vector.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        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);
+            }
+        }
+
+        /// <summary>
+        /// Bluestein generic FFT for arbitrary sized sample vectors.
+        /// </summary>
+        /// <param name="samples">Time-space sample vector.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        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);
+            }
+        }
+    }
+}

src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs

@@ -0,0 +1,222 @@
+// <copyright file="ManagedFourierTransformProvider.Radix2.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+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
+    {
+        /// <summary>
+        /// Radix-2 Reorder Helper Method
+        /// </summary>
+        /// <typeparam name="T">Sample type</typeparam>
+        /// <param name="samples">Sample vector</param>
+        private static void Radix2Reorder<T>(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);
+            }
+        }
+
+        /// <summary>
+        /// Radix-2 Step Helper Method
+        /// </summary>
+        /// <param name="samples">Sample vector.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <param name="levelSize">Level Group Size.</param>
+        /// <param name="k">Index inside of the level.</param>
+#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;
+            }
+        }
+
+        /// <summary>
+        /// Radix-2 Step Helper Method
+        /// </summary>
+        /// <param name="samples">Sample vector.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <param name="levelSize">Level Group Size.</param>
+        /// <param name="k">Index inside of the level.</param>
+#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;
+            }
+        }
+
+        /// <summary>
+        /// Radix-2 generic FFT for power-of-two sized sample vectors.
+        /// </summary>
+        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <exception cref="ArgumentException"/>
+        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);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Radix-2 generic FFT for power-of-two sized sample vectors.
+        /// </summary>
+        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <exception cref="ArgumentException"/>
+        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);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
+        /// </summary>
+        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <exception cref="ArgumentException"/>
+        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);
+                    }
+                });
+            }
+        }
+
+        /// <summary>
+        /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
+        /// </summary>
+        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
+        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <exception cref="ArgumentException"/>
+        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);
+                    }
+                });
+            }
+        }
+    }
+}

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 @@
 // </copyright>
 
 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
     {
         /// <summary>
         /// 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();
         }
+
+        /// <summary>
+        /// Rescale FFT-the resulting vector according to the provided convention options.
+        /// </summary>
+        /// <param name="options">Fourier Transform Convention Options.</param>
+        /// <param name="samples">Sample Vector.</param>
+        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;
+            }
+        }
+
+        /// <summary>
+        /// Rescale FFT-the resulting vector according to the provided convention options.
+        /// </summary>
+        /// <param name="options">Fourier Transform Convention Options.</param>
+        /// <param name="samples">Sample Vector.</param>
+        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;
+            }
+        }
+
+        /// <summary>
+        /// Rescale the iFFT-resulting vector according to the provided convention options.
+        /// </summary>
+        /// <param name="options">Fourier Transform Convention Options.</param>
+        /// <param name="samples">Sample Vector.</param>
+        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;
+            }
+        }
+
+        /// <summary>
+        /// Rescale the iFFT-resulting vector according to the provided convention options.
+        /// </summary>
+        /// <param name="options">Fourier Transform Convention Options.</param>
+        /// <param name="samples">Sample Vector.</param>
+        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;
+            }
+        }
     }
 }