/* VC 下运行通过 2007。4。26 把radix4中某些层for->while等效 14:51*/
/*****************fft programe*********************/
void rad4(double x[],double y[],double TAB[],int i,int n2,int adr1,int adr2,int adr3)
{
 int i1,i2,i3;
 double co1,co2,co3,si1,si2,si3;
 double r1,r2,r3,r4,s1,s2,s3,s4;

    i1=i+n2;i2=i1+n2;i3=i2+n2;                       /* 数据位置 */

 co1=TAB[2*adr1];co2=TAB[2*adr2];co3=TAB[2*adr3]; /* 旋转因子 */
 si1=TAB[2*adr1+1];si2=TAB[2*adr2+1];si3=TAB[2*adr3+1];

 r1=x[i]+x[i2];r3=x[i]-x[i2];
 s1=y[i]+y[i2];s3=y[i]-y[i2];
 r2=x[i1]+x[i3];r4=x[i1]-x[i3];
 s2=y[i1]+y[i3];s4=y[i1]-y[i3];

 x[i]=r1+r2; y[i]=s1+s2; 
 r2=r1-r2; r1=r3-s4;r3=r3+s4;
 s2=s1-s2; s1=s3+r4;s3=s3-r4;
 x[i1]=co1*r3+si1*s3;y[i1]=co1*s3-si1*r3;
 x[i2]=co2*r2+si2*s2;y[i2]=co2*s2-si2*r2;
 x[i3]=co3*r1+si3*s1;y[i3]=co3*s1-si3*r1;
}

void switchadr(double x[],double y[],int n) // 变换后,如果想顺序输出就调用这个
{
 double tempx,tempy;
 int i,j,k;
 for(i=0,j=0;i<=n-1;i++)
 {
  if(i<j)
  {
   tempx=x[j];  tempy=y[j]; x[j]=x[i];y[j]=y[i]; x[i]=tempx;  y[i]=tempy;
  }
  k=n/4;
  
  while(3*k<(j+1)&&(j>0) )
  {
   j=j-3*k; k=k/4;
  }
  j=j+k;
 }
}

void radix4(double x[],double y[],double TAB[],int n) // 核心变换代码
{
 int k,e,n1,n2;
 int a,b,c;
    n2=n;
 e=1;    /* 1 -- n 有效 */
 k=n;
 while(k>1) /*内部循环没有使用K值，只是用来控制循环次数而已,log4\n次*/
 {
  int j;
  n1=n2; n2=n2/4; 
  for(j=1;j<=n2;j++)
  {
   int i;
   a=(j-1)*e;b=a+a; c=a+b;
   a &= (n-1);  b &= (n-1); c &= (n-1); /* 旋转因子周期 N */

            i = (j-1);
   do {
    rad4(x,y,TAB,i,n2,a,b,c);
    i+=n1;
   }while(i<n);
  }
  e <<= 2; k >>= 2;
 }
}

评论