Merge pull request #80 from swharden/77

Experimental: Bluesteins algorithm

src/FftSharp.Tests/ExperimentalTests.cs

@@ -25,4 +25,24 @@ public void Test_Experimental_DFT()
             Assert.AreEqual(expectedImag[i], result[i].Imaginary, 1e-5);
         }
     }
+
+    [Test]
+    [System.Obsolete]
+    public void Test_Experimental_Bluestein()
+    {
+        double[] prime = { 1, 2, 3, 4, 5, 6, 7 };
+        System.Numerics.Complex[] result = FftSharp.Experimental.Bluestein(prime);
+
+        // tested with python
+        double[] expectedReal = { 28.00000, -3.50000, -3.50000, -3.50000, -3.50000, -3.50000, -3.50000 };
+        double[] expectedImag = { -0.00000, 7.26782, 2.79116, 0.79885, -0.79885, -2.79116, -7.26782 };
+
+        Assert.AreEqual(expectedReal.Length, result.Length);
+
+        for (int i = 0; i < result.Length; i++)
+        {
+            Assert.AreEqual(expectedReal[i], result[i].Real, 1e-5);
+            Assert.AreEqual(expectedImag[i], result[i].Imaginary, 1e-5);
+        }
+    }
 }

src/FftSharp/BluesteinOperations.cs

@@ -0,0 +1,145 @@
+using System;
+
+namespace FftSharp;
+
+internal static class BluesteinOperations
+{
+    /* 
+	 * Computes the discrete Fourier transform (DFT) or inverse transform of the given complex vector, storing the result back into the vector.
+	 * The vector can have any length. This is a wrapper function. The inverse transform does not perform scaling, so it is not a true inverse.
+	 */
+    public static void Transform(System.Numerics.Complex[] vec, bool inverse)
+    {
+        int n = vec.Length;
+        if (n == 0)
+            return;
+        else if ((n & (n - 1)) == 0)  // Is power of 2
+            TransformRadix2(vec, inverse);
+        else  // More complicated algorithm for arbitrary sizes
+            TransformBluestein(vec, inverse);
+    }
+
+
+    /* 
+	 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
+	 * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm.
+	 */
+    public static void TransformRadix2(System.Numerics.Complex[] vec, bool inverse)
+    {
+        // Length variables
+        int n = vec.Length;
+        int levels = 0;  // compute levels = floor(log2(n))
+        for (int temp = n; temp > 1; temp >>= 1)
+            levels++;
+        if (1 << levels != n)
+            throw new ArgumentException("Length is not a power of 2");
+
+        // Trigonometric table
+        System.Numerics.Complex[] expTable = new System.Numerics.Complex[n / 2];
+        double coef = 2 * Math.PI / n * (inverse ? 1 : -1);
+        for (int i = 0; i < n / 2; i++)
+            expTable[i] = System.Numerics.Complex.FromPolarCoordinates(1, i * coef);
+
+        // Bit-reversed addressing permutation
+        for (int i = 0; i < n; i++)
+        {
+            int j = ReverseBits(i, levels);
+            if (j > i)
+            {
+                System.Numerics.Complex temp = vec[i];
+                vec[i] = vec[j];
+                vec[j] = temp;
+            }
+        }
+
+        // Cooley-Tukey decimation-in-time radix-2 FFT
+        for (int size = 2; size <= n; size *= 2)
+        {
+            int halfsize = size / 2;
+            int tablestep = n / size;
+            for (int i = 0; i < n; i += size)
+            {
+                for (int j = i, k = 0; j < i + halfsize; j++, k += tablestep)
+                {
+                    System.Numerics.Complex temp = vec[j + halfsize] * expTable[k];
+                    vec[j + halfsize] = vec[j] - temp;
+                    vec[j] += temp;
+                }
+            }
+            if (size == n)  // Prevent overflow in 'size *= 2'
+                break;
+        }
+    }
+
+
+    /* 
+	 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
+	 * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function.
+	 * Uses Bluestein's chirp z-transform algorithm.
+	 */
+    public static void TransformBluestein(System.Numerics.Complex[] vec, bool inverse)
+    {
+        // Find a power-of-2 convolution length m such that m >= n * 2 + 1
+        int n = vec.Length;
+        if (n >= 0x20000000)
+            throw new ArgumentException("Array too large");
+        int m = 1;
+        while (m < n * 2 + 1)
+            m *= 2;
+
+        // Trigonometric table
+        System.Numerics.Complex[] expTable = new System.Numerics.Complex[n];
+        double coef = Math.PI / n * (inverse ? 1 : -1);
+        for (int i = 0; i < n; i++)
+        {
+            int j = (int)((long)i * i % (n * 2));  // This is more accurate than j = i * i
+            expTable[i] = System.Numerics.Complex.Exp(new System.Numerics.Complex(0, j * coef));
+        }
+
+        // Temporary vectors and preprocessing
+        System.Numerics.Complex[] avec = new System.Numerics.Complex[m];
+        for (int i = 0; i < n; i++)
+            avec[i] = vec[i] * expTable[i];
+        System.Numerics.Complex[] bvec = new System.Numerics.Complex[m];
+        bvec[0] = expTable[0];
+        for (int i = 1; i < n; i++)
+            bvec[i] = bvec[m - i] = System.Numerics.Complex.Conjugate(expTable[i]);
+
+        // Convolution
+        System.Numerics.Complex[] cvec = new System.Numerics.Complex[m];
+        Convolve(avec, bvec, cvec);
+
+        // Postprocessing
+        for (int i = 0; i < n; i++)
+            vec[i] = cvec[i] * expTable[i];
+    }
+
+
+    /* 
+	 * Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
+	 */
+    public static void Convolve(System.Numerics.Complex[] xvec, System.Numerics.Complex[] yvec, System.Numerics.Complex[] outvec)
+    {
+        int n = xvec.Length;
+        if (n != yvec.Length || n != outvec.Length)
+            throw new ArgumentException("Mismatched lengths");
+        xvec = (System.Numerics.Complex[])xvec.Clone();
+        yvec = (System.Numerics.Complex[])yvec.Clone();
+        Transform(xvec, false);
+        Transform(yvec, false);
+        for (int i = 0; i < n; i++)
+            xvec[i] *= yvec[i];
+        Transform(xvec, true);
+        for (int i = 0; i < n; i++)  // Scaling (because this FFT implementation omits it)
+            outvec[i] = xvec[i] / n;
+    }
+
+
+    private static int ReverseBits(int val, int width)
+    {
+        int result = 0;
+        for (int i = 0; i < width; i++, val >>= 1)
+            result = (result << 1) | (val & 1);
+        return result;
+    }
+}

src/FftSharp/Experimental.cs

@@ -70,5 +70,34 @@ public static System.Numerics.Complex[] DFT(double[] input, bool inverse = false
 
             return DFT(inputComplex, inverse);
         }
+
+        /// <summary>
+        /// Computes the discrete Fourier transform (DFT) of the given array.
+        /// The array may have any length.
+        /// Uses Bluestein's chirp z-transform algorithm.
+        /// </summary>
+        public static System.Numerics.Complex[] Bluestein(double[] real, bool inverse = false)
+        {
+            System.Numerics.Complex[] buffer = new System.Numerics.Complex[real.Length];
+
+            for (int i = 0; i < real.Length; i++)
+            {
+                buffer[i] = new System.Numerics.Complex(real[i], 0);
+            }
+
+            Bluestein(buffer, inverse);
+
+            return buffer;
+        }
+
+        /// <summary>
+        /// Computes the discrete Fourier transform (DFT) of the given complex vector in place.
+        /// The vector can have any length.
+        /// Uses Bluestein's chirp z-transform algorithm.
+        /// </summary>
+        public static void Bluestein(System.Numerics.Complex[] buffer, bool inverse = false)
+        {
+            BluesteinOperations.TransformBluestein(buffer, inverse);
+        }
     }
 }