下面是从书本的FORTRAN FFT程序转为的C程序：
1）原书的FORTRAN程序：
2）第一次改编：（VC6.0），由于FORTRAN的数组从1开始，所以这里的C的数组也是从1开始的
/*******************************************************
/* 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;
xi[j]=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[ip]*ur-xi[ip]*ui;
ti=xr[ip]*ui+xi[ip]*ur;
xr[ip]=xr-tr;
xi[ip]=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;
}
3）第二次改编：（VC6.0），这里的C的数组改为更常见的从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;
xi[j]=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[ip]*ur-xi[ip]*ui;
ti=xr[ip]*ui+xi[ip]*ur;
xr[ip]=xr-tr;
xi[ip]=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;
}
|
