rfft.cpp

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

#include "rfft.h"

using namespace spectral;

/// Real FFT constructor.  Creates real symmetric FFT object for sequence of length N.
/// \param N length of transform
/// \throws bad_parameter if N<2

rfft::rfft(int N)
{
  if(N<2)
    throw bad_parameter();
  n=N;
  work=new real[n];
  top=new node(n,1,1,0);
  bot=top->end();
}

/// Real FFT destructor.

rfft::~rfft()
{
  delete [] work;
  delete top;
}

/// Internal real FFT radix object.  
/// Recursively defined double-linked list object contains
/// butterfly rotation seeds for each factor.
/// \param m length of subsequence from factor down
/// \param L accumlation product of previous factors
/// \param parity even/odd position of factor
/// \param that pointer from calling node

rfft::node::node(int m,int L,bool parity,node *that)
{
  prev=that;
  odd=parity;
  r=factor(m);
  m=m/r;
  real angle=-2.0*pi/((real)(L*r));
  omega.re=-2.0*std::sin(0.5*angle)*std::sin(0.5*angle);
  omega.im=std::sin(angle);
  switch(r)
    {
    case 2:
      break;
    case 3:
      break;
    case 4:
      break;
    case 5:
      break;
    case 6:
      break;
    case 8:
      break;
    default:
      R=new complex[r];
      T=new complex[r];
      for(int k=0;k<r;k++)
	{
	  angle=-((real)k)*(2.0*pi/(real)r);
	  R[k].re=std::cos(angle);
	  R[k].im=std::sin(angle);
	}
      break;
    }
  L=L*r;
  if(m>1)
    next=new node(m,L,!parity,this);
  else
    next=0;
}

/// Internal method to that returns end of linked list.

rfft::node *rfft::node::end()
{
  if(next)
    return(next->end());
  else
    return(this);
}

/// Internal factorization routine.

int rfft::node::factor(int m)
{
  if(m%8==0&&m!=16)
    return(8); 
  else if(m%6==0)
    return(6);
  else if(m%5==0)
    return(5);
  else if(m%4==0)
    return(4);
  else if(m%3==0)
    return(3);
  else if(m%2==0)
    return(2);
  else
    for(int f=7;f*f<=m;f+=2)
      if(m%f==0)
	return(f);
  return(m);
}

/// Internal node destructor.

rfft::node::~node()
{
  if(next==0)
    switch(r)
      {
      case 2:
      case 3:
      case 4:
      case 5:
      case 6:
      case 8:
	break;
      default:
	delete [] T;
	delete [] R;
	break;
      }
  else
    delete next;
}

/// Forward transform.  Computes the Fourier coefficients
/// \f[ \hat{a}_k = \sum_{n=0}^{N-1} a_n e^{-2 \pi i k n/N} \f]
/// for \f$ k=0 \ldots \lfloor \frac{N}{2} \rfloor \f$ for a real sequence.
/// The Fourier coefficients are stored in the so-called half complex
/// format, packed into the original real array as
/// \f[ \left( \Re \hat{a}_0, \Re \hat{a}_1, \Im\hat{a}_1, \ldots, \left\{
/// \begin{array}{ll}  \Re \hat{a}_{N/2} & \textrm{$N$ even} \\ 
///  \Re \hat{a}_{(N-1)/2}, \Im \hat{a}_{(N-1)/2} & \textrm{$N$ odd} 
/// \end{array}  \right. \right) \f]
/// \param a real sequence to transform in-place

void rfft::analysis(real *a)
{
  top->analysis(a,work,1,n);
}

/// Backward transform.  Computes the real sequence values
/// \f[ a_n = \sum_{k=0}^{N-1} \hat{a}_k e^{2 \pi i k n/N} \f]
/// for \f$ n=0 \ldots N-1 \f$ using the fact that the original
/// real symmetric sequence has Fourier coefficients with the
/// property \f[ \hat{a}_k = \overline{\hat{a}}_{N-k} \f]
/// from the half-complex Fourier coefficients.
/// \param a sequence to transform in-place

void rfft::synthesis(real *a)
{
  bot->synthesis(a,work,n,1);
}

/// Multiple instance forward transform.
/// \param a array of pointers to sequences to transform
/// \param M number of sequences to transform
/// \param P offset (from the default zero) of first sequence in array 

void rfft::analysis(real **a,int M,int P)
{
  if((M<1)||(P<0))
    throw bad_parameter();
  for(int i=P;i<P+M;i++)
    top->analysis(a[i],work,1,n);
}

/// Multiple instance backward transform.
/// \param a array of pointers to sequences to transform
/// \param M number of sequences to transform
/// \param P offset (from the default zero) of first sequence in array 

void rfft::synthesis(real **a,int M,int P)
{
  if((M<1)||(P<0))
    throw bad_parameter();
  for(int i=P;i<P+M;i++)
    bot->synthesis(a[i],work,n,1);
}

/// Multiple instance forward transform.
/// \param a array of contiguous sequences to transform
/// \param M number of sequences to transform
/// \param stride distance between each sequence (must be at least n)

void rfft::analysis(real *a,int M,int stride)
{
  if((M<1)||(stride<n))
    throw bad_parameter();
  for(int i=0;i<M;i++)
    top->analysis(a+i*stride,work,1,n);
}

/// Multiple instance backward transform.
/// \param a array of contiguous sequences to transform
/// \param M number of sequences to transform
/// \param stride distance between each sequence (must be at least n)

void rfft::synthesis(real *a,int M,int stride)
{
  if((M<1)||(stride<n))
    throw bad_parameter();
  for(int i=0;i<M;i++)
    bot->synthesis(a+i*stride,work,n,1);
}

/// Internal recursive radix-node butterfly forward transform driver.
/// \param in input array
/// \param out output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis(real *in,real *out,int L,int N)
{
  N=N/r;
  switch(r)
    {
    case 2:
      analysis_radix2(in,out,L,N);
      break;
    case 3:
      analysis_radix3(in,out,L,N);
      break;
    case 4:
      analysis_radix4(in,out,L,N);
      break;
    case 5:
      analysis_radix5(in,out,L,N);
      break;
    case 6:
      analysis_radix6(in,out,L,N);
      break;
    case 8:
      analysis_radix8(in,out,L,N);
      break;
    default:
      analysis_radixg(in,out,L,N);
      break;
    }
  L=L*r;
  if(next)
    next->analysis(out,in,L,N);
  else
    if(odd)
      for(int i=0;i<L;i++)
	in[i]=out[i];
}

/// Internal recursive radix-node butterfly backward transform driver.
/// \param in input array
/// \param out output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::synthesis(real *in,real *out,int L,int N)
{
  L=L/r;
  switch(r)
    {
    case 2:
      synthesis_radix2(in,out,L,N);
      break;
    case 3:
      synthesis_radix3(in,out,L,N);
      break;
    case 4:
      synthesis_radix4(in,out,L,N);
      break;
    case 5:
      synthesis_radix5(in,out,L,N);
      break;
    case 6:
      synthesis_radix6(in,out,L,N);
      break;
    case 8:
      synthesis_radix8(in,out,L,N);
      break;
    default:
      synthesis_radixg(in,out,L,N);
      break;
    }
  N=N*r;
  if(prev)
    prev->synthesis(out,in,L,N);
  else
    if(odd)
      for(int i=0;i<N;i++)
	in[i]=out[i];
}

/// Internal radix-2 butterfly forward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis_radix2(real *x,real *y,int L,int N)
{
  complex w,T0,T1;
  real *x0,*x1,*y0,*y1;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      y0=y+j*2*L;
      y1=y+j*2*L+2*L;
      T0.re=x0[0];
      T1.re=x1[0];
      y0[ 0]=T0.re+T1.re;
      y1[-1]=T0.re-T1.re;
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  T0.re=x0[2*k-1];
	  T0.im=x0[2*k  ];
	  T1.re=x1[2*k-1]*w.re-x1[2*k  ]*w.im;
	  T1.im=x1[2*k  ]*w.re+x1[2*k-1]*w.im;
	  y0[ 2*k-1]=T0.re+T1.re;              
	  y0[ 2*k  ]=T0.im+T1.im;
	  y1[-2*k-1]=T0.re-T1.re;
	  y1[-2*k  ]=T1.im-T0.im;
	  T0.re=omega.re*w.re-omega.im*w.im+w.re;
	  T0.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=T0.re;
	  w.im=T0.im;
        }
      if(L%2==0)
        {
	  T0.re=x0[L-1];
	  T1.re=x1[L-1]*w.re;
	  T1.im=x1[L-1]*w.im;
	  y0[L-1]=T0.re+T1.re;
	  y0[L  ]=T1.im;
        }
    }
}

/// Internal radix-2 butterfly backward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::synthesis_radix2(real *y,real *x,int L,int N)
{
  complex w,a,b;
  complex T0,T1;
  real *x0,*x1,*y0,*y1;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      y0=y+j*2*L;
      y1=y+j*2*L+2*L;
      T0.re=y0[ 0];
      T1.re=y1[-1];
      x0[0]=T0.re+T1.re;
      x1[0]=T0.re-T1.re;
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  T0.re= y0[ 2*k-1];
	  T0.im= y0[ 2*k  ];
	  T1.re= y1[-2*k-1];
	  T1.im=-y1[-2*k  ];
	  a.re=T0.re+T1.re;
	  a.im=T0.im+T1.im;
	  b.re=T0.re-T1.re;
	  b.im=T0.im-T1.im;
	  x0[2*k-1]=a.re;
	  x0[2*k  ]=a.im;
	  x1[2*k-1]=b.re*w.re+b.im*w.im;
	  x1[2*k  ]=b.im*w.re-b.re*w.im;
	  a.re=omega.re*w.re-omega.im*w.im+w.re;
	  a.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=a.re;
	  w.im=a.im;
        }
      if(L%2==0)
        {
	  T0.re=y0[L-1]; 
	  T0.im=y0[L  ]; 
	  T1.re= T0.re;
	  T1.im=-T0.im;
	  a.re=T0.re+T1.re;
	  a.im=T0.im+T1.im;
	  b.re=T0.re-T1.re;
	  b.im=T0.im-T1.im;
	  x0[L-1]=a.re;
	  x1[L-1]=b.re*w.re+b.im*w.im;
        }
    }
}

/// Internal radix-3 butterfly forward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis_radix3(real *x,real *y,int L,int N)
{
  complex t0,t1,t2;
  complex w1,w2;
  complex z1,z2;
  real *x0,*x1,*x2,*y0,*y1;
  const real a=-0.5L;
  const real b= 0.5L*root3;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      y0=y+j*3*L;
      y1=y+j*3*L+2*L;
      t0.re=x0[0]; 
      t1.re=x1[0];
      t2.re=x2[0];
      z1.re=t1.re+t2.re;
      z1.im=t2.re-t1.re;
      y0[ 0]=t0.re+z1.re;
      y1[-1]=t0.re+a*z1.re;
      y1[ 0]=      b*z1.im;
      w1.re=omega.re+1.0;
      w1.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  t0.re=x0[2*k-1];
	  t0.im=x0[2*k  ];
	  z1.re=x1[2*k-1];
	  z1.im=x1[2*k  ];
	  z2.re=x2[2*k-1];
	  z2.im=x2[2*k  ];
	  t1.re=z1.re*w1.re-z1.im*w1.im;
	  t1.im=z1.im*w1.re+z1.re*w1.im;
	  t2.re=z2.re*w2.re-z2.im*w2.im;
	  t2.im=z2.im*w2.re+z2.re*w2.im;
	  z1.re=a*(t1.re+t2.re);
	  z1.im=a*(t1.im+t2.im);
	  z2.re=b*(t1.re-t2.re);
	  z2.im=b*(t1.im-t2.im);
	  y0[ 2*k-1]= (t0.re+t1.re+t2.re);
	  y0[ 2*k  ]= (t0.im+t1.im+t2.im);
	  y1[ 2*k-1]= (t0.re+z1.re+z2.im);
	  y1[ 2*k  ]= (t0.im+z1.im-z2.re);
	  y1[-2*k-1]= (t0.re+z1.re-z2.im);
	  y1[-2*k  ]=-(t0.im+z1.im+z2.re);
	  w2.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  w2.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=w2.re;
	  w1.im=w2.im;
        }
      if(L%2==0)
        {
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  t0.re=x0[L-1];
	  t0.im=0.0;
	  t1.re=x1[L-1]*w1.re; 
	  t1.im=x1[L-1]*w1.im;
	  t2.re=x2[L-1]*w2.re; 
	  t2.im=x2[L-1]*w2.im;
	  y0[L-1]=t0.re+t1.re+t2.re; 
	  y0[L  ]=t0.im+t1.im+t2.im; 
	  y1[L-1]=t0.re+a*(t1.re+t2.re)+b*(t1.im-t2.im);
        }
    }
}

/// Internal radix-3 butterfly backward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::synthesis_radix3(real *y,real *x,int L,int N)
{
  complex t0,t1,t2;
  complex w1,w2;
  complex z1,z2;
  complex a1,a2;
  real *x0,*x1,*x2,*y0,*y1;
  const real a=-0.5L;
  const real b= 0.5L*root3;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      y0=y+j*3*L;
      y1=y+j*3*L+2*L;
      t0.re=y0[ 0];
      t1.re=y1[-1];
      t1.im=y1[ 0];
      z1.re=a*t1.re;
      z1.im=b*t1.im;
      x0[0]=t0.re+2.0*t1.re;
      x1[0]=t0.re+2.0*(z1.re-z1.im);
      x2[0]=t0.re+2.0*(z1.re+z1.im);
      w1.re=omega.re+1.0;                   
      w1.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  t0.re= y0[ 2*k-1];
	  t0.im= y0[ 2*k  ];
	  t1.re= y1[ 2*k-1];
	  t1.im= y1[ 2*k  ];
	  t2.re= y1[-2*k-1];
	  t2.im=-y1[-2*k  ];
	  a1.re=a*(t1.re+t2.re);
	  a1.im=a*(t1.im+t2.im);
	  a2.re=b*(t1.re-t2.re);
	  a2.im=b*(t1.im-t2.im);
	  z1.re=t0.re+a1.re-a2.im;
	  z1.im=t0.im+a1.im+a2.re;
	  z2.re=t0.re+a1.re+a2.im;
	  z2.im=t0.im+a1.im-a2.re;
	  x0[2*k-1]=t0.re+t1.re+t2.re;
	  x0[2*k  ]=t0.im+t1.im+t2.im;
	  x1[2*k-1]=z1.re*w1.re+z1.im*w1.im;
	  x1[2*k  ]=z1.im*w1.re-z1.re*w1.im;
	  x2[2*k-1]=z2.re*w2.re+z2.im*w2.im;
	  x2[2*k  ]=z2.im*w2.re-z2.re*w2.im;
	  z1.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  z1.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=z1.re;
	  w1.im=z1.im;
        }      
      if(L%2==0)
        {
	  w2.re=w1.re*w1.re-w1.im*w1.im;
	  w2.im=2.0*w1.re*w1.im;
	  t0.re=y0[L-1]; 
	  t0.im=y0[L  ]; 
	  t1.re=y1[L-1]; 
	  a1.re=a*(t1.re+t0.re);
	  a1.im=b*(t1.re-t0.re);
	  a2.re=a*t0.im;
	  a2.im=b*t0.im;
	  z1.re=t0.re+a1.re-a2.im; 
	  z1.im=t0.im+a1.im-a2.re;