(完整word版)fft算法代码注释及流程框图

基2的DIT蝶形算法源代码及注释如下：/************FFT***********///整个程序输入和输出利用同一个空间x[N]，节约空间#includestdio.h#includemath.h#includestdlib.h#defineN1000//定义输入或者输出空间的最大长度typedefstruct{doublereal;doubleimg;}complex;//定义复数型变量的结构体voidfft();//快速傅里叶变换函数声明voidinitW();//计算W(0)~W(size_x-1)的值函数声明voidchange();//码元位置倒置函数函数声明voidadd(complex,complex,complex*);/*复数加法*/voidmul(complex,complex,complex*);/*复数乘法*/voidsub(complex,complex,complex*);/*复数减法*/voiddivi(complex,complex,complex*);/*复数除法*/voidoutput();/*输出结果*/complexx[N],*W;/*输出序列的值*/intsize_x=0;/*输入序列的长度，只限2的N次方*/doublePI;//pi的值intmain(){inti;system(cls);PI=atan(1)*4;printf(Pleaseinputthesizeofx:\n);/*输入序列的长度*/scanf(%d,&size_x);printf(Pleaseinputthedatainx[N]:(suchas:56)\n);/*输入序列对应的值*/for(i=0;isize_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的结果return0;}/*进行基-2FFT运算,蝶形算法。这个算法的思路就是，先把计算过程分为log(size_x)/log(2)-1级（用i控制级数）；然后把每一级蝶形单元分组（用j控制组的第一个元素起始下标）；最后算出某一级某一组每一个蝶形单元（用k控制个数，共l个）。*/voidfft(){inti=0,j=0,k=0,l=0;complexup,down,product;change();//实现对码位的倒置for(i=0;ilog(size_x)/log(2);i++)//循环算出fft的结果{l=1i;for(j=0;jsize_x;j+=2*l)//算出第m=i级的结果【i从0到（log(size_x)/log(2)）-1】{for(k=0;kl;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为蝶形单元右下方的值}}}}voidinitW()//计算W的实现函数{inti;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的虚部}}voidchange()//输入的码组码位倒置实现函数{complextemp;unsignedshorti=0,j=0,k=0;doublet;for(i=0;isize_x;i++){k=i;j=0;t=(log(size_x)/log(2));while((t--)0){j=j1;j|=(k&1);k=k1;}if(ji){temp=x[i];x[i]=x[j];x[j]=temp;}}}voidoutput()//输出结果实现函数{inti;printf(Theresultareasfollows\n);for(i=0;isize_x;i++){printf(%.4f,x[i].real);//输出实部if(x[i].img=0.0001)//如果虚部的值大于0.0001，输出+jx.img的形式printf(+j%.4f\n,x[i].img);elseif(fabs(x[i].img)0.0001)//如果虚部的绝对值小于0.0001，无需输出printf(\n);elseprintf(-j%.4f\n,fabs(x[i].img));//如果虚部的值小于-0.0001，输出-jx.img的形式}}voidadd(complexa,complexb,complex*c)//复数加法实现函数{c-real=a.real+b.real;//复数实部相加c-img=a.img+b.img;//复数虚部相加}voidmul(complexa,complexb,complex*c)//复数乘法实现函数{c-real=a.real*b.real-a.img*b.img;//获取相乘结果的实部c-img=a.real*b.img+a.img*b.real;//获取相乘结果的虚部}voidsub(complexa,complexb,complex*c)//复数减法实现函数{c-real=a.real-b.real;//复数实部相减c-img=a.img-b.img;//复数虚部相减}voiddivi(complexa,complexb,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，输入53）计算出前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]lX[j+k+l]*W[(size_x/2/l)*k]X[j+k+l]-1K加1