fft算法代码注释及流程框图

基2的DIT蝶形算法源代码及注释如下：
/************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1000 //定义输入或者输出空间的最大长度
typedef struct
{
double real;
double img;
}complex; //定义复数型变量的结构体
void fft(); //快速傅里叶变换函数声明
void initW(); //计算W(0)~W(size_x-1)的值函数声明
void change(); //码元位置倒置函数函数声明
void add(complex,complex,complex *); /*复数加法*/
void mul(complex,complex,complex *); /*复数乘法*/
void sub(complex,complex,complex *); /*复数减法*/
void divi(complex,complex,complex *); /*复数除法*/
void output(); /*输出结果*/
complex x[N],*W; /*输出序列的值*/
int size_x=0; /*输入序列的长度，只限2的N次方*/
double PI; //pi的值
int main()
{
int i;
system("cls");
PI=atan(1)*4;
printf("Please input the size of x:\n"); /*输入序列的长度*/
scanf("%d",&size_x);
printf("Please input the data in x[N]:(such as:5 6)\n"); /*输入序列对应的值*/
for(i=0;i<size_x;i++)
scanf("%lf %lf",&x[i].real,&x[i].img);
initW(); //计算W(0)~W(size_x-1)的值
fft(); //利用fft快速算法进行DFT变化output(); //顺序输出size_x个fft的结果
return 0;
}
/*进行基-2 FFT运算,蝶形算法。

这个算法的思路就是，
先把计算过程分为log(size_x)/log(2)-1级（用i控制级数）；
然后把每一级蝶形单元分组（用j控制组的第一个元素起始下标）；
最后算出某一级某一组每一个蝶形单元（用k控制个数，共l个）。

*/
void fft()
{
int i=0,j=0,k=0,l=0;
complex up,down,product;
change(); //实现对码位的倒置
for(i=0;i<log(size_x)/log(2);i++) //循环算出fft的结果
{
l=1<<i;
for(j=0;j<size_x;j+=2*l)
//算出第m=i级的结果【i从0到（log(size_x)/log(2)）-1】{
for(k=0;k<l;k++) //算出第i级内j组蝶形单元的结果
{ //算出j组中第k个蝶形单元mul(x[j+k+l],W[(size_x/2/l)*k],&product); /*size/2/l是该级W的相邻上标差，
l是该级该组取的W总个数*/ add(x[j+k],product,&up);
sub(x[j+k],product,&down);
x[j+k]=up; //up为蝶形单元右上方的值
x[j+k+l]=down; //down为蝶形单元右下方的值}
}
}
}
void initW() //计算W的实现函数
{
int i;
W=(complex *)malloc(sizeof(complex) * size_x); /*申请size_x个复数W的空间(这部
申请的空间有点多，实际上只要
申请size_x/2个即可)*/ for(i=0;i<(size_x/2);i++) /*预先计算出size_x/2个W的值，存
放,由于蝶形算法只需要前
size_x/2个值即可*/ {
W[i].real=cos(2*PI/size_x*i); //计算W的实部
W[i].img=-1*sin(2*PI/size_x*i); //计算W的虚部
}
}
void change() //输入的码组码位倒置实现函数{
complex temp;
unsigned short i=0,j=0,k=0;
double t;
for(i=0;i<size_x;i++)
{
k=i;
j=0;
t=(log(size_x)/log(2));
while((t--)>0)
{
j=j<<1;
j|=(k & 1);
k=k>>1;
}
if(j>i)
{
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
void output() //输出结果实现函数{
int i;
printf("The result are as follows\n");
for(i=0;i<size_x;i++)
{
printf("%.4f",x[i].real); //输出实部
if(x[i].img>=0.0001) //如果虚部的值大于0.0001，输出+jx.img的形式printf("+j%.4f\n",x[i].img);
else if(fabs(x[i].img)<0.0001) //如果虚部的绝对值小于0.0001，无需输出
printf("\n");
else
printf("-j%.4f\n",fabs(x[i].img));
//如果虚部的值小于-0.0001，输出-jx.img的形式}
}
void add(complex a,complex b,complex *c) //复数加法实现函数{
c->real = a.real + b.real; //复数实部相加
c->img = a.img + b.img; //复数虚部相加
}
void mul(complex a,complex b,complex *c) //复数乘法实现函数
{
c->real = a.real*b.real - a.img*b.img; //获取相乘结果的实部
c->img = a.real*b.img + a.img*b.real; //获取相乘结果的虚部
}
void sub(complex a,complex b,complex *c) //复数减法实现函数
{
c->real = a.real - b.real; //复数实部相减
c->img = a.img - b.img; //复数虚部相减
}
void divi(complex a,complex b,complex *c) //复数除法实现函数
{
c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real+b.img*b.img);
//获取相除结果的实部c->img = (a.img*b.real - a.real*b.img) / (b.real*b.real+b.img*b.img);
//获取相除结果的虚部
}
由于输入序列和输出序列存储（包括中间结果）在同一个空间数组x【size_x】，故每进行一个蝶形算法，相应的存储单元有两个值被替换。

该程序框图如下：
飞
是
否
是
否
是
否
输入序列长度size_x
输入序列对应值（例如5+j3，输入5 3）
计算出前size_x/2个 exp(-j*2π*k/size_x)个值，
即W 的值
级数i>=
？
输出fft 结果序列 结束
计算出该级需要的
W 的个数l
该级该组起始下标j>=？
级数i 加1
组起始下标加2*l
该级该组元素序数k>=？
X[j+k] X[j+k]l
X[j+k+l]*W[(size_x/2/l)*k] X[j+k+l] -1
K 加1。