cfft.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 "cfft.h"

using namespace spectral;

/// Complex FFT constructor.  Creates complex FFT object for sequence of length N.
/// \param N length of transform
/// \throws bad_parameter if N<2

cfft::cfft(int N)
{
  if(N<2)
    throw bad_parameter();
  n=N;
  work=new complex[n];
  top=new node(n,1,1);
}

/// Complex FFT destructor.

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

/// Internal complex FFT radix object.  Recursively defined 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

cfft::node::node(int m,int L,bool parity)
{
  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:
    case 3:
    case 4:
    case 5:
    case 6:
    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);
  else
    next=0;
}

/// Internal factorization routine.

int cfft::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 radix object destructor.

cfft::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;
}

/// 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 cfft::analysis(complex **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 cfft::synthesis(complex **a,int M,int P)
{
  if((M<1)||(P<0))
    throw bad_parameter();
  for(int i=P;i<P+M;i++)
    top->synthesis(a[i],work,1,n);
}

/// 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 cfft::analysis(complex *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 cfft::synthesis(complex *a,int M,int stride)
{
  if((M<1)||(stride<n))
    throw bad_parameter();
  for(int i=0;i<M;i++)
    top->synthesis(a+i*stride,work,1,n);
}

/// 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 N-1 \f$ for a sequence of complex numbers.
/// The wavenumbers are stored in the standard manner as
/// \f[ k=\{0,1,\ldots ,\lfloor n/2 \rfloor -1,- \lfloor n/2 \rfloor ,\ldots ,-1\}. \f]
/// \param a sequence to transform in-place
/// \sa wavenumber

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

/// Backward transform.  Computes the 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$ from the Fourier coefficients.
/// \param a sequence to transform in-place

void cfft::synthesis(complex *a)
{
  top->synthesis(a,work,1,n);
}

/// 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 cfft::node::analysis(complex *in,complex *out,int L,int N)
{
  N=N/r;
  if(!next&&odd)
    out=in;
  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);
}

/// 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 cfft::node::synthesis(complex *in,complex *out,int L,int N)
{
  N=N/r;
  if(!next&&odd)
    out=in;
  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;
    }
  L=L*r;
  if(next)
    next->synthesis(out,in,L,N);
}

/// 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 cfft::node::analysis_radix2(complex *x,complex *y,int L,int N)
{
  complex w;
  complex t0,t1;
  complex *x0,*x1;
  complex *y0,*y1;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+N*L+j*L;
      y0=y+j*2*L;
      y1=y+L+j*2*L;
      w.re=1.0;
      w.im=0.0;
      for(int i=0;i<L;i++)
        {
	  t0.re=x0[i].re;
	  t0.im=x0[i].im;
	  t1.re=w.re*x1[i].re-w.im*x1[i].im;
	  t1.im=w.re*x1[i].im+w.im*x1[i].re;
	  y0[i].re=t0.re+t1.re;
	  y0[i].im=t0.im+t1.im;
	  y1[i].re=t0.re-t1.re;
	  y1[i].im=t0.im-t1.im;
	  t1.re=omega.re*w.re-omega.im*w.im+w.re;
	  t1.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=t1.re;
	  w.im=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 cfft::node::synthesis_radix2(complex *x,complex *y,int L,int N)
{
  complex w;
  complex t0,t1;
  complex *x0,*x1;
  complex *y0,*y1;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+N*L+j*L;
      y0=y+j*2*L;
      y1=y+L+j*2*L;
      w.re=1.0;
      w.im=0.0;
      for(int i=0;i<L;i++)
        {
	  t0.re=x0[i].re;
	  t0.im=x0[i].im;
	  t1.re=w.re*x1[i].re-w.im*x1[i].im;
	  t1.im=w.re*x1[i].im+w.im*x1[i].re;
	  y0[i].re=t0.re+t1.re;
	  y0[i].im=t0.im+t1.im;
	  y1[i].re=t0.re-t1.re;
	  y1[i].im=t0.im-t1.im;
	  t1.re=omega.re*w.re+omega.im*w.im+w.re;
	  t1.im=omega.re*w.im-omega.im*w.re+w.im; 
	  w.re=t1.re;
	  w.im=t1.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 cfft::node::analysis_radix3(complex *x,complex *y,int L,int N)
{
  complex w,w2;
  complex t0,t1,t2;
  complex z0,z1;
  const real a=-0.5L;
  const real b= 0.5L*root3;
  complex *x0,*x1,*x2;
  complex *y0,*y1,*y2;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x0+N*L;
      x2=x1+N*L;
      y0=y+3*L*j;
      y1=y0+L;
      y2=y1+L;
      w.re=1.0;
      w.im=0.0;
      for(int i=0;i<L;i++)
        {
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  t0.re=x0[i].re;
	  t0.im=x0[i].im;
	  t1.re=w.re*x1[i].re-w.im*x1[i].im;
	  t1.im=w.re*x1[i].im+w.im*x1[i].re;
	  t2.re=w2.re*x2[i].re-w2.im*x2[i].im;
	  t2.im=w2.re*x2[i].im+w2.im*x2[i].re;
	  z0.re=t0.re+a*(t1.re+t2.re);
	  z0.im=      b*(t1.im-t2.im);
	  z1.re=t0.im+a*(t1.im+t2.im);
	  z1.im=      b*(t2.re-t1.re);
	  y0[i].re=t0.re+t1.re+t2.re;
	  y0[i].im=t0.im+t1.im+t2.im;
	  y1[i].re=z0.re+z0.im;
	  y1[i].im=z1.re+z1.im;
	  y2[i].re=z0.re-z0.im;
	  y2[i].im=z1.re-z1.im;
	  t1.re=omega.re*w.re-omega.im*w.im+w.re;
	  t1.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=t1.re;
	  w.im=t1.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 cfft::node::synthesis_radix3(complex *x,complex *y,int L,int N)
{
  complex w,w2;
  complex t0,t1,t2;
  complex z0,z1;
  const real a=-0.5L;
  const real b= 0.5L*root3;
  complex *x0,*x1,*x2;
  complex *y0,*y1,*y2;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x0+N*L;
      x2=x1+N*L;
      y0=y+3*L*j;
      y1=y0+L;
      y2=y1+L;
      w.re=1.0;
      w.im=0.0;
      for(int i=0;i<L;i++)
        {
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  t0.re=x0[i].re;
	  t0.im=x0[i].im;
	  t1.re=w.re*x1[i].re-w.im*x1[i].im;
	  t1.im=w.re*x1[i].im+w.im*x1[i].re;
	  t2.re=w2.re*x2[i].re-w2.im*x2[i].im;
	  t2.im=w2.re*x2[i].im+w2.im*x2[i].re;
	  z0.re=t0.re+a*(t1.re+t2.re);
	  z0.im=      b*(t2.im-t1.im);
	  z1.re=t0.im+a*(t1.im+t2.im);
	  z1.im=      b*(t1.re-t2.re);
	  y0[i].re=t0.re+t1.re+t2.re;
	  y0[i].im=t0.im+t1.im+t2.im;
	  y1[i].re=z0.re+z0.im;
	  y1[i].im=z1.re+z1.im;
	  y2[i].re=z0.re-z0.im;
	  y2[i].im=z1.re-z1.im;
	  t1.re=omega.re*w.re+omega.im*w.im+w.re;
	  t1.im=omega.re*w.im-omega.im*w.re+w.im; 
	  w.re=t1.re;
	  w.im=t1.im;
        }
    }
}  

/// Internal radix-4 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 cfft::node::analysis_radix4(complex *x,complex *y,int L,int N)
{
  complex w1,w2,w3;
  complex t0,t1,t2,t3;
  complex s0,s1,s2,s3;
  complex *x0,*x1,*x2,*x3;
  complex *y0,*y1,*y2,*y3;
    
  for(int j=0;j<N;j++)
    { 
      x0=x+j*L;
      x1=x0+L*N;
      x2=x1+L*N;
      x3=x2+L*N;
      y0=y+4*L*j;
      y1=y0+L;
      y2=y1+L;
      y3=y2+L;
      t0.re=x0[0].re;
      t0.im=x0[0].im;
      t1.re=x1[0].re;
      t1.im=x1[0].im;
      t2.re=x2[0].re;
      t2.im=x2[0].im;
      t3.re=x3[0].re;
      t3.im=x3[0].im;
      s0.re=t0.re+t2.re;
      s0.im=t0.im+t2.im;
      s1.re=t0.re-t2.re;
      s1.im=t0.im-t2.im;
      s2.re=t1.re+t3.re;
      s2.im=t1.im+t3.im;
      s3.re=t1.re-t3.re;
      s3.im=t1.im-t3.im;
      y0[0].re = s0.re+s2.re;
      y0[0].im = s0.im+s2.im;
      y1[0].re = s1.re+s3.im;
      y1[0].im = s1.im-s3.re;
      y2[0].re = s0.re-s2.re;
      y2[0].im = s0.im-s2.im;
      y3[0].re = s1.re-s3.im;
      y3[0].im = s1.im+s3.re;
      w1.re=omega.re+1.0;
      w1.im=omega.im;
      for(int i=1;i<L;i++)
        {
	  w2.re = w1.re*w1.re-w1.im*w1.im;
	  w2.im = 2.0*w1.re*w1.im;
	  w3.re = w2.re*w1.re-w2.im*w1.im;
	  w3.im = w2.re*w1.im+w2.im*w1.re; 
	  t0.re = x0[i].re;
	  t0.im = x0[i].im;
	  t1.re = w1.re*x1[i].re-w1.im*x1[i].im;
	  t1.im = w1.re*x1[i].im+w1.im*x1[i].re;
	  t2.re = w2.re*x2[i].re-w2.im*x2[i].im;
	  t2.im = w2.re*x2[i].im+w2.im*x2[i].re;
	  t3.re = w3.re*x3[i].re-w3.im*x3[i].im;
	  t3.im = w3.re*x3[i].im+w3.im*x3[i].re;
	  s0.re=t0.re+t2.re;
	  s0.im=t0.im+t2.im;
	  s1.re=t0.re-t2.re;
	  s1.im=t0.im-t2.im;
	  s2.re=t1.re+t3.re;
	  s2.im=t1.im+t3.im;
	  s3.re=t1.re-t3.re;
	  s3.im=t1.im-t3.im;
	  y0[i].re = s0.re+s2.re;
	  y0[i].im = s0.im+s2.im;
	  y1[i].re = s1.re+s3.im;
	  y1[i].im = s1.im-s3.re;
	  y2[i].re = s0.re-s2.re;
	  y2[i].im = s0.im-s2.im;
	  y3[i].re = s1.re-s3.im;
	  y3[i].im = s1.im+s3.re;
	  t1.re=omega.re*w1.re-w1.im*omega.im+w1.re;
	  t1.im=omega.re*w1.im+w1.re*omega.im+w1.im;
	  w1.re=t1.re;
	  w1.im=t1.im;
        }
    }
}  

/// Internal radix-4 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 cfft::node::synthesis_radix4(complex *x,complex *y,int L,int N)
{
  complex w1,w2,w3;
  complex t0,t1,t2,t3;
  complex s0,s1,s2,s3;
  complex *x0,*x1,*x2,*x3;
  complex *y0,*y1,*y2,*y3;
    
  for(int j=0;j<N;j++)
    { 
      x0=x+j*L;
      x1=x0+L*N;