//// 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<2cfft::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 factorcfft::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 wavenumbervoid 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-placevoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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;