rfft.cpp

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

	  y3[L-1] = (t6.re+t7.re);
	  y3[L  ] = (t6.im+t7.im);
        } 
    }
}

/// Internal radix-8 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_radix8(real *y,real *x,int L,int N)
{
  complex w1,w2,w3,w4,w5,w6,w7;
  complex s0,s1,s2,s3,s4,s5,s6,s7;
  complex t0,t1,t2,t3,t4,t5,t6,t7;
  real *x0,*x1,*x2,*x3,*x4,*x5,*x6,*x7;
  real *y0,*y1,*y2,*y3,*y4;
  const real c=0.5L*root2;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x0+N*L;
      x2=x1+N*L;
      x3=x2+N*L;
      x4=x3+N*L;  
      x5=x4+N*L;
      x6=x5+N*L;
      x7=x6+N*L;  
      y0=y+j*L*8;
      y1=y+2*L+j*L*8;
      y2=y+4*L+j*L*8;
      y3=y+6*L+j*L*8;
      y4=y+8*L+j*L*8;
      t0.re=y0[ 0];
      t1.re=y1[-1];
      t1.im=y1[ 0];
      t2.re=y2[-1];
      t2.im=y2[ 0];
      t3.re=y3[-1];
      t3.im=y3[ 0];
      t4.re=y4[-1];
      s0.re=t0.re+2.0*(t1.re+t2.re+t3.re)+t4.re;
      s1.re=t0.re+2.0*(c*(t1.re-t1.im)-t2.im-c*(t3.re+t3.im))-t4.re;
      s2.re=t0.re+2.0*(-t1.im-t2.re+t3.im)+t4.re;
      s3.re=t0.re+2.0*(c*(-t1.re-t1.im)+t2.im+c*(t3.re-t3.im))-t4.re;
      s4.re=t0.re+2.0*(-t1.re+t2.re-t3.re)+t4.re;
      s5.re=t0.re+2.0*(c*(-t1.re+t1.im)-t2.im+c*(t3.re+t3.im))-t4.re;
      s6.re=t0.re+2.0*(t1.im-t2.re-t3.im)+t4.re;
      s7.re=t0.re+2.0*(c*(t1.re+t1.im)+t2.im+c*(-t3.re+t3.im))-t4.re;
      x0[0]=s0.re;
      x1[0]=s1.re;
      x2[0]=s2.re;
      x3[0]=s3.re;
      x4[0]=s4.re;
      x5[0]=s5.re;
      x6[0]=s6.re;
      x7[0]=s7.re;
      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;
	  w3.re=w2.re*w1.re-w2.im*w1.im;
	  w3.im=w2.re*w1.im+w2.im*w1.re; 
	  w4.re=w3.re*w1.re-w3.im*w1.im;
	  w4.im=w3.re*w1.im+w3.im*w1.re; 
	  w5.re=w4.re*w1.re-w4.im*w1.im;
	  w5.im=w4.re*w1.im+w4.im*w1.re; 
	  w6.re=w5.re*w1.re-w5.im*w1.im;
	  w6.im=w5.re*w1.im+w5.im*w1.re; 
	  w7.re=w6.re*w1.re-w6.im*w1.im;
	  w7.im=w6.re*w1.im+w6.im*w1.re; 
	  t0.re= y0[ 2*k-1];
	  t0.im= y0[ 2*k  ];
	  t1.re= y1[ 2*k-1];
	  t1.im= y1[ 2*k  ];
	  t2.re= y2[ 2*k-1];
	  t2.im= y2[ 2*k  ];
	  t3.re= y3[ 2*k-1];
	  t3.im= y3[ 2*k  ];
	  t4.re= y1[-2*k-1];
	  t4.im=-y1[-2*k  ];
	  t5.re= y2[-2*k-1];
	  t5.im=-y2[-2*k  ];
	  t6.re= y3[-2*k-1];
	  t6.im=-y3[-2*k  ];
	  t7.re= y4[-2*k-1];
	  t7.im=-y4[-2*k  ];
	  s0.re=  ( t0.re+t2.re+t5.re+t7.re);
	  s0.im=  ( t0.im+t2.im+t5.im+t7.im);
	  s1.re=  ( t1.re+t3.re+t4.re+t6.re);
	  s1.im=  ( t1.im+t3.im+t4.im+t6.im);
	  s2.re=  ( t0.re-t2.im+t5.im-t7.re);
	  s2.im=  ( t0.im+t2.re-t5.re-t7.im);
	  s3.re=c*( t1.re-t1.im-t3.re-t3.im+t4.re+t4.im-t6.re+t6.im);
	  s3.im=c*( t1.im+t1.re-t3.im+t3.re+t4.im-t4.re-t6.im-t6.re);
	  s4.re=  ( t0.re-t2.re-t5.re+t7.re);
	  s4.im=  ( t0.im-t2.im-t5.im+t7.im);
	  s5.re=  (-t1.im+t3.im+t4.im-t6.im);
	  s5.im=  ( t1.re-t3.re-t4.re+t6.re);
	  s6.re=  ( t0.re+t2.im-t5.im-t7.re);
	  s6.im=  ( t0.im-t2.re+t5.re-t7.im);
	  s7.re=c*(-t1.re-t1.im+t3.re-t3.im-t4.re+t4.im+t6.re+t6.im);
	  s7.im=c*(-t1.im+t1.re+t3.im+t3.re-t4.im-t4.re+t6.im-t6.re);
	  t0.re=s0.re+s1.re;
	  t0.im=s0.im+s1.im;
	  t1.re=s2.re+s3.re;
	  t1.im=s2.im+s3.im;
	  t2.re=s4.re+s5.re;
	  t2.im=s4.im+s5.im;
	  t3.re=s6.re+s7.re;
	  t3.im=s6.im+s7.im;
	  t4.re=s0.re-s1.re;
	  t4.im=s0.im-s1.im;
	  t5.re=s2.re-s3.re;
	  t5.im=s2.im-s3.im;
	  t6.re=s4.re-s5.re;
	  t6.im=s4.im-s5.im;
	  t7.re=s6.re-s7.re;
	  t7.im=s6.im-s7.im;
	  x0[2*k-1]=t0.re;
	  x0[2*k  ]=t0.im;
	  x1[2*k-1]=t1.re*w1.re+t1.im*w1.im;
	  x1[2*k  ]=t1.im*w1.re-t1.re*w1.im;
	  x2[2*k-1]=t2.re*w2.re+t2.im*w2.im;
	  x2[2*k  ]=t2.im*w2.re-t2.re*w2.im;
	  x3[2*k-1]=t3.re*w3.re+t3.im*w3.im;
	  x3[2*k  ]=t3.im*w3.re-t3.re*w3.im;
	  x4[2*k-1]=t4.re*w4.re+t4.im*w4.im;
	  x4[2*k  ]=t4.im*w4.re-t4.re*w4.im;
	  x5[2*k-1]=t5.re*w5.re+t5.im*w5.im;
	  x5[2*k  ]=t5.im*w5.re-t5.re*w5.im;
	  x6[2*k-1]=t6.re*w6.re+t6.im*w6.im;
	  x6[2*k  ]=t6.im*w6.re-t6.re*w6.im;
	  x7[2*k-1]=t7.re*w7.re+t7.im*w7.im;
	  x7[2*k  ]=t7.im*w7.re-t7.re*w7.im;
	  t0.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  t0.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=t0.re;
	  w1.im=t0.im;
        }
      if(L%2==0)
        {
	  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; 
	  w4.re=w3.re*w1.re-w3.im*w1.im;
	  w4.im=w3.re*w1.im+w3.im*w1.re; 
	  w5.re=w4.re*w1.re-w4.im*w1.im;
	  w5.im=w4.re*w1.im+w4.im*w1.re; 
	  w6.re=w5.re*w1.re-w5.im*w1.im;
	  w6.im=w5.re*w1.im+w5.im*w1.re; 
	  w7.re=w6.re*w1.re-w6.im*w1.im;
	  w7.im=w6.re*w1.im+w6.im*w1.re; 
	  t0.re=y0[L-1];
	  t0.im=y0[L  ];
	  t1.re=y1[L-1];
	  t1.im=y1[L  ];
	  t2.re=y2[L-1];
	  t2.im=y2[L  ];
	  t3.re=y3[L-1];
	  t3.im=y3[L  ];
	  t4.re= t3.re;
	  t4.im=-t3.im;
	  t5.re= t2.re;
	  t5.im=-t2.im;
	  t6.re= t1.re;
	  t6.im=-t1.im;            
	  t7.re= t0.re;
	  t7.im=-t0.im;  
	  s0.re= t0.re+t2.re+t4.re+t6.re;
	  s0.im= t0.im+t2.im+t4.im+t6.im;
	  s1.re= t1.re+t3.re+t5.re+t7.re;
	  s1.im= t1.im+t3.im+t5.im+t7.im;
	  s2.re= t0.re-t2.im-t4.re+t6.im;
	  s2.im= t0.im+t2.re-t4.im-t6.re;
	  s3.re=c*(t1.re-t1.im-t3.re-t3.im-t5.re+t5.im+t7.re+t7.im);
	  s3.im=c*(t1.im+t1.re-t3.im+t3.re-t5.im-t5.re+t7.im-t7.re);
	  s4.re= t0.re-t2.re+t4.re-t6.re;
	  s4.im= t0.im-t2.im+t4.im-t6.im;
	  s5.re=-t1.im+t3.im-t5.im+t7.im;
	  s5.im= t1.re-t3.re+t5.re-t7.re;
	  s6.re= t0.re+t2.im-t4.re-t6.im;
	  s6.im= t0.im-t2.re-t4.im+t6.re;
	  s7.re=c*(-t1.re-t1.im+t3.re-t3.im+t5.re+t5.im-t7.re+t7.im);
	  s7.im=c*(-t1.im+t1.re+t3.im+t3.re+t5.im-t5.re-t7.im-t7.re);
	  t0.re=s0.re+s1.re;
	  t0.im=s0.im+s1.im;
	  t1.re=s2.re+s3.re;
	  t1.im=s2.im+s3.im;
	  t2.re=s4.re+s5.re;
	  t2.im=s4.im+s5.im;
	  t3.re=s6.re+s7.re;
	  t3.im=s6.im+s7.im;
	  t4.re=s0.re-s1.re;
	  t4.im=s0.im-s1.im;
	  t5.re=s2.re-s3.re;
	  t5.im=s2.im-s3.im;
	  t6.re=s4.re-s5.re;
	  t6.im=s4.im-s5.im;
	  t7.re=s6.re-s7.re;
	  t7.im=s6.im-s7.im;
	  x0[L-1]=t0.re;
	  x1[L-1]=t1.re*w1.re+t1.im*w1.im;
	  x2[L-1]=t2.re*w2.re+t2.im*w2.im;
	  x3[L-1]=t3.re*w3.re+t3.im*w3.im;  
	  x4[L-1]=t4.re*w4.re+t4.im*w4.im;  
	  x5[L-1]=t5.re*w5.re+t5.im*w5.im;  
	  x6[L-1]=t6.re*w6.re+t6.im*w6.im;     
	  x7[L-1]=t7.re*w7.re+t7.im*w7.im;  
        }
    }
}

/// Internal general radix 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_radixg(real *x,real *y,int L,int N)
{
  complex w,z,a,b;
  complex c,d,Z;
  real s,t;
  int l,k,j,q,K;
    
  for(j=0;j<N;j++)
    {
      for(q=0;q<r;q++)
        {
	  T[q].re=x[q*N*L+j*L];
        }
      s=T[0].re;           
      for(q=1;q<r;q++)
        {
	  s=s+T[q].re;
        }
      y[0+j*r*L]=s;
      for(l=1;l<(r+1)/2;l++)
        {
	  z.re=T[0].re;
	  z.im=0.0;
	  w.re=R[l].re;
	  w.im=R[l].im;
	  for(q=1;q<r;q++)
            {
	      z.re=z.re+T[q].re*w.re;
	      z.im=z.im+T[q].re*w.im;
	      a.re=w.re*R[l].re-w.im*R[l].im;
	      a.im=w.re*R[l].im+w.im*R[l].re;
	      w.re=a.re;
	      w.im=a.im;
            }
	  K=l*L;
	  y[2*K-1+j*r*L]=z.re;
	  y[2*K  +j*r*L]=z.im;
        }
      if(r%2==0)
        {
	  t=T[0].re-T[1].re;
	  for(q=2;q<r-1;q+=2)
            {
	      t=t+T[q].re-T[q+1].re;
            }
	  K=(r/2)*L;
	  y[2*K-1+j*r*L]=t;
        }
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(k=1;k<(L+1)/2;k++)
        {
	  z.re=1.0;
	  z.im=0.0;
	  for(q=0;q<r;q++)
            {
	      b.re=x[2*k-1+q*N*L+j*L];
	      b.im=x[2*k  +q*N*L+j*L];
	      T[q].re=b.re*z.re-b.im*z.im;
	      T[q].im=b.im*z.re+b.re*z.im;
	      a.re=z.re*w.re-z.im*w.im;
	      a.im=z.re*w.im+z.im*w.re;
	      z.re=a.re;
	      z.im=a.im;
            } 
	  z.re=T[0].re;              
	  z.im=T[0].im;
	  for(q=1;q<r;q++)
            {
	      z.re=z.re+T[q].re;
	      z.im=z.im+T[q].im;
            }
	  K=k;
	  y[2*K-1+j*r*L]=z.re;
	  y[2*K  +j*r*L]=z.im;
	  for(l=1;l<(r+1)/2;l++)
            {
	      z.re=z.im=0.0;
	      Z.re=Z.im=0.0;
	      a.re=1.0;
	      a.im=0.0;
	      for(q=0;q<r;q++)
                {
		  c.re=T[q].re*a.re;
		  c.im=T[q].re*a.im;
		  d.re=T[q].im*a.re;
		  d.im=T[q].im*a.im;
		  z.re+= c.re-d.im;
		  z.im+= c.im+d.re;
		  Z.re+= c.re+d.im;
		  Z.im+=-c.im+d.re;
		  b.re=a.re*R[l].re-a.im*R[l].im;
		  b.im=a.re*R[l].im+a.im*R[l].re;
		  a.re=b.re;
		  a.im=b.im;
                }
	      K=k+l*L;
	      y[2*K-1+j*r*L] = z.re;
	      y[2*K  +j*r*L] = z.im;
	      K=L-k+(l-1)*L;
	      y[2*K-1+j*r*L] = Z.re;
	      y[2*K  +j*r*L] =-Z.im;
            }
	  if(r%2==0)
            {
	      Z.re=0.0;
	      Z.im=0.0;
	      for(q=0;q<r-1;q+=2)
                {
		  Z.re=Z.re+T[q].re-T[q+1].re;
		  Z.im=Z.im+T[q].im-T[q+1].im;
                }
	      K=L-k+(r/2-1)*L;
	      y[2*K-1+j*r*L] = Z.re;
	      y[2*K  +j*r*L] =-Z.im;
            }
            
	  z.re=omega.re*w.re-omega.im*w.im+w.re;
	  z.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=z.re;
	  w.im=z.im;
        }
      if(L%2==0)
        {
	  z.re=1.0;
	  z.im=0.0;
	  for(q=0;q<r;q++)
            {
	      T[q].re=x[L-1+q*N*L+j*L]*z.re;
	      T[q].im=x[L-1+q*N*L+j*L]*z.im;
	      a.re=z.re*w.re-z.im*w.im;
	      a.im=z.re*w.im+z.im*w.re;
	      z.re=a.re;
	      z.im=a.im;
            } 
	  z.re=T[0].re;
	  z.im=T[0].im;
	  for(q=1;q<r;q++)
            {
	      z.re=z.re+T[q].re;
	      z.im=z.im+T[q].im;
            }
	  K=L/2;
	  y[2*K-1+j*r*L]=z.re;
	  y[2*K  +j*r*L]=z.im;
	  for(l=1;l<r/2;l++)
            {
	      z.re=z.im=0.0;
	      a.re=1.0;
	      a.im=0.0;
	      for(q=0;q<r;q++)
                {
		  z.re=z.re+T[q].re*a.re-T[q].im*a.im;
		  z.im=z.im+T[q].re*a.im+T[q].im*a.re;
		  b.re=a.re*R[l].re-a.im*R[l].im;
		  b.im=a.re*R[l].im+a.im*R[l].re;
		  a.re=b.re;
		  a.im=b.im;
                }
	      K=L/2+l*L;
	      y[2*K-1+j*r*L] = z.re;
	      y[2*K  +j*r*L] = z.im;
            }
	  if(r%2==1)
            {
	      s=0.0;
	      a.re=1.0;
	      a.im=0.0;
	      l=(r-1)/2;
	      for(q=0;q<r;q++)
                {
		  s=s+T[q].re*a.re-T[q].im*a.im;
		  b.re=a.re*R[l].re-a.im*R[l].im;
		  b.im=a.re*R[l].im+a.im*R[l].re;
		  a.re=b.re;
		  a.im=b.im;
                }
	      K=L/2+((r-1)/2)*L;
	      y[2*K-1+j*r*L] = s;
            }
        }
    }
}

/// Internal general radix 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_radixg(real *y,real *x,int L,int N)
{
  complex w,z,a,b,W;
  real s;
  int l,k,j,q;
    
  for(j=0;j<N;j++)
    {
      T[0].re=y[0+j*r*L];
      T[0].im=0.0;
      for(l=1;l<(r+1)/2;l++)
        {
	  T[l].re=y[2*l*L-1+j*r*L];
	  T[l].im=y[2*l*L  +j*r*L];
        }
      if(r%2==0)
        {
	  T[r/2].re=y[r*L-1+j*r*L];
	  T[r/2].im=0.0;
        }
      s=T[0].re;
      for(q=1;q<(r+1)/2;q++)
	s=s+2.0*T[q].re;
      if(r%2==0)
	s=s+T[r/2].re;
      x[j*L]=s;
      for(q=1;q<r;q++)
        {
	  s=T[0].re;
	  w.re=R[q].re;
	  w.im=R[q].im;
	  for(l=1;l<(r+1)/2;l++)
            {
	      s=s+2.0*(w.re*T[l].re+w.im*T[l].im);
	      a.re=w.re*R[q].re-w.im*R[q].im;
	      a.im=w.re*R[q].im+w.im*R[q].re;
	      w.re=a.re;
	      w.im=a.im;
            }
	  if(r%2==0)
	    if(q%2==0)
	      s=s+T[r/2].re;
	    else
	      s=s-T[r/2].re;
	  x[q*N*L+j*L]=s;
        }
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(k=1;k<(L+1)/2;k++)
        {
	  T[0].re=y[2*k-1+j*r*L];
	  T[0].im=y[2*k  +j*r*L];
	  for(l=1;l<(r+1)/2;l++)
            {
	      T[l].re=y[2*(k+l*L)-1+j*r*L];
	      T[l].im=y[2*(k+l*L)  +j*r*L];
	      T[l+(r+1)/2-1].re= y[2*(L-k+(l-1)*L)-1+j*r*L];
	      T[l+(r+1)/2-1].im=-y[2*(L-k+(l-1)*L)  +j*r*L];
            }
	  if(r%2==0)
            {
	      T[r-1].re= y[2*(L-k+(r/2-1)*L)-1+j*r*L];
	      T[r-1].im=-y[2*(L-k+(r/2-1)*L)  +j*r*L];
            }
	  W.re=1.0;
	  W.im=0.0;
	  for(q=0;q<r;q++)
            {
	      z.re=T[0].re;
	      z.im=T[0].im;
	      a.re=R[q].re;
	      a.im=R[q].im;
	      for(l=1;l<(r+1)/2;l++)
                {
		  z.re=z.re+T[l].re*a.re+T[l].im*a.im;
		  z.im=z.im+T[l].im*a.re-T[l].re*a.im;
		  b.re=a.re*R[q].re-a.im*R[q].im;
		  b.im=a.im*R[q].re+a.re*R[q].im;
		  a.re=b.re;
		  a.im=b.im;
                }
	      a.re=R[q].re;
	      a.im=R[q].im;
	      for(l=(r+1)/2;l<r;l++)
                {
		  z.re=z.re+T[l].re*a.re-T[l].im*a.im;
		  z.im=z.im+T[l].im*a.re+T[l].re*a.im;
		  b.re=a.re*R[q].re-a.im*R[q].im;
		  b.im=a.im*R[q].re+a.re*R[q].im;
		  a.re=b.re;
		  a.im=b.im;
                }
	      x[2*k-1+q*N*L+j*L]=z.re*W.re+z.im*W.im;
	      x[2*k  +q*N*L+j*L]=z.im*W.re-z.re*W.im;
	      z.re=W.re*w.re-W.im*w.im;
	      z.im=W.im*w.re+W.re*w.im;
	      W.re=z.re;
	      W.im=z.im;
            }
	  z.re=omega.re*w.re-omega.im*w.im+w.re;
	  z.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=z.re;
	  w.im=z.im;
        }
      if(L%2==0)
        {
	  for(l=0;l<r/2;l++)
            {
	      T[l].re=y[L+2*l*L-1+j*r*L];
	      T[l].im=y[L+2*l*L  +j*r*L];
            }
	  if(r%2==1)
            {
	      T[(r-1)/2].re=y[r*L-1+j*r*L];
	      T[(r-1)/2].im=0.0;
            }
	  for(l=1;l<=r/2;l++)
            {
	      T[r-l].re= T[l-1].re;
	      T[r-l].im=-T[l-1].im;
            }
	  W.re=1.0;
	  W.im=0.0;
	  for(q=0;q<r;q++)
            {
	      a.re=R[q].re;
	      a.im=R[q].im;
	      z.re=T[0].re;
	      z.im=T[0].im;
	      for(l=1;l<r;l++)
                {
		  z.re=z.re+T[l].re*a.re+T[l].im*a.im;
		  z.im=z.im+T[l].im*a.re-T[l].re*a.im;
		  b.re=a.re*R[q].re-a.im*R[q].im;
		  b.im=a.im*R[q].re+a.re*R[q].im;
		  a.re=b.re;
		  a.im=b.im;
                }
	      x[L-1+q*N*L+j*L]=z.re*W.re+z.im*W.im;
	      z.re=W.re*w.re-W.im*w.im;
	      z.im=W.im*w.re+W.re*w.im;
	      W.re=z.re;
	      W.im=z.im;
            }
        }
    }
}