然后还需判断次高位，这可与k=N\4相比较，若次高位为0，则需将它变为1<加N\4即可）其他位不变，既得到下一个倒位序数；若次高位是1，则需将它也变为0。
然后再判断下一位
。
2.复数运算
因为每一个蝶形结构完成的迭代运算为
算式中涉及到了复数的运算，而计算机是不能自己实现复数运算的，所以需要我们自己设计进行复数运算的程序。
迭代运算式中，=cos<2πr/N）-jsin<2πr/N）我们设=R（j>+jI（j>,=R（K>+jI（K>,
而我们最后希望得到的DFT结果是复数的模，根据它的模来绘制频谱，所以这里我们定义
X（k>=
相关程序我们编译为：
c.real=a.real*b.real-a.imag*b.imag。
c.imag=a.real*b.imag+a.imag*b.real。
根据迭代运算的式子，我们可以将其分解为：
R（K>+jI（K>=R（K>+jI（K>+[R（j>+jI（j>]*[cos<2πr/N）-jsin<2πr/N）]。
继续分解得到下列两式:
R（K>=R（K>+R（j>cos<2πr/N）+I（j>sin（2πr/N>I（K>=I（K>-R（j>sin<2πr/N）+I（j>cos（2πr/N>
同理
R（j>=R（K>-R（j>cos<2πr/N）-I（j>sin（2πr/N>I（j>=I（K>+R（j>sin<2πr/N）I（j>cos（2πr/N>
相关程序编译为：
xin[ip].real=xin[i].real-t.real。
xin[ip].imag=xin[i].imag-t.imag。
xin[i].real=xin[i].real+t.real。
xin[i].imag=xin[i].imag+t.imag。
3.节点距离运算
当输入为倒位序，输出为正常顺序时，第m级运算，每个蝶形的两节点距离为。
4.旋转因子的计算
这里主要解决的是迭代运算的每一项应乘的旋转因为中r应取多少的问题。
第L级的2^（L-1>个碟形因子WPN中的P，可表示为p=j*2^（m-L>，其中j=0，1，2，...，<2L-1-1）。
将上述方法整理一下，我们可以得到一个完整的程序的思路。
主要分为三个模块，复数计算DFT模函数，FFT函数，主函数。
其中复数计算DFT模函数的算法我们已经在前面介绍了。
主函数是编译整个程序进行的顺序，即我们先设定一个时域函数，从中抽取N点的值，然后进行FFT运算，然后求模长，最后将对应的频域的值输出。
这里比较复杂的是FFT函数。
下面进行解释。
这里我们共设三个循环来实现FFT这一函数的功能。
第一层：
首先我们知道，一个N=2^m点DFT，要进行m级计算,所以第一层循环是对运算的级数进行控制。
第二层：
因为第L级有2^（L-1>个蝶形因子，第二层循环根据乘数进行控制，保证对于每一个旋转因子第三层循环要执行一次，这样，第三层循环在第二层循环控制下，每一级要进行2^（L-1>次循环计算。
第三层：
因为第L级共有N/2^L个蝶形结构，并且同一级内不同蝶型结构的旋转因子分布相同，当第二层循环确定某一旋转因子后，第三层循环要将本级中每个蝶型结构中具有这一旋转引自的蝶形计算一次，即第三层循环每执行完一次要进行N/2^L个蝶形计算。
所以，在每一级中，第三层循环完成N/2^L个蝶形计算；第二层循环使得第三层循环进行2^（L-1>次，因此，第二层循环完成时，共进行2^（L-1>*N/2^L=N/2个碟形计算。
实质是：
第二、第三层循环完成了第L级的计算。
五、程序代码
1.C语言：
#include
#include
#include
#definePI3.14932384626433832795028841971
#defineFFT_N128//定义傅利叶变换的点数
structcompx{doublereal,imag。
}。
//定义一个复数结构
structcompxs[FFT_N]。
//从S[1]开始存放
structcompxEE（structcompxa,structcompxb>//求复数的模长函数
{
structcompxc。
c.real=a.real*b.real-a.imag*b.imag。
c.imag=a.real*b.imag+a.imag*b.real。
return（c>。
}
voidFFT（structcompx*xin>//FFT函数
{
intf,m,nv2,nm1,i,k,l,j=0。
structcompxu,w,t。
nv2=FFT_N/2。
//变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1。
for（i=0。
ii++>
{
if（i//如果i{
t=xin[j]。
xin[j]=xin[i]。
xin[i]=t。
}
k=nv2。
//求j的下一个倒位序
while（k<=j>//如果k<=j,表示j的最高位为1
{
j=j-k。
//把最高位变成0
k=k/2。
//比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k。
//把0改为1
}
{
intd,d1,ip。
f=FFT_N。
for（l=1。
（f=f/2>!
=1。
l++>//第一层循环，计算蝶形级数
。
for（m=1。
m<=l。
m++>//第二层循环，控制蝶形级数
{
d=2<<（m-1>。
//d蝶形结距离，即第m级蝶形的蝶形结构相距d点
d1=d/2。
//同一蝶形结中参加运算的两点的距离
u.real=1.0。
//u为蝶形结构运算系数，初始值为1
u.imag=0.0。
w.real=cos（PI/d1>。
//w为系数商，即当前系数与前一个系数的商
w.imag=-sin（PI/d1>。
for（j=0。
j<=d1-1。
j++>//第三层循环，进行旋转因子的计算完成蝶形计算。
这里控制计算系数不同的蝶形结
{
for（i=j。
i<=FFT_N-1。
i=i+d>//控制计算系数相同蝶形结
{
ip=i+d1。
//i，ip分别表示参加蝶形运算的两个节点
t=EE（xin[ip],u>。
//蝶形运算
xin[ip].real=xin[i].real-t.real。
xin[ip].imag=xin[i].imag-t.imag。
xin[i].real=xin[i].real+t.real。
xin[i].imag=xin[i].imag+t.imag。
}
u=EE（u,w>。
//改变系数，进行下一个蝶形运算
}
}
}
}
voidmain（>
{
inti。
for（i=0。
ii++>//给结构体赋值
{
s[i].real=sin（2*3.1493*i/FFT_N>。
s[i].imag=0。
}
longstart=clock（>。
FFT（s>。
for（i=0。
ii++>
{s[i].real=sqrt（s[i].real*s[i].real+s[i].imag*s[i].imag>。
longend=clock（>。
printf（"%.4f\n",s[i].real>。
}
longt=end-start。
printf（"%d\n",t>。
}
2.MATLAB：
n=0:
127。
b=sin（2*pi*n/128>。
fft（b>
六、实验结果
1.C语言
2.MATLAB
Columns1through5
0.0000-0.0000-64.0000i-0.0000-0.0000i-0.0000
升级会员 