快速傅立叶变换（FFT）的C++实现收藏标准的离散傅立叶DFT 变换形式如：y k=Σj=0n-1a jωn-kj = A (ωn-k).（ωn k为复数1 的第k 个n 次方根，且定义多项式A (x)=Σj=0n-1a j x j）而离散傅立叶逆变换IDFT (Inverse DFT)形式如：a j=(Σk=0n-1y kωn kj)/n .yk=Σj=0n-1 ajωn-kj = A (ωn-k).（ωnk 为复数1 的第k 个n 次方根，且定义多项式 A (x) = Σj=0n-1 ajxj ）而离散傅立叶逆变换IDFT (Inverse DFT)形式如：aj=(Σk=0n-1 ykωnkj)/n .以下不同颜色内容为引用并加以修正：快速傅立叶变换（Fast Fourier Transform，FFT）是离散傅立叶变换（Discrete Fourier transform，DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

设Xn 为N 项的复数序列，由DFT 变换，任一Xi 的计算都需要N 次复数乘法和N -1 次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的Xi ，即N 点DFT 变换大约就需要N2 次运算。

当N =1024 点甚至更多的时候，需要N2 = 1048576 次运算，在FFT 中，利用ωn 的周期性和对称性，把一个N 项序列（设N 为偶数），分为两个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% 的运算量。

而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT 运算单元，那么N 点的DFT 变换就只需要N * log2N 次的运算，N = 1024 点时，运算量仅有10240 次，是先前的直接算法的1% ，点数越多，运算量的节约就越大，这就是FFT 的优越性。

FFT 的实现可以自顶而下，采用递归，但是对于硬件实现成本高，对于软件实现都不够高效，改用迭代较好，自底而上地解决问题。

感觉和归并排序的迭代版很类似，不过先要采用“位反转置换”的方法把Xi 放到合适的位置，设i 和j 互为s = log2N 位二进制的回文数，假设s = 3, i = (110)2 = 6, j = (011)2 = 3, 如果i ≠j , 那么Xi 和Xj 应该互换位置。

（关于这个回文数的生成，是很有趣而且是很基本的操作，想当初偶初学C++ 的时候就有这样的习题。

）当“位反转置换”完成后，先将每一个Xi 看作是独立的多项式，然后两个两个地将它们合并成一个多项式（每个多项式有2 项），合并实际上是“蝶形运算”（Butterfly Operation, 参考《算法导论》吧^_^），继续合并（第二次的每个多项式有4 项），直到只剩下一个多项式（有N 项），这样，合并的层数就是log2N ，每层都有N 次操作，所以总共有N * log2N 次操作。

迭代过程如下图所示，自底而上。

C++ 源程序，如下：//// 快速傅立叶变换Fast Fourier Transform//********************.cn// 2007-07-20// 版本2.0// 改进了《算法导论》的算法，旋转因子取ωn-kj (ωnkj 的共轭复数)// 且只计算n / 2 次，而未改进前需要计算(n * lg n) / 2 次。

//#include <stdio.h>#include <stdlib.h>#include <math.h>const int N = 1024;const float PI = 3.1416;inline void swap (float &a, float &b){float t;t = a;a = b;b = t;}void bitrp (float xreal [], float ximag [], int n){// 位反转置换Bit-reversal Permutationint i, j, a, b, p;for (i = 1, p = 0; i < n; i *= 2){p ++;}for (i = 0; i < n; i ++){a = i;b = 0;for (j = 0; j < p; j ++){b = (b << 1) + (a & 1); // b = b * 2 + a % 2;a >>= 1; // a = a / 2;}if ( b > i){swap (xreal [i], xreal [b]);swap (ximag [i], ximag [b]);}}}void FFT(float xreal [], float ximag [], int n){// 快速傅立叶变换，将复数x 变换后仍保存在x 中，xreal, ximag 分别是x 的实部和虚部float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;int m, k, j, t, index1, index2;bitrp (xreal, ximag, n);// 计算1 的前n / 2 个n 次方根的共轭复数W'j = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1arg = - 2 * PI / n;treal = cos (arg);timag = sin (arg);wreal [0] = 1.0;wimag [0] = 0.0;for (j = 1; j < n / 2; j ++){wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;}for (m = 2; m <= n; m *= 2){for (k = 0; k < n; k += m){for (j = 0; j < m / 2; j ++){index1 = k + j;index2 = index1 + m / 2;t = n * j / m; // 旋转因子w 的实部在wreal [] 中的下标为ttreal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];ureal = xreal [index1];uimag = ximag [index1];xreal [index1] = ureal + treal;ximag [index1] = uimag + timag;xreal [index2] = ureal - treal;ximag [index2] = uimag - timag;}}}}void IFFT (float xreal [], float ximag [], int n){// 快速傅立叶逆变换float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;int m, k, j, t, index1, index2;bitrp (xreal, ximag, n);// 计算1 的前n / 2 个n 次方根Wj = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 arg = 2 * PI / n;treal = cos (arg);timag = sin (arg);wreal [0] = 1.0;wimag [0] = 0.0;for (j = 1; j < n / 2; j ++){wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;}for (m = 2; m <= n; m *= 2){for (k = 0; k < n; k += m){for (j = 0; j < m / 2; j ++){index1 = k + j;index2 = index1 + m / 2;t = n * j / m; // 旋转因子w 的实部在wreal [] 中的下标为ttreal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];ureal = xreal [index1];uimag = ximag [index1];xreal [index1] = ureal + treal;ximag [index1] = uimag + timag;xreal [index2] = ureal - treal;ximag [index2] = uimag - timag;}}}for (j=0; j < n; j ++){xreal [j] /= n;ximag [j] /= n;}}void FFT_test (){char inputfile [] = "input.txt"; // 从文件input.txt 中读入原始数据char outputfile [] = "output.txt"; // 将结果输出到文件output.txt 中float xreal [N] = {}, ximag [N] = {};int n, i;FILE *input, *output;if (!(input = fopen (inputfile, "r"))){printf ("Cannot open file. ");exit (1);}if (!(output = fopen (outputfile, "w"))){printf ("Cannot open file. ");exit (1);}i = 0;while ((fscanf (input, "%f%f", xreal + i, ximag + i)) != EOF){i ++;}n = i; // 要求n 为2 的整数幂while (i > 1){if (i % 2){fprintf (output, "%d is not a power of 2! ", n);exit (1);}i /= 2;}FFT (xreal, ximag, n);fprintf (output, "FFT: i real imag ");for (i = 0; i < n; i ++){fprintf (output, "%4d %8.4f %8.4f ", i, xreal [i], ximag [i]); }fprintf (output, "================================= ");IFFT (xreal, ximag, n);fprintf (output, "IFFT: i real imag ");for (i = 0; i < n; i ++){fprintf (output, "%4d %8.4f %8.4f ", i, xreal [i], ximag [i]); }if ( fclose (input)){printf ("File close error. ");exit (1);}if ( fclose (output)){printf ("File close error. ");exit (1);}}int main (){FFT_test ();return 0;}//////////////////////////////////////////////// // sample: input.txt////////////////////////////////////////////////0.5751 00.4514 00.0439 00.0272 00.3127 00.0129 00.3840 00.6831 00.0928 00.0353 00.6124 00.6085 00.0158 00.0164 00.1901 00.5869 0//////////////////////////////////////////////// // sample: output.txt//////////////////////////////////////////////// FFT:i real imag0 4.6485 0.00001 0.0176 0.31222 1.1114 0.04293 1.6776 -0.13534 -0.2340 1.38975 0.3652 -1.25896 -0.4325 0.20737 -0.1312 0.37638 -0.1949 0.00009 -0.1312 -0.376310 -0.4326 -0.207311 0.3652 1.258912 -0.2340 -1.389713 1.6776 0.135314 1.1113 -0.042915 0.0176 -0.3122================================= IFFT:i real imag0 0.5751 -0.00001 0.4514 0.00002 0.0439 -0.00003 0.0272 -0.00004 0.3127 -0.00005 0.0129 -0.00006 0.3840 -0.00007 0.6831 0.00008 0.0928 0.00009 0.0353 -0.000010 0.6124 0.000011 0.6085 0.000012 0.0158 0.000013 0.0164 0.000014 0.1901 0.000015 0.5869 -0.0000。