完整的快速傅立叶变换和逆变换（FFT IFFT）C语言算法

函数语句：

void kbfft(pr,pi,n,k,fr,fi,l,il)

pr－－－－-双精度实型一维数组，长度为n。当n=0时，存放n个采样输入的实部，返回时存放离散傅立叶变换的模；当l＝1 时，存放傅立叶变换的n个实部，返回时存放逆傅立叶变换的模。

pi－－－－-双精度实型一维数组，长度为n。当n=0时，存放n个采样输入的虚部，返回时存放离散傅立叶变换的幅角；当l＝ 1时，存放傅立叶变换的n个虚部，返回时存放逆傅立叶变换的幅角。其中幅角的单位为度。

n－－－－－－整形变量。输入的点数。

k－－－－－－整形变量。满足n=2~k（2的k次方）。

fr－－－－-双精度实型一维数组，长度为n。当l=0时，返回傅立叶变换的实部；当l＝1时，返回逆傅立叶变换的实部。

fi－－－－-双精度实型一维数组，长度为n。当l=0时，返回傅立叶变换的虚部；当l＝1时，返回逆傅立叶变换的虚部。

l－－－－－－整形变量。若l＝0，表示要求本函数计算傅立叶变换；若l＝1，表示要求本函数计算逆傅立叶变换。

il－－－－-整形变量。若il＝0，表示不要求本函数计算傅立叶变换或逆傅立叶变换的模与幅角；若il＝1，表示要求本函数 计算傅立叶变换或逆傅立叶变换的模与幅角。

＃i nclude"math.h"

void kbfft(pr,pi,n,k,fr,fi,l,il)

int n,k,l,il;
double pr【】,pi【】,fr【】,fi【】;

{int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;

for(it=0;it<=n-1;it++)
{m=it;is=0;
for(i=0;i<=k-1;i++)
{j=m/2;is=2*is+(m-2*j);m=j;}
fr【it】=pr【is】;fi【it】=pi【is】;
}

pr【0】=1.0;pi【0】=0.0;
p=6.283185306/(1.0*n);
pr【1】=cos(p);pi【1】=-sin(p);

if(l!=0) pi【1】=-pi【1】;

for(i=2;i<n-1;i++)
{p=pr【i-1】*pr【1】;q=pi【i-1】*pi【1】;
s=(pr【i-1】+pi【i-1】)*(pr【1】+pi【1】);
pr【i】=p-q;pi【i】=s-p-q;
}

for(it=0;it<=n-2;it=it+2)
{vr=fr【it】;vi=fi【it】;
fr【it】=vr+fr【it+1】;fi【it】=vi+fi【it+1】;
fr【it+1】=vr-fr【it+1】;fi【it+1】=vi-fi【it+1】;
}

m=n/2;nv=2;

for(l0=k-2;l0>=0;l0－－)
{m=m/2;nv=2*nv;
for(it=0;it<=(m-1)*nv;it=it+nv)
for(j=0;j<=(nv/2)-1;j++)
{p=pr【m*j】*fr【it+j+nv/2】;
q=pi【m*j】*fi【it+j+nv/2】;
s=pr【m*j】+pi【m*j】;
s=s*(fr【it+j+nv/2】+fi【it+j+nv/2】);
poddr=p-q;poddi=s-p-q;
fr【it+j+nv/2】=fr【it+j】-poddr;
fi【it+j+nv/2】=fi【it+j】-poddi;
fr【it+j】=fr【it+j】+poddr;
fi【it+j】=fi【it+j】+poddi;
}
}

if(l!=0)
for(i=0;i<=n-1;i++)
{fr【i】=fr【i】/(1.0*n);
fi【i】=fi【i】/(1.0*n);
}

if(il!=0)
for(i=0;i<=n-1;i++)
{pr【i】=sqrt(fr【i】*fr【i】+fi【i】*fi【i】);
if(fabs(fr【i】)<0.000001*fabs(fi【i】))
{if((fi【i】*fr【i】)>0) pi【i】=90.0;
else pi【i】=-90.0;
}
else
pi【i】=atan(fi【i】/fr【i】)*360.0/6.283185306;
}

return;
}