cfft.cpp

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

      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-5 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_radix5(complex *x,complex *y,int L,int N)
{
  complex w1,w2,w3,w4;
  complex t0,t1,t2,t3,t4;
  complex s0,s1,s2,s3,s4;
  complex *x0,*x1,*x2,*x3,*x4;
  complex *y0,*y1,*y2,*y3,*y4;
  const real a1=0.3090169943749474241022934171830L;  // cos(-2*pi/5)
  const real b1=-0.9510565162951535721164393333793L; // sin(-2*pi/5)
  const real a2=-0.8090169943749474241022934171825L; // cos(-4*pi/5)
  const real b2=-0.5877852522924731291687059546394L; // sin(-2*pi/5)
    
  for(int j=0;j<N;j++)
    { 
      x0=x+j*L;
      x1=x0+L*N;
      x2=x1+L*N;
      x3=x2+L*N;
      x4=x3+L*N;
      y0=y+5*L*j;
      y1=y0+L;
      y2=y1+L;
      y3=y2+L;
      y4=y3+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;
      t4.re=x4[0].re;
      t4.im=x4[0].im;
      s0.re=t0.re;
      s0.im=t0.im;
      s1.re=t1.re+t4.re;
      s1.im=t1.im+t4.im;
      s2.re=t2.re+t3.re;
      s2.im=t2.im+t3.im;
      s3.re=t1.re-t4.re;
      s3.im=t1.im-t4.im;
      s4.re=t2.re-t3.re;
      s4.im=t2.im-t3.im;
      t0.re=s0.re+s1.re+s2.re;
      t0.im=s0.im+s1.im+s2.im;
      t1.re=s0.re+a1*s1.re+a2*s2.re;
      t1.im=s0.im+a1*s1.im+a2*s2.im;
      t2.re=s0.re+a2*s1.re+a1*s2.re;
      t2.im=s0.im+a2*s1.im+a1*s2.im;
      t3.re=b1*s3.re+b2*s4.re;
      t3.im=b1*s3.im+b2*s4.im;
      t4.re=b2*s3.re-b1*s4.re;
      t4.im=b2*s3.im-b1*s4.im;
      y0[0].re=t0.re;
      y0[0].im=t0.im;
      y1[0].re=t1.re-t3.im;
      y1[0].im=t1.im+t3.re;
      y2[0].re=t2.re-t4.im;
      y2[0].im=t2.im+t4.re;
      y3[0].re=t2.re+t4.im;
      y3[0].im=t2.im-t4.re;
      y4[0].re=t1.re+t3.im;
      y4[0].im=t1.im-t3.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; 
	  w4.re = w3.re*w1.re-w3.im*w1.im;
	  w4.im = w3.re*w1.im+w3.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;
	  t4.re = w4.re*x4[i].re-w4.im*x4[i].im;
	  t4.im = w4.re*x4[i].im+w4.im*x4[i].re;
	  s0.re=t0.re;
	  s0.im=t0.im;
	  s1.re=t1.re+t4.re;
	  s1.im=t1.im+t4.im;
	  s2.re=t2.re+t3.re;
	  s2.im=t2.im+t3.im;
	  s3.re=t1.re-t4.re;
	  s3.im=t1.im-t4.im;
	  s4.re=t2.re-t3.re;
	  s4.im=t2.im-t3.im;
	  t0.re=s0.re+s1.re+s2.re;
	  t0.im=s0.im+s1.im+s2.im;
	  t1.re=s0.re+a1*s1.re+a2*s2.re;
	  t1.im=s0.im+a1*s1.im+a2*s2.im;
	  t2.re=s0.re+a2*s1.re+a1*s2.re;
	  t2.im=s0.im+a2*s1.im+a1*s2.im;
	  t3.re=b1*s3.re+b2*s4.re;
	  t3.im=b1*s3.im+b2*s4.im;
	  t4.re=b2*s3.re-b1*s4.re;
	  t4.im=b2*s3.im-b1*s4.im;
	  y0[i].re=t0.re;
	  y0[i].im=t0.im;
	  y1[i].re=t1.re-t3.im;
	  y1[i].im=t1.im+t3.re;
	  y2[i].re=t2.re-t4.im;
	  y2[i].im=t2.im+t4.re;
	  y3[i].re=t2.re+t4.im;
	  y3[i].im=t2.im-t4.re;
	  y4[i].re=t1.re+t3.im;
	  y4[i].im=t1.im-t3.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-5 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_radix5(complex *x,complex *y,int L,int N)
{
  complex w1,w2,w3,w4;
  complex t0,t1,t2,t3,t4;
  complex s0,s1,s2,s3,s4;
  complex *x0,*x1,*x2,*x3,*x4;
  complex *y0,*y1,*y2,*y3,*y4;
  const real a1=0.3090169943749474241022934171830L;
  const real b1=-0.9510565162951535721164393333793L;
  const real a2=-0.8090169943749474241022934171825L;
  const real b2=-0.5877852522924731291687059546394L;
    
  for(int j=0;j<N;j++)
    { 
      x0=x+j*L;
      x1=x0+L*N;
      x2=x1+L*N;
      x3=x2+L*N;
      x4=x3+L*N;
      y0=y+5*L*j;
      y1=y0+L;
      y2=y1+L;
      y3=y2+L;
      y4=y3+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;
      t4.re=x4[0].re;
      t4.im=x4[0].im;
      s0.re=t0.re;
      s0.im=t0.im;
      s1.re=t1.re+t4.re;
      s1.im=t1.im+t4.im;
      s2.re=t2.re+t3.re;
      s2.im=t2.im+t3.im;
      s3.re=t1.re-t4.re;
      s3.im=t1.im-t4.im;
      s4.re=t2.re-t3.re;
      s4.im=t2.im-t3.im;
      t0.re=s0.re+s1.re+s2.re;
      t0.im=s0.im+s1.im+s2.im;
      t1.re=s0.re+a1*s1.re+a2*s2.re;
      t1.im=s0.im+a1*s1.im+a2*s2.im;
      t2.re=s0.re+a2*s1.re+a1*s2.re;
      t2.im=s0.im+a2*s1.im+a1*s2.im;
      t3.re=b1*s3.re+b2*s4.re;
      t3.im=b1*s3.im+b2*s4.im;
      t4.re=b2*s3.re-b1*s4.re;
      t4.im=b2*s3.im-b1*s4.im;
      y0[0].re=t0.re;
      y0[0].im=t0.im;
      y1[0].re=t1.re+t3.im;
      y1[0].im=t1.im-t3.re;
      y2[0].re=t2.re+t4.im;
      y2[0].im=t2.im-t4.re;
      y3[0].re=t2.re-t4.im;
      y3[0].im=t2.im+t4.re;
      y4[0].re=t1.re-t3.im;
      y4[0].im=t1.im+t3.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; 
	  w4.re = w3.re*w1.re-w3.im*w1.im;
	  w4.im = w3.re*w1.im+w3.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;
	  t4.re = w4.re*x4[i].re-w4.im*x4[i].im;
	  t4.im = w4.re*x4[i].im+w4.im*x4[i].re;
	  s0.re=t0.re;
	  s0.im=t0.im;
	  s1.re=t1.re+t4.re;
	  s1.im=t1.im+t4.im;
	  s2.re=t2.re+t3.re;
	  s2.im=t2.im+t3.im;
	  s3.re=t1.re-t4.re;
	  s3.im=t1.im-t4.im;
	  s4.re=t2.re-t3.re;
	  s4.im=t2.im-t3.im;
	  t0.re=s0.re+s1.re+s2.re;
	  t0.im=s0.im+s1.im+s2.im;
	  t1.re=s0.re+a1*s1.re+a2*s2.re;
	  t1.im=s0.im+a1*s1.im+a2*s2.im;
	  t2.re=s0.re+a2*s1.re+a1*s2.re;
	  t2.im=s0.im+a2*s1.im+a1*s2.im;
	  t3.re=b1*s3.re+b2*s4.re;
	  t3.im=b1*s3.im+b2*s4.im;
	  t4.re=b2*s3.re-b1*s4.re;
	  t4.im=b2*s3.im-b1*s4.im;
	  y0[i].re=t0.re;
	  y0[i].im=t0.im;
	  y1[i].re=t1.re+t3.im;
	  y1[i].im=t1.im-t3.re;
	  y2[i].re=t2.re+t4.im;
	  y2[i].im=t2.im-t4.re;
	  y3[i].re=t2.re-t4.im;
	  y3[i].im=t2.im+t4.re;
	  y4[i].re=t1.re-t3.im;
	  y4[i].im=t1.im+t3.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-6 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_radix6(complex *x,complex *y,int L,int N)
{
  complex w1,w2,w3,w4,w5;
  complex t0,t1,t2,t3,t4,t5;
  complex s0,s1,s2,s3,s4,s5;
  complex *x0,*x1,*x2,*x3,*x4,*x5;
  complex *y0,*y1,*y2,*y3,*y4,*y5;
  const real a= 0.5L;
  const real b=-0.5L*root3;
    
  for(int j=0;j<N;j++)
    { 
      x0=x+j*L;
      x1=x0+L*N;
      x2=x1+L*N;
      x3=x2+L*N;
      x4=x3+L*N;
      x5=x4+L*N;
      y0=y+6*L*j;
      y1=y0+L;
      y2=y1+L;
      y3=y2+L;
      y4=y3+L;
      y5=y4+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;
      t4.re=x4[0].re;
      t4.im=x4[0].im;
      t5.re=x5[0].re;
      t5.im=x5[0].im;
      s0.re=t0.re+t3.re;
      s0.im=t0.im+t3.im;
      s1.re=t0.re-t3.re;
      s1.im=t0.im-t3.im;
      s2.re=t1.re+t2.re+t4.re+t5.re;
      s2.im=t1.im+t2.im+t4.im+t5.im;
      s3.re=t1.re-t2.re-t4.re+t5.re;
      s3.im=t1.im-t2.im-t4.im+t5.im;
      s4.re=t1.re+t2.re-t4.re-t5.re;
      s4.im=t1.im+t2.im-t4.im-t5.im;
      s5.re=t1.re-t2.re+t4.re-t5.re;
      s5.im=t1.im-t2.im+t4.im-t5.im;
      t0.re=s0.re+s2.re;
      t0.im=s0.im+s2.im;
      t1.re=s1.re+a*s3.re;
      t1.im=s1.im+a*s3.im;
      t2.re=s0.re-a*s2.re;
      t2.im=s0.im-a*s2.im;
      t3.re=s1.re-s3.re;
      t3.im=s1.im-s3.im;
      t4.re=b*s4.re;
      t4.im=b*s4.im;
      t5.re=b*s5.re;
      t5.im=b*s5.im;
      y0[0].re=t0.re;
      y0[0].im=t0.im;
      y1[0].re=t1.re-t4.im;
      y1[0].im=t1.im+t4.re;
      y2[0].re=t2.re-t5.im;
      y2[0].im=t2.im+t5.re;
      y3[0].re=t3.re;
      y3[0].im=t3.im;
      y4[0].re=t2.re+t5.im;
      y4[0].im=t2.im-t5.re;
      y5[0].re=t1.re+t4.im;
      y5[0].im=t1.im-t4.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; 
	  w4.re = w2.re*w2.re-w2.im*w2.im;
	  w4.im = 2.0*w2.re*w2.im;
	  w5.re = w4.re*w1.re-w4.im*w1.im;
	  w5.im = w4.re*w1.im+w4.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;
	  t4.re = w4.re*x4[i].re-w4.im*x4[i].im;
	  t4.im = w4.re*x4[i].im+w4.im*x4[i].re;
	  t5.re = w5.re*x5[i].re-w5.im*x5[i].im;
	  t5.im = w5.re*x5[i].im+w5.im*x5[i].re;
	  s0.re=t0.re+t3.re;
	  s0.im=t0.im+t3.im;
	  s1.re=t0.re-t3.re;
	  s1.im=t0.im-t3.im;
	  s2.re=t1.re+t2.re+t4.re+t5.re;
	  s2.im=t1.im+t2.im+t4.im+t5.im;
	  s3.re=t1.re-t2.re-t4.re+t5.re;
	  s3.im=t1.im-t2.im-t4.im+t5.im;
	  s4.re=t1.re+t2.re-t4.re-t5.re;
	  s4.im=t1.im+t2.im-t4.im-t5.im;
	  s5.re=t1.re-t2.re+t4.re-t5.re;
	  s5.im=t1.im-t2.im+t4.im-t5.im;
	  t0.re=s0.re+s2.re;
	  t0.im=s0.im+s2.im;
	  t1.re=s1.re+a*s3.re;
	  t1.im=s1.im+a*s3.im;
	  t2.re=s0.re-a*s2.re;
	  t2.im=s0.im-a*s2.im;
	  t3.re=s1.re-s3.re;
	  t3.im=s1.im-s3.im;
	  t4.re=b*s4.re;
	  t4.im=b*s4.im;
	  t5.re=b*s5.re;
	  t5.im=b*s5.im;
	  y0[i].re=t0.re;
	  y0[i].im=t0.im;
	  y1[i].re=t1.re-t4.im;
	  y1[i].im=t1.im+t4.re;
	  y2[i].re=t2.re-t5.im;
	  y2[i].im=t2.im+t5.re;
	  y3[i].re=t3.re;
	  y3[i].im=t3.im;
	  y4[i].re=t2.re+t5.im;
	  y4[i].im=t2.im-t5.re;
	  y5[i].re=t1.re+t4.im;
	  y5[i].im=t1.im-t4.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-6 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_radix6(complex *x,complex *y,int L,int N)
{
  complex w1,w2,w3,w4,w5;
  complex t0,t1,t2,t3,t4,t5;
  complex s0,s1,s2,s3,s4,s5;
  complex *x0,*x1,*x2,*x3,*x4,*x5;
  complex *y0,*y1,*y2,*y3,*y4,*y5;
  const real a= 0.5L;
  const real b=-0.5L*root3;
    
  for(int j=0;j<N;j++)
    { 
      x0=x+j*L;
      x1=x0+L*N;
      x2=x1+L*N;
      x3=x2+L*N;
      x4=x3+L*N;
      x5=x4+L*N;
      y0=y+6*L*j;
      y1=y0+L;
      y2=y1+L;
      y3=y2+L;
      y4=y3+L;
      y5=y4+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;
      t4.re=x4[0].re;
      t4.im=x4[0].im;
      t5.re=x5[0].re;
      t5.im=x5[0].im;
      s0.re=t0.re+t3.re;
      s0.im=t0.im+t3.im;
      s1.re=t0.re-t3.re;
      s1.im=t0.im-t3.im;
      s2.re=t1.re+t2.re+t4.re+t5.re;
      s2.im=t1.im+t2.im+t4.im+t5.im;
      s3.re=t1.re-t2.re-t4.re+t5.re;
      s3.im=t1.im-t2.im-t4.im+t5.im;
      s4.re=t1.re+t2.re-t4.re-t5.re;
      s4.im=t1.im+t2.im-t4.im-t5.im;
      s5.re=t1.re-t2.re+t4.re-t5.re;
      s5.im=t1.im-t2.im+t4.im-t5.im;
      t0.re=s0.re+s2.re;
      t0.im=s0.im+s2.im;
      t1.re=s1.re+a*s3.re;
      t1.im=s1.im+a*s3.im;
      t2.re=s0.re-a*s2.re;
      t2.im=s0.im-a*s2.im;
      t3.re=s1.re-s3.re;
      t3.im=s1.im-s3.im;
      t4.re=b*s4.re;
      t4.im=b*s4.im;
      t5.re=b*s5.re;
      t5.im=b*s5.im;
      y0[0].re=t0.re;
      y0[0].im=t0.im;
      y1[0].re=t1.re+t4.im;
      y1[0].im=t1.im-t4.re;
      y2[0].re=t2.re+t5.im;
      y2[0].im=t2.im-t5.re;
      y3[0].re=t3.re;
      y3[0].im=t3.im;
      y4[0].re=t2.re-t5.im;
      y4[0].im=t2.im+t5.re;
      y5[0].re=t1.re-t4.im;
      y5[0].im=t1.im+t4.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; 
	  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; 
	  t0.re=x0[i].re;
	  t0.im=x0[i].im;