原创 时域抽取基2FFT算法C程序注解

                    a=-SIN_TAB[n-FFT_N/2];
  return a;
}
/******************************************************************
函数原型：void cos_tab(float pi)
函数功能：采用查表的方法计算一个数的余弦值
输入参数：pi 所要计算余弦值弧度值，范围0--2*PI，不满足时需要转换
输出参数：输入值pi的余弦值
******************************************************************/
float cos_tab(float pi)
{
   float a,pi2;
   pi2=pi+PI/2;
   if(pi2>2*PI)
     pi2-=2*PI;
   a=sin_tab(pi2);             //直接利用前述正弦表，简化处理
   return a;
}
/*****************************************************************
函数原型：void FFT(struct compx *xin)
函数功能：对输入的复数组进行快速傅里叶变换（FFT）
输入参数：*xin复数结构体组的首地址指针，struct型
输出参数：无
*****************************************************************/
void FFT(struct compx *xin)
{
  int f,m,nv2,nm1,i,k,l,j=0;
  struct compx u,w,t;
  
   nv2=FFT_N/2;             //变址运算，即把自然顺序变成倒位序，采用雷德算法
   nm1=FFT_N-1; 
   for(i=0;i

          {
           if(i)>;i++)>

                   {
                   t=xin[j];          
                   xin[j]=xin;
                   xin=t;
                  }
           k=nv2;                    //求j 的下一个倒位序
           while(k<=j)               //如果k<=j,表示j的最高位为1  
                 {          
                 j=j-k;                 //把最高位变成0
                 k=k/2;                 //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
                }
            j=j+k;                   //如果k>j，把0改为1
           }
                        
  {
    int le,lei,ip;                            //FFT运算核，使用蝶形运算完成FFT运算
    f=FFT_N;
    for(l=1;(f=f/2)!=1;l++)              //计算l的值，即计算蝶形级数
        for(m=1;m<=l;m++)                   // 控制蝶形结级数，没有第0级说法
            {                               //m表示第m级蝶形，l为蝶形级总数l=log（2）N
            le=2<<(m-1);         //le蝶形结距离，2^m，即第m级蝶形的蝶形结相距le点
            lei=le/2;                         //2^m-1，同一蝶形结中参加运算的两点的距离
            u.real=1.0;                        //u为蝶形结运算系数，初始值为1
            u.imag=0.0;
            //w.real=cos(PI/lei);                  //不适用查表法计算sin值和cos值
            // w.imag=-sin(PI/lei);
            w.real=cos_tab(PI/lei);   //w(值取决lei)为系数商，即当前系数u与前一个系数的商
            w.imag=-sin_tab(PI/lei);
            for(j=0;j<=lei-1;j++)              //控制计算不同种蝶形结，即计算系数不同的蝶形结
                    {
                    for(i=j;i<=FFT_N-1;i=i+le)      //控制同一蝶形结运算，即计算系数相同蝶形结
                            {
                            ip=i+lei;                          //i，ip分别表示参加蝶形运算的两个节点
                            t=EE(xin[ip],u);              //蝶形运算，详见公式，一次复数乘法
                            xin[ip].real=xin.real-t.real;    //两次复数加法，其中一次为减法
                            xin[ip].imag=xin.imag-t.imag;  //若用函数分装起来更简洁
                            xin.real=xin.real+t.real;
                            xin.imag=xin.imag+t.imag;
                             }
                    u=EE(u,w);                          //改变系数，进行下一个蝶形运算
                   }
            }
  }
 
}

/************************************************************
函数原型：void main()
函数功能：测试FFT变换，演示函数使用方法
输入参数：无
输出参数：无
************************************************************/
void main()  
{ 
  int i;
  create_sin_tab(SIN_TAB);
  for(i=0;i;i++)>)>