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

基2的DIT蝶形算法源代码及注释如下：

/************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间

#include

#include

#include

#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

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

{

l=1<

for(j=0;j

//算出第m=i级的结果【i从0到（log(size_x)/log(2)）-1】{

for(k=0;k

{ //算出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

{

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

{

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的形式}

}