rfft.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 __rfft__
#define __rfft__

#include "spectral.h"

namespace spectral
{

  /// Real symmetric Fast Fourier Transform. FFT object for real symmetric sequences.
  /// There are no restrictions on sequence length, although best performance is attained when the length is a 
  /// composite number of small primes.  
  /// Uses a direct real symmetric transform based on Edson's method in 
  /// Swartztrauber, P.N., Symmetric FFTs, <I>Mathematics of Computation</I> <B>47</B> (1986), pp. 323-346,
  /// applied to the same transposed Stockham variant used in cfft.
  ///
  /// The following example uses the real FFT to compute the derivative of cos(x).  Note the method used in
  /// multiplying the spectral coefficients by ik for the half-complex packed array.
  /// \code
  /// #include <rfft.h>
  /// #include <alloc.h>
  /// #include <iostream>
  /// #include <math.h>
  ///
  /// using namespace spectral;
  /// 
  /// int main()
  /// {
  ///   int n=100;
  ///   real dx=2.0*pi/(real)n;
  ///   rfft FFT(n);           
  ///   real *c=alloc<real>(n);
  ///   for(int i=0;i<n;i++)             
  ///     c[i]=cos(-pi+(real)i*dx);   
  ///   FFT.analysis(c);
  ///   for(int k=0;k<=n/2;k++)
  ///     {
  ///       if(k==0)
  ///         {
  ///           c[0]=0.0;
  ///         }
  ///       else if((k==n/2)&&(n%2==0))
  ///         {
  ///           c[n-1]=0.0;
  ///         }
  ///       else
  ///         {
  ///           real cr=c[2*k-1];
  ///           real ci=c[2*k];
  ///           c[2*k-1]=-ci*(real)k/(real)n;
  ///           c[2*k]  = cr*(real)k/(real)n;
  ///         }
  ///      }
  ///   FFT.synthesis(c);            
  ///   real error=0.0;              
  ///   real l_inf;
  ///   for(int i=0;i<n;i++)
  ///     {
  ///       l_inf=fabs(c[i]+sin(-pi+(real)i*dx));
  ///       error=max(error,l_inf);
  ///     }
  ///   std::cout << "l-infinity error=" << error << std::endl;
  ///   return(0);
  /// }
  /// \endcode
  ///
  /// \sa cfft

  class rfft
  {
  public:
    rfft(int Length);
    ~rfft();
    void analysis(real **a,int m,int p=0);
    void synthesis(real **a,int m,int p=0);
    void analysis(real *a,int m,int s);
    void synthesis(real *a,int m,int s);
    void analysis(real *a);
    void synthesis(real *a);
  private:
    /// internal recursive radix object for rfft        
    class node
    {
    public:
      node(int,int,bool,node*);
      ~node();
      node *end();
      void analysis(real*,real*,int,int);
      void synthesis(real*,real*,int,int);
    private:
      /// rotation
      complex omega;
      /// temp vector for general radix r
      complex *T;
      /// first column of butterfuly rotation for general radix r
      complex *R;
      /// pointer to forward node in linked list
      node *next;
      /// pointer to backward node in linked list
      node *prev;
      /// radix
      int r;
      /// flag for copy on last factor
      bool odd;
      int factor(int);
      void analysis_radix2(real*,real*,int,int);
      void analysis_radix3(real*,real*,int,int);
      void analysis_radix4(real*,real*,int,int);
      void analysis_radix5(real*,real*,int,int);
      void analysis_radix6(real*,real*,int,int);
      void analysis_radix8(real*,real*,int,int);
      void analysis_radixg(real*,real*,int,int);
      void synthesis_radix2(real*,real*,int,int);
      void synthesis_radix3(real*,real*,int,int);
      void synthesis_radix4(real*,real*,int,int);
      void synthesis_radix5(real*,real*,int,int);
      void synthesis_radix6(real*,real*,int,int);
      void synthesis_radix8(real*,real*,int,int);
      void synthesis_radixg(real*,real*,int,int);
    };
    /// temp vector for ping-ping
    real *work;
    /// top of linked list
    node *top;
    /// bottom of linked list
    node *bot;
    /// length of sequence
    int n;
        
  };
}

#endif

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