DFT（离散傅里叶变换）与FFT（快速傅里叶变换）运算效率对比

#include "stdio.h"
#include "math.h"
#include "time.h"
#define N 2048                     /* 测试时可修改N大小，如2048、4096等 */
#define PI 3.1415926

//#include "dft_sub.c"  /* 从dft_sub.c读入DFT子函数*/
//#include "fft_sub.c"  /* 从fft_sub.c读入FFT子函数*/
main()
{ void dft(float xr[],float xi[],int flag),fft(float sr[],float sx[],int m0,int inv);
  float xr[N],xi[N]={0};
  int i,n,k,Loop;
  double tt1,tt2,tt3,tt4;
  FILE *fp;
  k=log(N*1.0)/log(2.0);
  printf("%d points DFT and FFT, and k=%d\n",N,k);
  printf("Input Loop Number\n");   /* 测试时循环次数可以设大一点 */
  scanf("%d",&Loop);

  for(n=0;n<=N/2-1;n++)
 { xr[n]=0.0;}
  for(n=3*N/4;n<=N-1;n++)
 { xr[n]=0.0;}
  for(n=N/2;n<=3*N/4-1;n++)
 { xr[n]=1.0;}

  // DFT Calculating time.

  tt1=(double)clock();
  for(i=0;i<Loop;i++)
  { dft(xr,xi,1);
    dft(xr,xi,-1);
  }
  tt2=(double)clock();

  printf("DFT calculating time: %12.3f (ms)\n",1000*(tt2-tt1)/CLK_TCK);

// FFT Calculating time.
 
  tt3=(double)clock();
  printf("TT3: %12.3f (ms)\n",1000*tt3/CLK_TCK);

  for(i=0;i<Loop;i++)
  { fft(xr,xi,k,1);
    fft(xr,xi,k,-1);
  }
  tt4=(double)clock();

  printf("TT4: %12.3f (ms)\n",1000*tt4/CLK_TCK);

  printf("FFT calculating time: %12.3f (ms)\n",1000*(tt4-tt3)/CLK_TCK);

  printf("FFT/DFT Calculating Time Ratio: %12.3f (ms)\n",(tt2-tt1)/(tt4-tt3));


  fp=fopen("test.txt","w");
  for(n=0;n<N;n++)
   fprintf(fp,"%12d %12.4f\n",n,xr[n]);
  fclose(fp);
    return 0;
}

// fft_sub.c

void fft(float sr[],float sx[],int m0,int inv)
{
 int i,j,lm,li,k,lmx,np,lix,mm2;
 double t1,t2,c,s,cv,st,ct;
 if(m0<0)
  return;
 lmx=1;
 for(i=1;i<=m0;++i)
        lmx+=lmx;   //form 2**m0
    cv=2.0*PI/(double)lmx;
 ct=cos(cv); st=-inv*sin(cv);
 np=lmx;mm2=m0-2;
 /* fft butterfly numeration */
 for(i=1;i<=mm2;++i)
 {
  lix=lmx;lmx/=2;
  c=ct;s=st;
  for(li=0;li<np;li+=lix)
  {
   j=li;k=j+lmx;
   t1=sr[j]-sr[k];t2=sx[j]-sx[k];
   sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;
   ++j;++k;
   t1=sr[j]-sr[k];t2=sx[j]-sx[k];
   sr[j]+=sr[k];sx[j]+=sx[k];
   sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
  }
  for(lm=2;lm<lmx;++lm)
  {
   cv=c;c=ct*c-st*s;s=st*cv+ct*s;
   for(li=0;li<np;li+=lix)
   {
    j=li+lm;k=lmx+j;
    t1=sr[j]-sr[k];t2=sx[j]-sx[k];
    sr[j]+=sr[k];sx[j]+=sx[k];
    sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
   }
  }
  cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;
 }
 /* 4 points DFT */
 if(m0>=2)
 for(li=0;li<np;li+=4)
 {
  j=li;k=j+2;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];
  sx[j]+=sx[k];
  sr[k]=t1;sx[k]=t2;
  ++j;++k;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];sx[j]+=sx[k];
  sr[k]=inv*t2;sx[k]=-inv*t1;
 }
 /* 2 points DFT */
 for(li=0;li<np;li+=2)
 {
  j=li;k=j+1;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];sx[j]+=sx[k];
  sr[k]=t1;sx[k]=t2;
 }
 /* sort according to bit reversal */
 lmx=np/2;j=0;
 for(i=1;i<np-1;++i)
 {
  k=lmx;
  while(k<=j)
  {
   j-=k;k/=2;
  }
  j+=k;
  if(i<j)
  {
   t1=sr[j];sr[j]=sr[i];sr[i]=t1;
   t1=sx[j];sx[j]=sx[i];sx[i]=t1;
  }
 }
 /*  if Inverse FFT, multiply 1.0/np  */
 if(inv!=-1)
  return;
 t1=1.0/np;
 for(i=0;i<np;++i)
 {
  sr[i]*=t1;sx[i]*=t1;
 }
}
//dft_sub.c

 

void dft(float xr[],float xi[],int flag)
//float xr[N],xi[N];
//int flag;
{ float XR[N],XI[N];
  int k,n;
  float sum1,sum2,cita;
  for(k=0;k<=N-1;k++)
  { sum1=0.0;
 sum2=0.0;
    { for(n=0;n<=N-1;n++)
     { cita=2.0*3.1415926/N*n*k;
   sum1=sum1+xr[n]*cos(cita)+flag*xi[n]*sin(cita);
   sum2=sum2-flag*xr[n]*sin(cita)+xi[n]*cos(cita);
  }
   XR[k]=sum1;
      XI[k]=sum2;
 }
  }
  if(flag==1)
   for(k=0;k<=N-1;k++)
  {
   xr[k]=XR[k];
   xi[k]=XI[k];
  }
  else
   for(k=0;k<=N-1;k++)
  {
   xr[k]=XR[k]/N;
   xi[k]=XI[k]/N;
  }
}

显然，FFT比DFT运算效率高得多。

如果程序无法自己编辑出，可以以N=8为例理解验证程序的正确性。