fft变换c源代码22.txt

快速傅利叶变换即FFT变换C源代码,可能实现数据的快速傅利叶变换

发信人: horluk (不才), 信区: Signal       
标  题: Re: 哪里有fft2的C＋＋源代码？谢谢
发信站: BBS 水木清华站 (Mon Jul  1 23:07:58 2002)

前几天写的，基2实数序列FFT程序，测试过了
double * Xin : 需要变换的实数序列
double * SineTab : SineTab[i]=exp(2*PAI*i/N); i : 0 ～ N/2-1
unsigned int N : 序列长度,须为 2 的整指数幂
unsigned int M : M=log2(N)
输出：在 Xin[] 中，排列次序－》  R[0],R[1],...R[N/2],I[N/2-1],...,I[1]

void FFT_R(double * Xin ,double * SineTab, unsigned int N, unsigned int M){

 unsigned int m , i , j ; /* 当前级数，元素位置或者块首元素位置，元素在块中的
位置 */
 unsigned int uiCoeGap_BlkNum=1;  /*相邻两蝶行系数指数的差，也等于本级中的块
数 */
 unsigned int BN=1;  /*，每块中蝶行的数目，也等于每个蝶行两元素之间的位置差 
*/
 unsigned int k; /*内循环实部所在蝶行系数指数*/
 unsigned int N2=N>>1; /*N/2*/
 unsigned int N4=N>>2; /*N/4*/
 double dTemp, dTemp1,dTemp2;/* 临时变量 */
 //位码倒序排列
 for(i=0 ; i<N ; i++){
  j = BitInverse(i , M); /*获得i的低M位位码倒序排列*/
  if(i<j){dTemp= Xin[j]; Xin[j] = Xin[i]; Xin[i] = dTemp; }
  }
 //第一级
 for(i=0 ; i<N ; i+=2){
  dTemp=Xin[i]+Xin[i+1]; Xin[i+1]=Xin[i]-Xin[i+1]; Xin[i]=dTemp;
 }
 //第二级
 for(i=0 ; i<N ; i+=4){
  dTemp=Xin[i]+Xin[i+2];  Xin[i+2]=Xin[i]-Xin[i+2];  Xin[i]=dTemp;
  Xin[i+3] = -Xin[i+3];
 }
 //主循环
 m=3;//从第三级开始
 uiCoeGap_BlkNum <<= (M-m);     /* pow(2, (M-m)) */
 BN = 4; /* pow(2, (m-1)) */
 for(m=3 ; m<=M ; m++){//重复 M-2 次
     for(i=0 ; i<N ; i += BN<<1 ){ /*Block*/
      /* Xin[i] is the first element of this block */
      //第j=0个蝶形，Xin[i+j] 和 Xin[i+j+BN] 全是实数，系数为Wn0＝1；
      dTemp=Xin[i]+Xin[i+BN];
      Xin[i+BN]=Xin[i]-Xin[i+BN];
      Xin[i]=dTemp;
      //第 j=1 到 j= (BN>>1)-1 及其虚部所对应的第j=(BN-j)个双蝶行
      k = uiCoeGap_BlkNum; //系数初始化
      for(j=1 ; j<BN/2 ; j++ ){ /*Butterflys*/
         /* 第 j 和 (-j+BN)双蝶行*/
         /* Xin[i+j],    Xin[i+j+BN],   Wn[k]＝SineTab[k+N4]-jsin[k] */
         /* Xin[i-j+BN], Xin[i-j+2*BN], Wn[kI]-SineTab[k+N4]-jsin[k] */
         dTemp1 = SineTab[k+N4]*Xin[i+j+BN] +    SineTab[k]*Xin[i-j+2*BN];
         dTemp2 =  - SineTab[k]*Xin[i+j+BN] + SineTab[k+N4]*Xin[i-j+2*BN];
         dTemp = Xin[i-j+BN] ;
         Xin[i+j+BN]   = dTemp2 - dTemp;  /*Xin[i+j+BN]=I[i-j+BN]*/
         Xin[i-j+2*BN] = dTemp2 + dTemp;  /*Xin[i-j+2*BN]=I[i+j]*/
         dTemp = Xin[i+j];
         Xin[i+j]    =  dTemp + dTemp1;  /*Xin[i+j]＝R[i+j]*/
         Xin[i-j+BN] =  dTemp - dTemp1;  /*Xin[i-j+BN]=R[i+j+BN]*/
         k +=uiCoeGap_BlkNum;
      }
      //第j=(BN>>1)个蝶行，k=(uiCoeGap_BlkNum)(BN>>1)=pow(2,M-2)
      //Xin[i+j] 和 Xin[i+j+BN] 都是实数
      Xin[i+j+BN] = - Xin[i+j+BN];//Wn[k]=-j, k=N/4
      //一块结束
   }//一级结束
   BN <<=1; //每块中蝶行数增倍，蝶行两元素位置差增倍
   uiCoeGap_BlkNum >>=1; //块数减半，系数间距减半
 }
}
unsigned int BitInverse(unsigned int unIput, unsigned int unBitNum){
 unsigned int unOutPut=0;
 int i=0;
 int j=0;
 for(i=unBitNum-1 , j=0 ; i>=0 ; i--, j++){
  if( unIput & (1<<i) ) unOutPut += 1<<j;
 }
 return unOutPut;
}