fanfft.txt

将共轭对称性复序列进行快速傅里叶反变换,得出变换结果为实序列。

//将共轭对称性复序列进行快速傅里叶反变换,得出变换结果为实序列。
void fan_fft(double *x,int n,int m)
//x--双精度实型一维数组，长度为n，开始存放具有共轭对称性的复序列X（K）的前n/2+1个值，存储顺序为*/
//re(0),re(1)...re(N/2),im(N/2-1),....im(1))，其中re(0)=X(0),re(n/2)=X(n/2),最后存放变换结果实序列x(i)(i=0,n-1)*/
//n--整型变量，数据长度，必须是2的整数次幂,
//m-整型变量,n的2整数次幂*/
{
  int i,j,k,i1,i2,i3,i4,i5,i6,i7,i8,n2,n4,n8,id,is;
  double a,e,a3,t1,t2,t3,t4,t5,cc1,cc3,ss1,ss3;
  n2=2*n;
  for(k=1;k<m;k++)
  { is=0;
    id=n2;
    n2=n2/2;
    n4=n2/4;
    n8=n4/2;
   e=6.28318530718/n2;
  do{
   for(i=is;i<n;i+=id) //i控制进入下一组，对应的旋转因子相同
		{i1=i;
                 i2=i1+n4;
		  i3=i2+n4;
		  i4=i3+n4;
		  t1=x[i1]-x[i3];
		  x[i1]=x[i1]+x[i3];
                  x[i2]=2*x[i2];
		  x[i3]=t1-2*x[i4];
		  x[i4]=t1+2*x[i4];
                  if (n4==1) continue;
                   i1+=n8;
                   i2+=n8;
                   i3+=n8;
                   i4+=n8;
                   t1=(x[i2]-x[i1])/sqrt(2.0);
                   t2=(x[i4]+x[i3])/sqrt(2.0);
                  x[i1]=x[i1]+x[i2];
                  x[i2]=x[i4]-x[i3];
                  x[i3]=2*(-t2-t1);
                  x[i4]=2*(-t2+t1);
                  }
                  is=2*id-n2;
                  id=4*id;
                  } while (is<(n-1));
                  a=e;
                   for (j=1;j<n8;j++)
                   {a3=3*a;
                   cc1=cos(a);
                   ss1=sin(a);
                   cc3=cos(a3);
                   ss3=sin(a3);
                   a=(j+1)*e;
                   is=0;
                   id=2*n2;
                   do
                   { for (i=is;i<=(n-1);i=i+id)
                   {i1=i-j;
                   i2=i1+n4;
                   i3=i2+n4;
                   i4=i3+n4;
                   i5=i+n4-j;
                   i6=i5+n4;
                   i7=i6+n4;
                   i8=i7+n4;
                   t1=x[i1]-x[i6];
                   x[i1]=x[i1]+x[i6];
                   t2=x[i5]-x[i2];
                   x[i5]=x[i2]+x[i5];
                   t3=x[i8]+x[i3];
                   x[i6]=x[i8]-x[i3];
                   t4=x[i4]+x[i7];
                   x[i2]=x[i4]-x[i7];
                   t5=t1-t4;
                   t1=t1+t4;
                   t4=t2-t3;
                   t2=t2+t3;
                   x[i3]=t5*cc1+t4*ss1;
                   x[i7]=-t4*cc1+t5*ss1;
                   x[i4]=t1*cc3-t2*ss3;
                   x[i8]=t2*cc3+t1*ss3;
                   }
                   is=2*id-n2;
                   id=4*id;
                   }
                   while (is<(n-1));
                   }
                   }
                   is=0;
                   id=4;
                   do
                   {for (i=is;i<n;i=i+id)
                   {i1=i+1;
                   t1=x[i];
                   x[i]=t1+x[i1];
                   x[i1]=t1-x[i1];
                   }
                   is=2*id-2;
                   id=4*id;
                   } while (is<(n-1));
                   for (j=0,i=0;i<(n-1);i++)
                   {if (i<j)
                   {t1=x[j];
                   x[j]=x[i];
                   x[j]=t1;
                   }
                   k=n/2;
                   while (k<(j+1))
                   {j=j-k;
                   k=k/2;
                   }
                   j=j+k;
                   }
                   for (i=0;i<n;i++)
                   x[i]=x[i]/n;
                   }