//// 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<2rfft::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 noderfft::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-placevoid 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-placevoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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 downvoid 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;