﻿ Chris Lomont's C# Fast Fourier Transform code
``` // Code to implement decently performing FFT for complex and real valued                                          // signals. See www.lomont.org for a derivation of the relevant algorithms                                        // from first principles. Copyright Chris Lomont 2010-2012.                                                       // This code and any ports are free for all to use for any reason as long                                         // as this header is left in place.                                                                               // Version 1.1, Sept 2011                                                                                         using System;                                                                                                                                                                                                                       /* History:                                                                                                        * Sep 2011 - v1.1 - added parameters to support various sign conventions                                          *                   set via properties A and B.                                                                   *                 - Removed dependencies on LINQ and generic collections.                                         *                 - Added unit tests for the new properties.                                                      *                 - Switched UnitTest to static.                                                                  * Jan 2010 - v1.0 - Initial release                                                                               */                                                                                                               namespace Lomont                                                                                                  {                                                                                                                     /// <summary>                                                                                                     /// Represent a class that performs real or complex valued Fast Fourier                                           /// Transforms. Instantiate it and use the FFT or TableFFT methods to                                             /// compute complex to complex FFTs. Use FFTReal for real to complex                                              /// FFTs which are much faster than standard complex to complex FFTs.                                             /// Properties A and B allow selecting various FFT sign and scaling                                               /// conventions.                                                                                                  /// </summary>                                                                                                    public class LomontFFT                                                                                            {                                                                                                                     /// <summary>                                                                                                     /// Compute the forward or inverse Fourier Transform of data, with                                                /// data containing complex valued data as alternating real and                                                   /// imaginary parts. The length must be a power of 2. The data is                                                 /// modified in place.                                                                                            /// </summary>                                                                                                    /// <param name="data">The complex data stored as alternating real                                                /// and imaginary parts</param>                                                                                   /// <param name="forward">true for a forward transform, false for                                                 /// inverse transform</param>                                                                                     public void FFT(double[] data, bool forward)                                                                      {                                                                                                                     var n = data.Length;                                                                                              // checks n is a power of 2 in 2's complement format                                                              if ((n & (n - 1)) != 0)                                                                                               throw new ArgumentException(                                                                                          "data length " + n + " in FFT is not a power of 2");                                                      n /= 2;    // n is the number of samples                                                                                                                                                                                            Reverse(data, n); // bit index data reversal                                                                                                                                                                                        // do transform: so single point transforms, then doubles, etc.                                                   double sign = forward ? B : -B;                                                                                   var mmax = 1;                                                                                                     while (n > mmax)                                                                                                  {                                                                                                                     var istep = 2 * mmax;                                                                                             var theta = sign * Math.PI / mmax;                                                                                double wr = 1, wi = 0;                                                                                            var wpr = Math.Cos(theta);                                                                                        var wpi = Math.Sin(theta);                                                                                        for (var m = 0; m < istep; m += 2)                                                                                {                                                                                                                     for (var k = m; k < 2 * n; k += 2 * istep)                                                                        {                                                                                                                     var j = k + istep;                                                                                                var tempr = wr * data[j] - wi * data[j + 1];                                                                      var tempi = wi * data[j] + wr * data[j + 1];                                                                      data[j] = data[k] - tempr;                                                                                        data[j + 1] = data[k + 1] - tempi;                                                                                data[k] = data[k] + tempr;                                                                                        data[k + 1] = data[k + 1] + tempi;                                                                            }                                                                                                                 var t = wr; // trig recurrence                                                                                    wr = wr * wpr - wi * wpi;                                                                                         wi = wi * wpr + t * wpi;                                                                                      }                                                                                                                 mmax = istep;                                                                                                 }                                                                                                                                                                                                                                   // perform data scaling as needed                                                                                 Scale(data,n, forward);                                                                                       }                                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Compute the forward or inverse Fourier Transform of data, with data                                           /// containing complex valued data as alternating real and imaginary                                              /// parts. The length must be a power of 2. This method caches values                                             /// and should be slightly faster on than the FFT method for repeated uses.                                       /// It is also slightly more accurate. Data is transformed in place.                                              /// </summary>                                                                                                    /// <param name="data">The complex data stored as alternating real                                                /// and imaginary parts</param>                                                                                   /// <param name="forward">true for a forward transform, false for                                                 /// inverse transform</param>                                                                                     public void TableFFT(double[] data, bool forward)                                                                 {                                                                                                                     var n = data.Length;                                                                                              // checks n is a power of 2 in 2's complement format                                                              if ((n & (n - 1)) != 0)                                                                                               throw new ArgumentException(                                                                                          "data length " + n + " in FFT is not a power of 2"                                                                );                                                                                                        n /= 2;    // n is the number of samples                                                                                                                                                                                            Reverse(data, n); // bit index data reversal                                                                                                                                                                                        // make table if needed                                                                                           if ((cosTable == null) || (cosTable.Length != n))                                                                     Initialize(n);                                                                                                                                                                                                                  // do transform: so single point transforms, then doubles, etc.                                                   double sign = forward ? B : -B;                                                                                   var mmax = 1;                                                                                                     var tptr = 0;                                                                                                     while (n > mmax)                                                                                                  {                                                                                                                     var istep = 2 * mmax;                                                                                             for (var m = 0; m < istep; m += 2)                                                                                {                                                                                                                     var wr = cosTable[tptr];                                                                                          var wi = sign * sinTable[tptr++];                                                                                 for (var k = m; k < 2 * n; k += 2 * istep)                                                                        {                                                                                                                     var j = k + istep;                                                                                                var tempr = wr * data[j] - wi * data[j + 1];                                                                      var tempi = wi * data[j] + wr * data[j + 1];                                                                      data[j] = data[k] - tempr;                                                                                        data[j + 1] = data[k + 1] - tempi;                                                                                data[k] = data[k] + tempr;                                                                                        data[k + 1] = data[k + 1] + tempi;                                                                            }                                                                                                             }                                                                                                                 mmax = istep;                                                                                                 }                                                                                                                                                                                                                                                                                                                                                     // perform data scaling as needed                                                                                 Scale(data, n, forward);                                                                                      }                                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Compute the forward or inverse Fourier Transform of data, with                                                /// data containing real valued data only. The output is complex                                                  /// valued after the first two entries, stored in alternating real                                                /// and imaginary parts. The first two returned entries are the real                                              /// parts of the first and last value from the conjugate symmetric                                                /// output, which are necessarily real. The length must be a power                                                /// of 2.                                                                                                         /// </summary>                                                                                                    /// <param name="data">The complex data stored as alternating real                                                /// and imaginary parts</param>                                                                                   /// <param name="forward">true for a forward transform, false for                                                 /// inverse transform</param>                                                                                     public void RealFFT(double[] data, bool forward)                                                                  {                                                                                                                                                                                                                                       var n = data.Length; // # of real inputs, 1/2 the complex length                                                  // checks n is a power of 2 in 2's complement format                                                              if ((n & (n - 1)) != 0)                                                                                               throw new ArgumentException(                                                                                          "data length " + n + " in FFT is not a power of 2"                                                                );                                                                                                                                                                                                                          var sign = -1.0; // assume inverse FFT, this controls how algebra below works                                     if (forward)                                                                                                      { // do packed FFT. This can be changed to FFT to save memory                                                         TableFFT(data, true);                                                                                             sign = 1.0;                                                                                                       // scaling - divide by scaling for N/2, then mult by scaling for N                                                if (A != 1)                                                                                                       {                                                                                                                     var scale = Math.Pow(2.0, (A - 1) / 2.0);                                                                         for (var i = 0; i < data.Length; ++i)                                                                                 data[i] *= scale;                                                                                         }                                                                                                             }                                                                                                                                                                                                                                   var theta = B * sign * 2 * Math.PI / n;                                                                           var wpr = Math.Cos(theta);                                                                                        var wpi = Math.Sin(theta);                                                                                        var wjr = wpr;                                                                                                    var wji = wpi;                                                                                                                                                                                                                      for (var j = 1; j <= n/4; ++j)                                                                                    {                                                                                                                     var k = n / 2 - j;                                                                                                var tkr = data[2 * k];    // real and imaginary parts of t_k  = t_(n/2 - j)                                       var tki = data[2 * k + 1];                                                                                        var tjr = data[2 * j];    // real and imaginary parts of t_j                                                      var tji = data[2 * j + 1];                                                                                                                                                                                                          var a = (tjr - tkr) * wji;                                                                                        var b = (tji + tki) * wjr;                                                                                        var c = (tjr - tkr) * wjr;                                                                                        var d = (tji + tki) * wji;                                                                                        var e = (tjr + tkr);                                                                                              var f = (tji - tki);                                                                                                                                                                                                                // compute entry y[j]                                                                                             data[2 * j] = 0.5 * (e + sign * (a + b));                                                                         data[2 * j + 1] = 0.5 * (f + sign * (d - c));                                                                                                                                                                                       // compute entry y[k]                                                                                             data[2 * k] = 0.5 * (e - sign * (b + a));                                                                         data[2 * k + 1] = 0.5 * (sign * (d - c) - f);                                                                                                                                                                                       var temp = wjr;                                                                                                   // todo - allow more accurate version here? make option?                                                          wjr = wjr * wpr - wji * wpi;                                                                                      wji = temp * wpi + wji * wpr;                                                                                 }                                                                                                                                                                                                                                   if (forward)                                                                                                      {                                                                                                                     // compute final y0 and y_{N/2}, store in data[0], data[1]                                                        var temp = data[0];                                                                                               data[0] += data[1];                                                                                               data[1] = temp - data[1];                                                                                     }                                                                                                                 else                                                                                                              {                                                                                                                     var temp = data[0]; // unpack the y0 and y_{N/2}, then invert FFT                                                 data[0] = 0.5 * (temp + data[1]);                                                                                 data[1] = 0.5 * (temp - data[1]);                                                                                 // do packed inverse (table based) FFT. This can be changed to regular inverse FFT to save memory                 TableFFT(data, false);                                                                                            // scaling - divide by scaling for N, then mult by scaling for N/2                                                //if (A != -1) // todo - off by factor of 2? this works, but something seems weird                                {                                                                                                                     var scale = Math.Pow(2.0, -(A + 1) / 2.0)*2;                                                                      for (var i = 0; i < data.Length; ++i)                                                                                 data[i] *= scale;                                                                                         }                                                                                                             }                                                                                                             }                                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Determine how scaling works on the forward and inverse transforms.                                            /// For size N=2^n transforms, the forward transform gets divided by                                              /// N^((1-a)/2) and the inverse gets divided by N^((1+a)/2). Common                                               /// values for (A,B) are                                                                                          ///     ( 0, 1)  - default                                                                                        ///     (-1, 1)  - data processing                                                                                ///     ( 1,-1)  - signal processing                                                                              /// Usual values for A are 1, 0, or -1                                                                            /// </summary>                                                                                                    public int A { get; set; }                                                                                                                                                                                                          /// <summary>                                                                                                     /// Determine how phase works on the forward and inverse transforms.                                              /// For size N=2^n transforms, the forward transform uses an                                                      /// exp(B*2*pi/N) term and the inverse uses an exp(-B*2*pi/N) term.                                               /// Common values for (A,B) are                                                                                   ///     ( 0, 1)  - default                                                                                        ///     (-1, 1)  - data processing                                                                                ///     ( 1,-1)  - signal processing                                                                              /// Abs(B) should be relatively prime to N.                                                                       /// Setting B=-1 effectively corresponds to conjugating both input and                                            /// output data.                                                                                                  /// Usual values for B are 1 or -1.                                                                               /// </summary>                                                                                                    public int B { get; set; }                                                                                                                                                                                                          public LomontFFT()                                                                                                {                                                                                                                     A = 0;                                                                                                            B = 1;                                                                                                        }                                                                                                                                                                                                                                   #region Internals                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Scale data using n samples for forward and inverse transforms as needed                                       /// </summary>                                                                                                    /// <param name="data"></param>                                                                                   /// <param name="n"></param>                                                                                      /// <param name="forward"></param>                                                                                void Scale(double[] data, int n, bool forward)                                                                    {                                                                                                                     // forward scaling if needed                                                                                      if ((forward) && (A != 1))                                                                                        {                                                                                                                     var scale = Math.Pow(n, (A - 1) / 2.0);                                                                           for (var i = 0; i < data.Length; ++i)                                                                                 data[i] *= scale;                                                                                         }                                                                                                                                                                                                                                   // inverse scaling if needed                                                                                      if ((!forward) && (A != -1))                                                                                      {                                                                                                                     var scale = Math.Pow(n, -(A + 1) / 2.0);                                                                          for (var i = 0; i < data.Length; ++i)                                                                                 data[i] *= scale;                                                                                         }                                                                                                             }                                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Call this with the size before using the TableFFT version                                                     /// Fills in tables for speed. Done automatically in TableFFT                                                     /// </summary>                                                                                                    /// <param name="size">The size of the FFT in samples</param>                                                     void Initialize(int size)                                                                                         {                                                                                                                     // NOTE: if you port to non garbage collected languages                                                           // like C# or Java be sure to free these correctly                                                                cosTable = new double[size];                                                                                      sinTable = new double[size];                                                                                                                                                                                                        // forward pass                                                                                                   var n = size;                                                                                                     int mmax = 1, pos = 0;                                                                                            while (n > mmax)                                                                                                  {                                                                                                                     var istep = 2 * mmax;                                                                                             var theta = Math.PI / mmax;                                                                                       double wr = 1, wi = 0;                                                                                            var wpi = Math.Sin(theta);                                                                                        // compute in a slightly slower yet more accurate manner                                                          var wpr = Math.Sin(theta / 2);                                                                                    wpr = -2 * wpr * wpr;                                                                                             for (var m = 0; m < istep; m += 2)                                                                                {                                                                                                                     cosTable[pos] = wr;                                                                                               sinTable[pos++] = wi;                                                                                             var t = wr;                                                                                                       wr = wr * wpr - wi * wpi + wr;                                                                                    wi = wi * wpr + t * wpi + wi;                                                                                 }                                                                                                                 mmax = istep;                                                                                                 }                                                                                                             }                                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Swap data indices whenever index i has binary                                                                 /// digits reversed from index j, where data is                                                                   /// two doubles per index.                                                                                        /// </summary>                                                                                                    /// <param name="data"></param>                                                                                   /// <param name="n"></param>                                                                                      static void Reverse(double [] data, int n)                                                                        {                                                                                                                     // bit reverse the indices. This is exercise 5 in section                                                         // 7.2.1.1 of Knuth's TAOCP the idea is a binary counter                                                          // in k and one with bits reversed in j                                                                           int j = 0, k = 0; // Knuth R1: initialize                                                                         var top = n / 2;  // this is Knuth's 2^(n-1)                                                                      while (true)                                                                                                      {                                                                                                                     // Knuth R2: swap - swap j+1 and k+2^(n-1), 2 entries each                                                        var t = data[j + 2];                                                                                              data[j + 2] = data[k + n];                                                                                        data[k + n] = t;                                                                                                  t = data[j + 3];                                                                                                  data[j + 3] = data[k + n + 1];                                                                                    data[k + n + 1] = t;                                                                                              if (j > k)                                                                                                        { // swap two more                                                                                                    // j and k                                                                                                        t = data[j];                                                                                                      data[j] = data[k];                                                                                                data[k] = t;                                                                                                      t = data[j + 1];                                                                                                  data[j + 1] = data[k + 1];                                                                                        data[k + 1] = t;                                                                                                  // j + top + 1 and k+top + 1                                                                                      t = data[j + n + 2];                                                                                              data[j + n + 2] = data[k + n + 2];                                                                                data[k + n + 2] = t;                                                                                              t = data[j + n + 3];                                                                                              data[j + n + 3] = data[k + n + 3];                                                                                data[k + n + 3] = t;                                                                                          }                                                                                                                 // Knuth R3: advance k                                                                                            k += 4;                                                                                                           if (k >= n)                                                                                                           break;                                                                                                        // Knuth R4: advance j                                                                                            var h = top;                                                                                                      while (j >= h)                                                                                                    {                                                                                                                     j -= h;                                                                                                           h /= 2;                                                                                                       }                                                                                                                 j += h;                                                                                                       } // bit reverse loop                                                                                         }                                                                                                                                                                                                                                   /// <summary>                                                                                                     /// Pre-computed sine/cosine tables for speed                                                                     /// </summary>                                                                                                    double [] cosTable;                                                                                               double [] sinTable;                                                                                                                                                                                                                 #endregion                                                                                                                                                                                                                      }                                                                                                             }                                                                                                                 // end of file                                                                                                                                                                                                                                                                                                                                        ```