RecursiveFFT.cs

Go to the documentation of this file.

using System;
using System.Diagnostics;

namespace FFT {

   /// 
   /// Recursive Fast Fourier Transformation. 
   ///
   /// Time complexity is \f$ O(n\log_2 n) \f$.
   ///
   public class RecursiveFFT : IDFT
   {
      /// 
      /// \copydoc FFT::IDFT::DFT(Complex[] a)
      ///
      public Complex[] DFT(Complex[] a)
      {
         Debug.Assert(a != null, "a is null");
         Debug.Assert(Utils.IsPowerOf2(a.Length),
                      "length of a is not a power of 2!");
         
         int n = a.Length;
         
         if (n == 1) {
            return a;
         }
         
         // First, we divide the input in two
         Complex[] a0 = new Complex[n/2];
         Complex[] a1 = new Complex[n/2];
         for (int k = 0 ; k < n/2; k++) {
            a0[k] = a[2*k];
            a1[k] = a[2*k + 1];
         }
         

         // Now, recursively handling the two pieces
         Complex[] y0 = DFT(a0);
         Complex[] y1 = DFT(a1);
         
         // Finally, constructing the result
         Complex[] y = new Complex[n];
         for (int k = 0 ; k < n/2; k++) {
            Complex twiddle = Complex.FromPolar(1,2*k*Math.PI/n);

            y[k]     = y0[k] + twiddle * y1[k];
            y[k+n/2] = y0[k] - twiddle * y1[k];
         }
         
         return y;

      } //DFT

      ///
      /// \copydoc FFT::IDFT::InverseDFT(Complex[])
      ///
      public Complex[] InverseDFT(Complex[] y)
      {
         Debug.Assert(y != null, "y is null");
         Debug.Assert(Utils.IsPowerOf2(y.Length),
                      "length of y is not a power of 2!");

         int n = y.Length;
         
         if (n == 1) {
            return y;
         }

         // First, spliting the input into a0 and a1
         Complex[] y0 = new Complex[n/2];
         Complex[] y1 = new Complex[n/2];
         for (int k = 0 ; k < n/2; k++) {
            y0[k] = y[2*k];
            y1[k] = y[2*k + 1];
         }

         // Recursive call 
         Complex[] a0 = InverseDFT(y0);
         Complex[] a1 = InverseDFT(y1);

         // Constructing the final result
         Complex[] a = new Complex[n];
         for (int k = 0; k < n/2; k++) {
            Complex twiddle = Complex.FromPolar(1,-2*k*Math.PI/n);

            a[k]     = (a0[k] + twiddle*a1[k])/2;
            a[k+n/2] = (a0[k] - twiddle*a1[k])/2;
         }

         return a;

      } //InverseFFT

      ///
      /// \copybrief FFT::IDFT::ToString()
      ///
      /// Returns the string <tt>"Recursive"</tt>.
      ///
      public override string ToString()
      {
         return "Recursive";
      }
   }
}

---