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

基2的DIT蝶形算法源代码及注释如下： /*FFT*/ /整个程序输入和输出利用同一个空间xN，节约空间 #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 xN,*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 xN:(such as:5 6)n); /*输入序列对应的值*/ for(i=0;isize_x;i+) scanf(%lf %lf,&xi.real,&xi.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;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(xj+k+l,W(size_x/2/l)*k,&product); /*size/2/l是该级W的相邻上标差，l是该级该组取的W总个数*/ add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; /up为蝶形单元右上方的值 xj+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个值即可*/ Wi.real=cos(2*PI/size_x*i); /计算W的实部 Wi.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;i0) j=j1; if(ji) temp=xi; xi=xj; xj=temp; void output() /输出结果实现函数 int i; printf(The result are as followsn); for(i=0;i=0.0001) /如果虚部的值大于0.0001，输出+jx.img的形式 printf(+j%.4fn,xi.img); else if(fabs(xi.img)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/2个exp(-j*2*k/size_x)个值，即W的值输入序列对应值（例如5+j3，输入5 3）输入序列长度size_x飞级