void FFT1(float a_r1[],float a_im[],int ex,bool inv)
{
int i,length=1;
float *sin_tbl;
float *cos_tbl;
float *buf;
for(i=0;i<ex;i++)
{
length*=2;
}
sin_tbl=(float *)malloc(length*sizeof(float));
cos_tbl=(float *)malloc(length*sizeof(float));
buf=(float *)malloc(length*sizeof(float));
cstb(length,inv,sin_tbl,cos_tbl);
fft1core(a_r1,a_im,length,ex,sin_tbl,cos_tbl,buf);
free(sin_tbl);
free(cos_tbl);
free(buf);
}
二維的：
void FFT2(float *a_rl,float *a_im,bool inv,int xsize,int ysize)
{
float *b_rl,*b_im;
float *hsin_tbl,*hcos_tbl,*vsin_tbl,*vcos_tbl;
float *buf_x,*buf_y;
int i;
b_rl=(float *)calloc(xsize*ysize,sizeof(float));
b_im=(float *)calloc(xsize*ysize,sizeof(float));
hsin_tbl=(float *)calloc(xsize,sizeof(float));
hcos_tbl=(float *)calloc(xsize,sizeof(float));
vsin_tbl=(float *)calloc(ysize,sizeof(float));
vcos_tbl=(float *)calloc(ysize,sizeof(float));
buf_x=(float *)malloc(xsize*sizeof(float));
buf_y=(float *)malloc(ysize*sizeof(float));
cstb(xsize,inv,hsin_tbl,hcos_tbl);
cstb(ysize,inv,vsin_tbl,vcos_tbl);
int x_exp=(int)(log(double(xsize))/log(2.0));
for(i=0;i<ysize;i++)
{
fft1core(&(*(a_rl+i*xsize)),&(*(a_im+i*xsize)),xsize,x_exp,hsin_tbl,hcos_tbl,buf_x);
}
rvmtx1(a_rl,b_rl,xsize,ysize);
rvmtx1(a_im,b_im,xsize,ysize);
int y_exp=(int)(log(double(ysize))/log(2.0));
for(i=0;i<xsize;i++)
{
fft1core(&(*(b_rl+i*ysize)),&(*(b_im+i*ysize)),ysize,y_exp,vsin_tbl,vcos_tbl,buf_y);
}
rvmtx2(b_rl,a_rl,xsize,ysize);
rvmtx2(b_im,a_im,xsize,ysize);
free(b_rl);
free(b_im);
free(hsin_tbl);
free(hcos_tbl);
free(vsin_tbl);
free(vcos_tbl);
free(buf_x);
free(buf_y);
}
用到的函數：
void fft1core(float a_rl[],float a_im[],int length,int ex,float sin_tbl[],float cos_tbl[],float buf[])
{
int i,j,k,w,j1,j2;
int numb,lenb,timb;
float xr,xi,yr,yi,nrml;
for(i=1;i<length;i+=2)
{
a_rl[i]=-a_rl[i];
a_im[i]=-a_im[i];
}
numb=1;
lenb=length;
for(i=0;i<ex;i++)
{
lenb/=2;
timb=0;
for(j=0;j<numb;j++)
{
w=0;
for(k=0;k<lenb;k++)
{
j1=timb+k;
j2=j1+lenb;
xr=a_rl[j1];
xi=a_im[j1];
yr=a_rl[j2];
yi=a_im[j2];
a_rl[j1]=xr+yr;
a_im[j1]=xi+yi;
xr=xr-yr;
xi=xi-yi;
a_rl[j2]=xr*cos_tbl[w]-xi*sin_tbl[w];
a_im[j2]=xr*sin_tbl[w]+xi*cos_tbl[w];
w+=numb;
}
timb+=(2*lenb);
}
numb*=2;
}
birv(a_rl,length,ex,buf);
birv(a_im,length,ex,buf);
for(i=1;i<length;i+=2)
{
a_rl[i]=-a_rl[i];
a_im[i]=-a_im[i];
}
nrml=1/sqrt(float(length));
for(i=0;i<length;i++)
{
a_rl[i]*=nrml;
a_im[i]*=nrml;
}
}
void cstb(int length,bool inv,float sin_tbl[],float cos_tbl[])
{
int i;
float xx,arg;
xx=float((-PI)*2/length);
if(!inv)
{
xx=-xx;
}
for(i=0;i<length;i++)
{
arg=i*xx;
sin_tbl[i]=sin(arg);
cos_tbl[i]=cos(arg);
}
}
void birv(float a[],int length,int ex,float b[])
{
int i,ii,k,bit;
for(i=0;i<length;i++)
{
for(k=0,ii=i,bit=0;;bit<<=1,ii>>=1)
{
bit=(ii&1)|bit;
if(++k==ex)
{
break;
}
}
b[i]=a[bit];
}
for(i=0;i<length;i++)
{
a[i]=b[i];
}
}
void rvmtx1(float *a,float *b,int xsize,int ysize)
{
int i,j;
for(j=0;j<ysize;j++)
{
for(i=0;i<xsize;i++)
{
*(b+i*ysize+j)=*(a+j*xsize+i);
}
}
}
void rvmtx2(float *a,float *b,int xsize,int ysize)
{
int i,j;
for(j=0;j<ysize;j++)
{
for(i=0;i<xsize;i++)
{
*(b+j*xsize+i)=*(a+i*ysize+j);
}
}
}