快速傅里叶变换法(FFT)

    快速傅里叶变换法（FFT）是离散傅立叶变换的一种快速计算方法，它能使N点DFT的乘法计算量由N2 次降为次。下图是采样点数为8点FFT时间抽取算法信号流图，本程序也是以这种形式设计的。

程序设计的基本思路是在程序开始时刻要求输入采样点数，如果采样点数是2的整数次方（不包括0次方），则要求输入采样点的数值，并根据采样点数分配响应的数组大小，计算迭代次数。然后对采样点进行逆二进制排序，再按上图所示的算法进行计算，程序流程图如下图所示：

本程序运用VC语言对程序进行设计，下面分别对程序设计中复数类的应用，判断和求迭代次数，逆二进制排序，蝶形运算进行具体说明。

1．复数类的应用

   C语言本身并不包含复数数据类型，但C语言可以根据需要定义自己的数据类型，本程序定义了一个复数结构体complex，包括实部real和虚部img两部分，代码如下：

typedef struct 

     {

       double real;

       double img;

     }complex;

在FFT程序设计中，复数类主要被用来计算两复数的加法和乘法以及旋转因子，其中，式中N=2的m+1次方，m代表计算流图的第m级，而k代表第k次蝶形运算，由于C中的math.h函数库中没有带参数的复数运算函数，所以本程序编写了带参数的复数运算cw(double x,double y)，用于计算，设计的基本思路，首先把e的次幂用欧拉公式化成三角函数，然后化复数乘法和除法运算为几个复数基本单元的加法运算和除法运算，其中运算的次数由函数输入参量double x决定。函数mul(complex x1, complex x2)用于计算复数的乘运算。

2．判断和求迭代次数

本程序编写iternumb(int numb)函数对采样点数进行判断，如果采样点数不符合2的整数次方或采样点数为1或0，则跳出程序，程序设计基本思路是对输入采样点数的十进制形式进行模2运算和除法运算，在除法运算结果大于1之前，一旦模2运算的结果等于1，则说明输入采样点数不符合要求，而如果符合要求，则把出发结果存入数组当中，函数代码如下：

int iternumb(int numb)

    {

     int iternumb1=0;

     if((numb==0)||(numb==1))

       {

          printf("numb error!\\n");

           exit(0);

        }

      while ((numb!=0)&&(numb!=1))

     {

        if (numb%2)

        {

           printf("numb error!\\n");

           exit(0);

        }

        numb=numb/2;

        iternumb1=iternumb1+1;

     }

     return iternumb1;

  }

3．码位倒置

  在逆二进制排序程序中，设置for循环分别将输入数据数组input[i]的索引号i进行模2运算，所得的结果按逆序存入inverse[ ]数组（存入inverse[ ]数组的顺序是从数组尾部开始）。

然后再将inverse[ ]中的二进制数转换为十进制数，并以此数为A[j]的索引号，令A[j]= input[i]，从而实现了逆二进制排序。

4．蝶形运算

蝶形运算是FFT中最基本的一个运算单元。在FFT程序设计中要找到蝶形运算地址与第几级迭代，第几组之间的关系。在第m级，有如下关系：

式中，m为第级数，可以用函数cw(double x,double y)计算出来，p、q分别为单元地址，p、q之间的距离为2m，即。根据上述关系可运用3次for循环实现上述计算，最外层为级数m的循环，m的取值为0～log2N,中间为组数i循环，每一级i的取值为1～，最内层为旋转因子中的r循环，每一组中r的取值为0～2m -1，循环的部分代码如下：

itm=iternumb(N);

       sign1=itm;

for(m=0;m      {

        X=pow(2,m);

for(i=1;i        {

           p=2*X*(i-1);

           q=X+p;

for(r=0;r             {

                   w=cw(r,X);

                   tempA3=mul(w,A[q]);

                   tempA1.real=A[p].real+tempA3.real;

                   tempA1.img=A[p].img+tempA3.img;

                   tempA2.real=A[p].real-tempA3.real;

                   tempA2.img=A[p].img-tempA3.img;

                   A[p].real=tempA1.real;

                   A[p].img=tempA1.img;

                   A[q].real=tempA2.real;

                   A[q].img=tempA2.img;

                   p=p+1;

                   q=X+p;

              }

          }

      }

5．执行并检验程序

数入采样点数为6，程序运行如下：

输入采样点数为8，分别为1，2，3，4，5，6，7，8，程序运行如下：

在Matlab中运算的结果如下：

结果与所设计程序一致。

6. 结论

通过对程序的设计以及对运行结果的验证可以说明所设计的程序运行结果正确，完成了对采样点进行FFT变换的要求。

附录

FFT程序代码：

#include

#include

#include

typedef struct 

   {

    double real;

    double img;

   }complex;

int iternumb();

complex cw();

complex mul();

#define pi 4*atan(1)

int iternumb(int numb)

  {

     int iternumb1=0;

     if((numb==0)||(numb==1))

       {

          printf("numb error!\\n");

           exit(0);

        }

      while ((numb!=0)&&(numb!=1))

     {

        if (numb%2)

        {

           printf("numb error!\\n");

           exit(0);

        }

        numb=numb/2;

        iternumb1=iternumb1+1;

     }

     return iternumb1;

  }

complex cw(double x,double y)

{

  complex z,z1,z2;

   int i;

   double t,t1;

   t=2*pi/(2*y);

   t1=x;

   z.real=cos(t);

   z.img=-sin(t);

      z1.real=cos(t);

    z1.img=-sin(t);

         if (t1==0)

         {

             z.real=1;

             z.img=0;

         }

         if(t1==1)

        {

            z.real=cos(t);

            z.img=-sin(t);

        }

        if ((t1!=0)&&(t1!=1))

        {

         for(i=2;i            { 

                z2.real=z.real*z1.real-z.img*z1.img;

                z2.img=z.real*z1.img+z.img*z1.real;

                z.real=z2.real;

                 z.img=z2.img;

            }

      }

    return z;

 }

complex mul(complex x1, complex x2)

 {

      complex  z;

      z.real=x1.real*x2.real-x1.img*x2.img;

      z.img=x1.real*x2.img+x1.img*x2.real;

      return z;

  }

int main()

{

  int a,b,c,N,X,m,i,k,p,q,itm,temp2,temp3,temp4,temp1,sign,*inverse,*input;

  double temp,r,sign1;

  complex tempA1,tempA2,tempA3,tempA4,w;

  complex *A;

  printf("please input number N:\\n");

  scanf("%d",&N);

  input=(int*)malloc(sizeof(int)* N);

  A=(complex*)malloc(sizeof(complex)*N);

  itm=iternumb(N);

  inverse =(int*)malloc(sizeof(int)* itm);

  printf("please input digital:\\n");

for(a=0;a     {

     scanf("%d",input+a);

     }

for(a=0;a  {

       temp2=a;

for(b=0;b          {

            temp4=temp2%2;

            inverse[itm-1-b]=temp4;

            temp2=temp2/2;

          }

        temp1=0;

        sign=1;

for(c=0;c           {

              temp1=sign*inverse[c]+ temp1;

                sign=sign*2;

            }

           A[temp1].real=input[a];

           A[temp1].img=0;

}

       itm=iternumb(N);

       sign1=itm;

for(m=0;m      {

        X=pow(2,m);

for(i=1;i        {

           p=2*X*(i-1);

           q=X+p;

for(r=0;r             {

                   w=cw(r,X);

                   tempA3=mul(w,A[q]);

                   tempA1.real=A[p].real+tempA3.real;

                   tempA1.img=A[p].img+tempA3.img;

                   tempA2.real=A[p].real-tempA3.real;

                   tempA2.img=A[p].img-tempA3.img;

                   A[p].real=tempA1.real;

                   A[p].img=tempA1.img;

                   A[q].real=tempA2.real;

                   A[q].img=tempA2.img;

                   p=p+1;

                   q=X+p;

              }

          }

      }

for(k=0;k       {

          printf("A[%d]=%f+(%f)i\\n",k,A[k].real, A[k].img);

        }

      return 0;