cfft.h

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// spectral toolkit 
// copyright (c) 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//

#ifndef __cfft__
#define __cfft__

#include "spectral.h"

namespace spectral
{
  /// Complex Fast Fourier Transform. Constructs FFT object for a particular length sequence.  There are no restrictions on
  /// sequence length, although best performance is attained when the length is a composite number of small primes.
  /// Uses a transposed variant of the Stockham algorithm in
  /// Swartztrauber, P.N., FFT algorithms for vector computers, <I>Parallel Computing</I>, <B>1</B> (1984), pp. 45-63.
  ///
  /// The following example uses the cfft to compute the derivative of cos(x) and checks the result.
  /// \code
  /// #include <cfft.h>
  /// #include <alloc.h>
  /// #include <iostream>
  /// #include <math.h>
  /// 
  /// using namespace spectral;
  /// 
  /// int main()
  /// {
  ///   int n=10;
  ///   real dx=2.0*pi/(real)n;
  ///   cfft FFT(n);           
  ///   complex *c=alloc<complex>(n);
  ///   for(int i=0;i<n;i++)             
  ///     c[i]=complex(cos(-pi+(real)i*dx),0.0);   
  ///   FFT.analysis(c);             
  ///   for(int i=0;i<n;i++)
  ///     c[i]*=complex(0.0,wavenumber(n,i)/n);
  ///   FFT.synthesis(c);            
  ///   real error=0.0;              
  ///   real l_inf;
  ///   for(int i=0;i<n;i++)
  ///     {
  ///       l_inf=abs(c[i]+sin(-pi+(real)i*dx));
  ///       error=max(error,l_inf);
  ///     }
  ///   std::cout << "l-infinity error=" << error << std::endl;
  ///   return(0);
  /// }
  /// \endcode
  ///
  /// \sa rfft
  
  class cfft
  {
  public:
    cfft(int Length);
    ~cfft();
    void analysis(complex **a,int m,int p=0);
    void synthesis(complex **a,int m,int p=0);
    void analysis(complex *a,int m,int s);
    void synthesis(complex *a,int m,int s);
    void analysis(complex *a);
    void synthesis(complex *a);
  private:
    /// internal recursive radix object for cfft
    class node
    {
    public:
      node(int,int,bool);
      ~node();
      void analysis(complex*,complex*,int,int);
      void synthesis(complex*,complex*,int,int);
    private:
      /// rotation
      complex omega;
      /// temp vector for general radix r
      complex *T;
      /// first column of butterfly rotation for general radix r
      complex *R;
      /// pointer to next node in linked list
      node *next;
      /// factor
      int r;
      /// flag for odd number of factors trick to overwrite on last factor
      bool odd;
      int factor(int);
      void analysis_radix2(complex*,complex*,int,int);
      void analysis_radix3(complex*,complex*,int,int);
      void analysis_radix4(complex*,complex*,int,int);
      void analysis_radix5(complex*,complex*,int,int);
      void analysis_radix6(complex*,complex*,int,int);
      void analysis_radix8(complex*,complex*,int,int);
      void analysis_radixg(complex*,complex*,int,int);
      void synthesis_radix2(complex*,complex*,int,int);
      void synthesis_radix3(complex*,complex*,int,int);
      void synthesis_radix4(complex*,complex*,int,int);
      void synthesis_radix5(complex*,complex*,int,int);
      void synthesis_radix6(complex*,complex*,int,int);
      void synthesis_radix8(complex*,complex*,int,int);
      void synthesis_radixg(complex*,complex*,int,int);
    };
    /// temp vector for ping-pong
    complex *work;
    /// length of sequence
    int n;
    /// top of linked list
    node *top;
  };
}

#endif


// Local Variables:
// mode:C++
// End: