//// 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: