快速傅里叶变换FFT

大青梅

求用vc6.0B编写FFT的程序，新人求指导

TestGuru

原著：（FORTRAN)：

TestGuru

第一改编：（VC6.0）
/*******************************************************
/*             FFT                                     *
/* xr,xi-----input & output array
/* n-----must be integer exponent of 2                 *
/* inv---0:FFT   1---IFFT note the 1/N difference      *
/* arry of fft, e.g. xr,xi, is started with 1,not 0    *
/*******************************************************/
/*
   int FFT(float * xr,float * xi,long n,int inv)  
   {
         float tr,ti,ur,ui,wr,wi,ur1,ui1;
     float pi;
     long m,nv2,nm1,i,j,k,l,le,le1,ip;
         //  n point fft  
     m=log(n)/log(2)+0.1;
         nv2=n/2;
     nm1=n-1;
     j=1;
     for(i=1;i<=nm1;i++)
     {  
                 if(i>=j) goto lop10;
                 tr=xr[j];
                 ti=xi[j];
                 xr[j]=xr[i];
                 xi[j]=xi[i];
                 xr[i]=tr;
                 xi[i]=ti;
                 lop10:        k=nv2;
                 lop20:        if(k>=j) goto lop30;
                 j=j-k;
                 k=k/2;
                 goto lop20;
                 lop30:        j=j+k;
      }
          pi=4.0*atan(1.0);
          for(l=1;l<=m;l++)
          {
                  le=pow(2,l)+0.1;
              le1=le/2;
              ur=1.0;
                  ui=0.0;
                  wr=cos(pi/(float)le1);
              wi=-sin(pi/(float)le1);
              if(inv!=0) wi=-wi;
                  for(j=1;j<=le1;j++)
                        {
                          for(i=j;i<=n;i+=le)
                                {
                                  ip=i+le1;
                                  tr=xr[ip]*ur-xi[ip]*ui;
                                  ti=xr[ip]*ui+xi[ip]*ur;
                                  xr[ip]=xr[i]-tr;
                                  xi[ip]=xi[i]-ti;
                                  xr[i]=xr[i]+tr;
                                  xi[i]=xi[i]+ti;
                                }
                    ur1=ur*wr-ui*wi;
                   ui1=ur*wi+ui*wr;
                   ur=ur1;
                   ui=ui1;
                        }
          }
          if(inv==0) return 0;
          for(i=1;i<=n;i++)
          {
            xr[i]=xr[i]/(float)n;
            xi[i]=xi[i]/(float)n;
          }
          return 0;
   }

*/

TestGuru

第二改编：（VC6.0）
/*******************************************************
/*             FFT                                     *
/* xr,xi-----input & output array
/* n-----must be integer exponent of 2                 *
/* inv---0:FFT   1---IFFT note the 1/N difference      *
/* arry of fft, e.g. xr,xi, is started with 0,not 1    *
/*******************************************************/

   int FFT(float * xr,float * xi,long n,int inv)  
   {
         float tr,ti,ur,ui,wr,wi,ur1,ui1;
     float pi;
     long m,nv2,nm1,i,j,k,l,le,le1,ip;
         //  n point fft  
     m=log(n)/log(2)+0.1;
         nv2=n/2;
     nm1=n-1;
     j=0;//j=1;
     for(i=0;i<nm1;i++) //for(i=1;i<=nm1;i++)
     {  
                 if(i>=j) goto lop10;
                 tr=xr[j];
                 ti=xi[j];
                 xr[j]=xr[i];
                 xi[j]=xi[i];
                 xr[i]=tr;
                 xi[i]=ti;
                 lop10:        k=nv2;
                 lop20:        if(k>j) goto lop30;//if(k>=j) goto lop30;
                 j=j-k;
                 k=k/2;
                 goto lop20;
                 lop30:        j=j+k;
      }
          pi=4.0*atan(1.0);
          for(l=0;l<m;l++)  //for(l=1;l<=m;l++)  
          {
                  le=pow(2,l+1)+0.1;//le=pow(2,l)+0.1;
              le1=le/2;
              ur=1.0;
                  ui=0.0;
                  wr=cos(pi/(float)le1);
              wi=-sin(pi/(float)le1);
              if(inv!=0) wi=-wi;
                  for(j=0;j<le1;j++)  //for(j=1;j<=le1;j++)
                        {
                          for(i=j;i<n;i+=le)  // for(i=j;i<=n;i+=le)
                                {
                                  ip=i+le1;
                                  tr=xr[ip]*ur-xi[ip]*ui;
                                  ti=xr[ip]*ui+xi[ip]*ur;
                                  xr[ip]=xr[i]-tr;
                                  xi[ip]=xi[i]-ti;
                                  xr[i]=xr[i]+tr;
                                  xi[i]=xi[i]+ti;
                                }
                    ur1=ur*wr-ui*wi;
                   ui1=ur*wi+ui*wr;
                   ur=ur1;
                   ui=ui1;
                        }
          }
          if(inv==0) return 0;
          for(i=0;i<n;i++) //for(i=1;i<=n;i++)
          {
            xr[i]=xr[i]/(float)n;
            xi[i]=xi[i]/(float)n;
          }
          return 0;
   }