IterativeFFT.cs

Go to the documentation of this file.

using System;
using System.Diagnostics;

namespace FFT
{
   ///
   /// Iterative Discrete Fourier Transformation.
   ///
   public class IterativeFFT : IDFT
   {
      ///
      /// Direction of DFT.
      ///
      protected enum Direction : int {
         Forward, ///< Indicating the standard DFT.
         Backward ///< Indicating the inverse DFT.
      }

      ///
      /// \copydoc FFT::IDFT::DFT(Complex[])
      ///
      public Complex[] DFT(Complex[] a)
      {
         return FFT(a, Direction.Forward);
      }

      ///
      /// \copydoc FFT::IDFT::InverseDFT(Complex[])
      ///
      public Complex[] InverseDFT(Complex[] y)
      {
         return FFT(y, Direction.Backward);
      }

      ///
      /// Copy table and rearrange the elements. 
      ///
      /// The length <tt>n</tt> of the input table must be a power of 2.
      /// The index of an element <tt>a[i]</tt> in the output table 
      /// is the bit reverse of the binary representation <tt>i</tt> 
      /// (relative to the number of bits of <tt>n-1</tt>). For instance, if the 
      /// array <tt>a</tt> has 16 elements indexed with numbers 
      /// from <tt>0=0000B</tt> to <tt>15=1111B</tt>, then the 
      /// index of the element <tt>a[5]</tt> in the output is <tt>10</tt> 
      /// since <tt>5=0101B</tt> and <tt>10=1010B</tt>.
      ///
      /// \param a table whose length is a power of 2.
      ///
      /// Example: 
      /// \code
      ///  Complex[] a = {0,1,2,3,4,5,6,7};
      ///  Complex[] b = BitReverseCopy(a); // {0,4,2,6,1,5,3,7}
      ///  Complex[] c = {5,3,4,6,7,2,9,1};
      ///  Complex[] d = BitReverseCopy(c); // {5,7,4,9,3,2,6,1}
      /// \endcode
      ///
      protected internal Complex[] BitReverseCopy(Complex[] a)
      {
         Debug.Assert(a != null, "a is null");
         Debug.Assert(Utils.IsPowerOf2(a.Length),"length of a is not a power of 2!");

         Complex[] result = new Complex[a.Length];
         
         // Iteratively perform an increment of counter c in the reversed bit order
         // Start with the counter equal to 0. After each increment, copy an element of a
         // An example of sequence of counter increments for n=8:
         // i = 0 = 000B; c = 0 = 000B; result[0] = a[0];
         // i = 1 = 001B; c = 4 = 100B; result[4] = a[1];
         // i = 2 = 010B; c = 2 = 010B; result[2] = a[2];
         // i = 3 = 110B; c = 6 = 110B; result[6] = a[3];
         // i = 4 = 001B; c = 1 = 001B; result[1] = a[4];
         // i = 5 = 101B; c = 5 = 101B; result[5] = a[5];
         // i = 6 = 110B; c = 3 = 011B; result[3] = a[6];
         // i = 7 = 111B; c = 7 = 111B; result[7] = a[7];


         int n = a.Length;
         int c = 0;  // the counter increased in the reversed bit order
         int i = 0;  // the current index (the number of times counter has been increased)

         result[c] = a[i];
         while (c < n - 1) {

            int bit = n / 2;          //start with the "least" significant bit
            
            while ((bit & c) != 0) {  //while the current bit is set to 1, 
                                      //carry the increment to the following bit, i.e. 
               
               c -= bit;              // 1. set the bit to 0
               bit /= 2;              // 2. move to the "more" significant bit

            }                         //found the "least" significant bit set to 0

            c += bit;                 //set this bit to 1
            
            i++;
            
            result[c] = a[i];
         }
         
         return result;
      } // BitReverseCopy
      
      ///
      /// Computes DFT in the given direction.
      ///
      protected virtual Complex[] FFT(Complex[] a, Direction d)
      {
         int n = a.Length;
         int logN = Utils.Log2(n);
         int sgn = (d == Direction.Forward) ? (1) : (-1);
         
         a = BitReverseCopy(a);

         for (int s = 1; s < logN + 1; s++) {
            int m = Utils.Exp2(s);

            for (int j = 0; j < m / 2; j++) {
               Complex twiddle = Complex.FromPolar(1, sgn * 2 * j * Math.PI / m);
               for (int k = j; k < n; k += m) {
                  Complex t = twiddle * a[k + m / 2];
                  Complex u = a[k];
                  a[k] = u + t;
                  a[k + m / 2] = u - t;
               }
            }
         }

         if (d == Direction.Backward) {
            for (int i = 0; i < n; i++) {
               a[i] = a[i] / n;
            }
         }
         
         return a;
      } //FFT

      ///
      /// \copybrief FFT::IDFT::ToString()
      ///
      /// Returns the string <tt>"Iterative"</tt>
      public override string ToString()
      {
         return "Iterative";
      }
      
   } //class IterativeDFT
   
} //namespace FFT

---