利用C语言实现FFT、IFFT运算

一、FFT算法理论

上述分别为FFT、IFFT公式。下面首先讨论FFT的算法实现。 本文采用输入逆序、输出顺序的FFT计算方法。实质上就是在时域对x(n)进行“奇偶分类”、在频域上对X(k)进行“前后分类”。 值得说明的是，这里的“奇”和“偶”是相对的概念，并不完全是通常我们所理解的“奇”和“偶”。下面将给出一个例子进行说明： ·图中，总共可以分为 “三级” ，每一级包含若干 “蝶形计算单元”。每一个蝶形单元的结构如下： 因此，进行FFT运算只需考虑三件事情。一、将输入序列倒序排列，也就是进行奇偶分类。二、实现最基本的蝶形计算单元。三、根据规律，划分出FFT运算的级数，并且确定出每一级中参与蝶形运算的数据分组。

二、FFT运算的代码实现

输入倒序排列（奇偶分类）。事实上，将输入倒序排列就是进行了奇偶分类。如果FFT的点数为8，那么这8个数的二进制表示分别为000，001，010，011，100，101，110，111。进行二进制层面倒序排列就变为000，100，010，110，001，101，011，111；重新转变为十进制就是0，4，2，6，1，5，3，7。 上面的计算方式对人类来说是很方便的，但是对于电脑来说，效果并不是太好，因此需要另外发现一种规律。输入奇偶分类（倒序排列）。上面已经提到，倒序排列实际上就是在进行奇偶分类，而进行奇偶分类对电脑来说会简单很多。总结1中提到的规律，三位二进制数可以理解为FFT运算的级数。我们可以发现，在倒序排列完成后，“偶数”总是排在“奇数”前面。 这里需要提出“偶中的偶”，“奇中的奇”的概念。即每次进行“余2”运算，结果是0的都排在前面。现在看来，电脑实现算法的方式就是每次进行“余2”运算，进行奇偶分类；然后进行“除2”运算，再进行“余2”运算，进行奇偶分类……直到最后排序完成。 实现代码如下： void reverse() //输入倒序 { complex temp; char* flag = (char*)malloc(N * sizeof(char)); //用于标记是否被换过，若为1，则已经换过 int i, j, current_n, target_n; for (i = 0;i < N;i++) { flag[i] = 0; } for (i = 0;i < N;i++) { current_n = i; target_n = 0; //获取两个互为倒序的标号，换种思路 for (j = 0;j < bits;j++) { target_n = target_n + (int)((current_n % 2) * pow(2, bits - j - 1)); current_n /= 2; } current_n = i; //对应标号值互换 if (current_n != target_n && flag[current_n] != 1) { temp = x[current_n]; x[current_n] = x[target_n]; x[target_n] = temp; flag[target_n] = 1; } } free(flag); }

说明：complex为自己定义的复数数据结构。再次强调，奇偶分类就是倒序排序。二者是一样的。

蝶形运算单元的实现。分析蝶形运算单元，两输入、两输出，因此蝶形运算单元需要注意的就是输入、输出采用同一片内存存储即可，也就是让输出数据覆盖输入数据，这样可以有效节约存储资源。 void butterfly(int x1_point, int x2_point, complex wn) //蝶形计算单元 { complex result1, result2, T; T = multiplication(x[x2_point], wn); result1 = add(x[x1_point], T); result2 = sub(x[x1_point], T); x[x1_point] = result1; x[x2_point] = result2; }

说明：上面multiplication、add、sub方法均为自己定义，wn为旋转因子。

进行单级FFT运算。需要注意的就是总结规律。总结出两对偶结点的“结点跨距”、“分组间隔”、“分组数”……有了这些才能够编写将参与蝶形运算的数据分类配对。 void single_fft(int m) //进行单级的fft运算 { //首先进行分组 int point_distance = (int)pow(2, m - 1); //结点跨距 int group_distance = 2 * point_distance; //分组间隔 int group_count = N / group_distance; //分组数 int group_header = 0; int x1, x2; //参与计算两结点的标号 complex w; int i, j; //循环控制符 //分组 for (i = 0;i < group_count;i++) { //获取一组的标号范围 group_header = i * group_distance; //将每组分成两半，前一半均为x1,后一半均为x2 for (j = 0;j < group_distance / 2;j++) { w = w_builder(m, j); x1 = j + group_header; x2 = x1 + point_distance; butterfly(x1, x2, w); } } } 有了单级FFT运算以后，将每一级组合就可以实现FFT运算了。 complex* fft(int fft_point, int xn_length, complex* xn) //fft最终封装 { int i; //初始化 fft_init(fft_point, xn_length, xn); //输入倒序 reverse(); //fft运算 for (i = 1;i <= bits;i++) { single_fft(i); } return x; }

三、IFFT算法理论

从上面的等式可以看出，括号里面的内容实际上就是FFT运算。也就是说只需要对输入、输出序列进行一定的处理，就可以利用刚刚已经做好的FFT算法实现IFFT的运算！

四、IFFT运算的代码实现

complex* ifft(int xk_length, complex* xk) //ifft最终封装 { int i; //对输入序列进行处理 ifft_init(xk_length, xk); //套用fft算法 reverse(); for (i = 1;i <= bits;i++) { single_fft(i); } //对输出序列进行处理 for (i = 0;i < N;i++) { x[i] = conjugate(x[i]); //求共轭 x[i].real /= N; x[i].imaginary /= N; } return x; }

五、总结

本文也程序也考虑到了FFT自动补零问题的处理，大家可以参考完整代码（FFT、IFFT均验证正确）。FFT.zip