大青梅 发表于 2015-5-19 11:30

快速傅里叶变换FFT

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

TestGuru 发表于 2015-5-26 10:44

原著：（FORTRAN)：


TestGuru 发表于 2015-5-26 10:45

第一改编：（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;
               ti=xi;
               xr=xr;
               xi=xi;
               xr=tr;
               xi=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*ur-xi*ui;
                                  ti=xr*ui+xi*ur;
                                  xr=xr-tr;
                                  xi=xi-ti;
                                  xr=xr+tr;
                                  xi=xi+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=xr/(float)n;
          xi=xi/(float)n;
          }
          return 0;
   }

*/

TestGuru 发表于 2015-5-26 10:46

第二改编：（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;
               ti=xi;
               xr=xr;
               xi=xi;
               xr=tr;
               xi=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*ur-xi*ui;
                                  ti=xr*ui+xi*ur;
                                  xr=xr-tr;
                                  xi=xi-ti;
                                  xr=xr+tr;
                                  xi=xi+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=xr/(float)n;
          xi=xi/(float)n;
          }
          return 0;
   }
页: [1]
查看完整版本: 快速傅里叶变换FFT