作者：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 ) 為輸出的蝶形運算的另一側，這兩個輸出分別存放在正方形上下兩邊對應的行向量進行矩陣運算結果位置。
我的一個新思路是以矩陣對角線元素為主導變量形成的算法實現可能更為簡單，但目前還沒嘗試。