作者：Abraham 转载请标明出处，谢谢！
找到一个通俗易懂并且神奇并且有趣的FFT算法C语言实现教程：http://www.katjaas.nl/FFTimplement/FFTimplement.html
只要对矩阵比较熟悉就能在教程的辅助下很快实现FFT算法的C代码。
这个教程的第二部分 “Bitwiser FFT” 是一个基于位运算的FFT优化代码，我花了一段时间反复看才理解，不过理解之后印象深刻。难理解的地方主要是作者所说的 lo - set 和 odd - even 两对概念。作者开始介绍的时候这似乎是两个不同的概念，后来在程序中又变成了统一的 odd - even 命名，然而其实质是两个不同的概念。一个是为了跳过矩阵中的子空矩阵，通过与 span 求“或”实现（ odd |= span ），另一个是为了求蝶形运算对中另一个元素的下标（横坐标），通过与 span 求“异或”实现（ even = odd ^ span ）。还需要注意的是，在每个 stage 内的每个 loop 都是以 odd 变量（这个名字太容易引起歧义，我索性不去管它的字面意思，而把它理解为蝶形运算对里靠后的那个元素，或者把它理解为二进制数 span 中为 "1" 的那一位的奇偶）为主导变量的，即 odd 变导致 even 变然后进行其他运算。
图形化理解 “Bitwiser FFT” 一节的方法在于，把每个 stage 内的每个 loop 理解为一个正方形的四个顶点的运算（如 “Bitwiser FFT” 一节的图2、4、6），正方形的边长是 span + 1 。这个算法中的主导变量 odd 即代表正方形右上角顶点元素的纵坐标（其纵坐标指向蝶形运算输入的元素下标，其横坐标指向蝶形运算输出的元素下标）。整个 stage 内的过程可以描述为三步：
1. odd 自增（包括跳过空白子矩阵）。
2. 计算 even 。
3. 蝶形运算产生新的 odd 和 even 。
这个算法是莱昂斯的《数字信号处理》中提到的按位反转输出的频域抽取FFT（见中文版第96页图4.12）。其蝶形运算公式为：x'' = x + y 和 y'' = WkN * x - WkN * y = WkN * ( x - y ) 。图4.12和 “Bitwiser FFT” 中相同的一点就是在这种算法内每次蝶形运算的输入和输出总是平行的，即当蝶形运算的输入是第 0 个和第 4 个元素时，其输出也是第 0 个和第 4 个元素。这个蝶形运算就对应上一段提到的正方形。正方形的上边（左上、右上两个顶点）代表以 x 和 y 为输入， 以x + y 为输出的蝶形运算一侧，下边（左下、右下）两个顶点代表以 x 和 y 为输入， 以 WkN * ( x - y ) 为输出的蝶形运算的另一侧，这两个输出分别存放在正方形上下两边对应的行向量进行矩阵运算结果位置。
我的一个新思路是以矩阵对角线元素为主导变量形成的算法实现可能更为简单，但目前还没尝试。