1、串行 FFT递归算法（蝶式递归计算原理）求傅里叶变换摘要FFT，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设 x(n)为 N 项的复数序列，由 DFT 变换，任一 X（m）的计算都需要 N 次复数乘法和 N-1 次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算” （四次实数乘法和四次实数加法） ，那么求出 N 项复
2、数序列的 X（m）,即 N 点 DFT 变换大约就需要 N2 次运算。当 N=1024 点甚至更多的时候，需要 N2=1048576 次运算，在 FFT 中，利用 WN 的周期性和对称性，把一个 N 项序列（设 N=2k,k 为正整数） ，分为两个 N/2 项的子序列，每个 N/2 点 DFT 变换需要（N/2）2 次运算，再用 N 次运算把两个 N/2 点的DFT 变换组合成一个 N 点的 DFT 变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024 时，总的运算次数就变成了 525312次，节省了大约 50%的运算量。而如果我们将这种“一分为
3、二”的思想不断进行下去，直到分成两两一组的 DFT 运算单元，那么 N 点的 DFT 变换就只需要 Nlog(2)(N)次的运算，N 在 1024 点时，运算量仅有 10240 次，是先前的直接算法的 1%，点数越多，运算量的节约就越大，这就是 FFT 的优越性。关键字：FFT 蝶式计算 傅里叶变换目录一题目及要求 .11.1题目 .1二设计算法、算法原理 .12.1算法原理与设计 .12.2设计步骤 .2三算法描述、设计流程 .43.1算法描述 .43.2流程图 .5四 源程序代码及运行结果 .74.1源程序代码 .74.2运行结果 .12五 算法分析、优缺点 .135.1算法分析 .135
4、.2优缺点 .14六 总结 .15七 参考文献 .161一题目及要求1.1题目对给定的 ，利用串行 FFT 递归算法（蝶式递归计算)23,7168,30742,1(原理）计算其傅里叶变换的结果。1.2要求利用串行递归与蝶式递归原理，对给定的向量求解傅里叶变换的结果。二设计算法、算法原理2.1算法原理与设计蝶式递归计算原理：令 为 n/2 次单位元根，则有 ，将 b 向量的偶数项 和奇数项 分别记 为 和 。注意推导中反复使用： 。图 2.1 公式图形)2/(2nie 2Tnb),.(220Tnb),.(131 Tnb),.(1102n,.(02 psnn ,ln/22.2设计步骤 DFTaab
5、nTn nnn 的是因 此 ， 向 量 ),.,(),.( 11020 222 DFTaaabnTn nnn 的是因 此 ， 向 量 )(),.(),(),.( 1110131 222 对于以上的分析可画出如图 2.2 所示的离散傅里叶变换递归计算流图。图 2.3 就是一个按此递归方法计算的 n=8 的 FFT 蝶式计算图。 1,0)( )()( )()( 2111 220 112 22420 12412 1201222 22222222 nkkkl nl llnl ll nlll lllnkklll laaaaaaabn nnn nnnnnnn 偶 数 时 ： 1,0)( )()()()(
6、210 12210 12421121 112420 )(1)2(1)1( 1)2(1)(201)(22 2222 222 22 2 nkkkl lll nlll nll lll nlnll lllkklll laaaaaaaabnn nnnnn nnnnn nnn nnnn nn 奇 数 时 ：3FFT 的蝶式递归计算图（有计算原理推出）：图 2.2 递归计算流图特别的，的 FFT 蝶式计算图（展开的）：图 2.3 蝶式计算图按输入元素 展开，前面将输出序列之元素 按其偶下标（ ）和（kajblj2）展开，导出 和 递归计算式，按此构12lj 1022)(nnkkla1022)(nnkkkla
7、造出了如图 1 所示的 FFT 递归计算流程图。事实上，我们也可以将输入序列之元素.a0a1an-1DFT(n/2)DFT(n/2).+-1 -a (a)n2-1n21n-1-a () n2+11-an20+aan2-1n-1+an20+an2+11an2-1an2an2+1 . n2-14按其偶下标（ ） 和几下标（ ）展开，则导出另一种形式的 FFT 递归计算式 kak212k。10nkkjjwb三算法描述、设计流程3.1算法描述SISD 上的 FFT 分治递归算法：输入: a=(a 0,a1,an-1); 输出: B=(b 0,b1,bn-1)Procedure RFFT(a,b)beg
8、inif n=1 then b0=a0 else(1)RFFT(a0,a2,an-2, u0,u1,un/2-1)(2)RFFT(a1,a3,an-1, v0,v1,vn/2-1)(3)z=1(4)for j=0 to n-1 do(4.1)bj=uj mod n/2+zvj mod n/2(4.2)z=z endfor endifend 注: (1)算法时间复杂度 t(n)=2t(n/2)+O(n)t(n)=O(nlogn)5n=8 的 FFT 蝶式计算图:图 3.1 FFT 蝶式计算图n=6 的 FFT 递归计算流程图:图 3.2 FFT 递归计算流程图3.2流程图6飞是 否 是否是否输入
9、序列长度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= ？Xj+k Xj+klXj+k+l*W(size_x/2/l)*k Xj+k+l-1K 加 17四 源程序代码及运行结果4.1源程序代码/*FFT*/#include /整个程序输入和输出利用同一个空间 xN，节约空间#include #include #define N 1000 /定
10、义输入或者输出空间的最大长度typedefstructdouble real;doubleimg;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
11、 divi(complex,complex,complex *); /*复数除法*/void output(); /*输出结果*/complex xN,*W; /*输出序列的值*/intsize_x=0; /*输入序列的长度，只限 2 的 N 次方*/double PI; /pi 的值int main()inti;system(“cls“);PI=atan(1)*4;8printf(“Please input the size of x:n“); /*输入序列的长度*/scanf(“%d“,printf(“Please input the data in xN:(such as:5 6)n“)
12、;for(i=0;isize_x;i+) /*输入序列对应的值*/scanf(“%lf %lf“,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 个） 。*/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)for(k=0;kl;k+) /算出第 i 级内 j 组蝶形单元的结果 /算出 j 组中第 k 个蝶形单元mul(xj+k+l,W(size_x/2/l)*k, /*size/2/l 是该级 W 的相邻上标差，l 是该级该组取的 W 总个数*/add(xj+k,product,sub(xj+k,product,
