FFT: precompute twiddle factors

src/Benchmark/Transforms/FFT.cs

@@ -11,15 +11,17 @@ namespace Benchmark.Transforms
     {
         readonly Dictionary<int, Complex[]> _data = new Dictionary<int, Complex[]>();
 
-        [Params(32, 64, 128, 1024, 8192, 65536)]
+        // 64=1.5ms, 128=2.9ms, 2048=46.4ms, 8192=185.8ms, 65536=1.5s, 524288=11.9s (44.1 kHz)
+        //[Params(32, 64, 128, 1024, 2048, 4096, 8192, 65536, 1048576)]
+        [Params(128, 1024, 1048576)]
         public int N { get; set; }
 
         public FFT()
         {
-            var realSinusoidal = Generate.Sinusoidal(65536, 32, -2.0, 2.0);
+            var realSinusoidal = Generate.Sinusoidal(1048576, 32, -2.0, 2.0);
-            var imagSawtooth = Generate.Sawtooth(65536, 32, -20.0, 20.0);
+            var imagSawtooth = Generate.Sawtooth(1048576, 32, -20.0, 20.0);
             var signal = Generate.Map2(realSinusoidal, imagSawtooth, (r, i) => new Complex(r, i));
-            foreach (var n in new[] { 32, 64, 128, 1024, 8192, 65536 })
+            foreach (var n in new[] { 32, 64, 128, 1024, 2048, 4096, 8192, 65536, 1048576 })
             {
                 var s = new Complex[n];
                 Array.Copy(signal, 0, s, 0, n);

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.
         /// </summary>
         /// <param name="samples">Time-space sample vector.</param>
-        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <param name="positiveExponentSign">Fourier series exponent sign: true for positive, false for negative.</param>
-        internal static void Bluestein(Complex[] samples, int exponentSign)
+        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);
             }

src/Numerics/IntegralTransforms/Fourier.Naive.cs

@@ -45,11 +45,11 @@ namespace MathNet.Numerics.IntegralTransforms
         /// Naive generic DFT, useful e.g. to verify faster algorithms.
         /// </summary>
         /// <param name="samples">Time-space sample vector.</param>
-        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <param name="positiveExponentSign">Fourier series exponent sign: true for positive, false for negative.</param>
         /// <returns>Corresponding frequency-space vector.</returns>
-        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) =>

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.
         /// </summary>
         /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
-        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <param name="positiveExponentSign">Fourier series exponent sign: true for positive, false for negative.</param>
         /// <exception cref="ArgumentException"/>
-        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).
         /// </summary>
         /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
-        /// <param name="exponentSign">Fourier series exponent sign.</param>
+        /// <param name="positiveExponentSign">Fourier series exponent sign: true for positive, false for negative.</param>
         /// <exception cref="ArgumentException"/>
-        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;
+        }
     }
 }

src/Numerics/IntegralTransforms/Fourier.cs

@@ -392,7 +392,7 @@ namespace MathNet.Numerics.IntegralTransforms
         /// <returns>Corresponding frequency-space vector.</returns>
         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
         /// <returns>Corresponding time-space vector.</returns>
         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
         /// <exception cref="ArgumentException"/>
         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
         /// <exception cref="ArgumentException"/>
         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
         /// <param name="options">Fourier Transform Convention Options.</param>
         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
         /// <param name="options">Fourier Transform Convention Options.</param>
         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
         /// </summary>
         /// <param name="options">Fourier Transform Convention Options.</param>
         /// <returns>Fourier series exponent sign.</returns>
-        static int SignByOptions(FourierOptions options)
+        static bool PositiveSignByOptions(FourierOptions options)
         {
-            return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent ? 1 : -1;
+            return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent;
         }
 
         /// <summary>

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.Radix2(samples, false));
-            Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Parallel(samples, -1));
+            Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Parallel(samples, false));
         }
     }
 }