快速傅里葉算法

**快速傅立葉變換**（英語：**Fast Fourier Transform, FFT**），是快速計算序列的離散傅立葉變換（DFT）或其反轉換的方法。傅立葉分析將訊號從原始域（通常是時間或空間）轉換到頻域的表示或者逆過來轉換。FFT 會通過把 DFT 矩陣分解為稀疏（大多為零）因子之積來快速計算此類變換。 因此，它能夠將計算 DFT 的複雜度從只用 DFT 定義計算需要的 $O(n^{2})$ ，降低到 $O(n\log n)$ ，其中 $n$ 為資料大小。

快速傅立葉變換廣泛的應用於工程、科學和數學領域。這裡的基本思想在 1965 年才得到普及，但早在 1805 年就已推導出來。 1994 年美國數學家吉爾伯特・斯特朗把 FFT 描述為「我們一生中最重要的數值演算法」，它還被 IEEE 科學與工程計算期刊列入 20 世紀十大演算法。

## 定義和速度

用 FFT 計算 DFT 會得到與直接用 DFT 定義計算相同的結果；最重要的區別是 FFT 更快。（由於捨入誤差的存在，許多 FFT 演算法還會比直接運用定義求值精確很多，後面會討論到這一點。）

令 _x_0, ...., _x_N-1 為複數。DFT 由下式定義

　
: $X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{i2\pi k{\frac {n}{N}}}}\qquad k=0,\dots ,N-1.$

直接按這個定義求值需要 O (_N_2) 次運算：_X_k 共有 _N_ 個輸出，每個輸出需要 _N_ 項求和。直接使用 DFT 運算需使用 N 個複數乘法 (4N 個實數乘法) 與 N-1 個複數加法 (2N-2 個實數加法)，因此，計算使用 DFT 所有 N 點的值需要 _N_2 複數乘法與 _N_2-N 個複數加法。FFT 則是能夠在 O (_N_ log _N_) 次操作計算出相同結果的任何方法。更準確的說，所有已知的 FFT 演算法都需要 O (_N_ log _N_) 次運算（技術上 O 只標記上界），雖然還沒有已知的證據證明更低的複雜度是不可能的。

要說明 FFT 節省時間的方式，就得考慮複數相乘和相加的次數。直接計算 DFT 的值涉及到 _N_2 次複數相乘和 _N_(_N_−1) 次複數相加（可以通過削去瑣碎運算（如乘以 1）來節省 O (_N_) 次運算）。眾所周知的基 2 庫利 - 圖基演算法，_N_ 為 2 的冪，可以只用 (_N_/2)log2(_N_) 次複數乘法（再次忽略乘以 1 的簡化）和 _N_log2(_N_) 次加法就可以得到相同結果。在實際中，現代電腦通常的實際效能通常不受算術運算的速度和對複雜主體的分析主導，但是從 O (_N_2) 到 O (_N_ log _N_) 的母體改進仍然能夠體現出來。

## 一般的簡化理論

假設一個 M\*N 型矩陣 **S** 可分解成列向量以及行向量相乘：

$$
\mathbf {S} ={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{m}\end{bmatrix}}{\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}
$$

若 ${\begin{bmatrix}a_{1}&a_{2}&\cdots &a_{m}\end{bmatrix}}^{T}$ 有 $M_{0}$ 個相異的非平凡值（ $a_{m}\neq \pm 2^{k},a_{m}\neq \pm 2^{k}a_{n}$ where $m\neq n$ ）　

 ${\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}$ 有 $N_{0}$ 個相異的非平凡值

則 **S** 共需要 $M_{0}*N_{0}$ 個乘法。

$$
{\begin{bmatrix}Z[1]\\Z[2]\\\vdots \\Z[N]\end{bmatrix}}=\mathbf {S} {\begin{bmatrix}X[1]\\X[2]\\\vdots \\X[N]\end{bmatrix}}={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{m}\end{bmatrix}}{\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}{\begin{bmatrix}X[1]\\X[2]\\\vdots \\X[N]\end{bmatrix}}
$$

Step 1： $Z_{a}=b_{1}X[1]+b_{2}X[2]+\cdots +b_{n}X[N]$

Step 2： $Z[1]=a_{1}Z_{a},Z[2]=a_{2}Z_{a},\cdots ,Z[N]=a_{m}Z_{a}$

簡化理論的變型：

$$
\mathbf {S} ={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{m}\end{bmatrix}}{\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}+\mathbf {S} _{1}
$$

 $S_{1}$ 也是一個 M\*N 的矩陣。

若 $S_{1}$ 有 $P_{1}$ 個值不等於 0，則 $\mathbf {S} $ 的乘法量上限為 $M_{0}+N_{0}+P_{1}$ 。

## 快速傅立葉變換乘法量的計算

假設 $N=P_{1}\times P_{2}\times \cdots \times P_{k}$ ，其中 $P_{1},P_{2},\cdots ,P_{k}$ 彼此互質

 $\mathbf {P_{k}} $ 點 DFT 的乘法量為 $\mathbf {B_{k}} $ ，則 $\mathbf {N} $ 點 DFT 的乘法量為：

　
: ${\frac {N}{P_{1}}}B_{1}+{\frac {N}{P_{2}}}B_{2}+\cdots \cdots +{\frac {N}{P_{k}}}B_{k}$

假設 $\mathbf {N=P^{c}} $ ，P 是一個質數。

若 $\mathbf {N_{1}=P^{a}} $ 點的 DFT 需要的乘法量為 $\mathbf {B_{1}} $

且 $n_{1}\times n_{2}$ 當中（ $n_{1}=0,1,\cdots ,N_{1}-1,\quad n_{2}=0,1,\cdots ,N_{2}-1$ ）

有 $D_{1}$ 個值不為 ${\frac {N}{12}}$ 及 ${\frac {N}{8}}$ 的倍數，

有 $D_{2}$ 個值為 ${\frac {N}{12}}$ 及 ${\frac {N}{8}}$ 的倍數，但不為 ${\frac {N}{4}}$ 的倍數，

則 N 點 DFT 的乘法量為：

　
: $\mathbf {N_{2}B_{1}+N_{1}B_{2}+3D_{1}+2D_{2}} $

## 庫利 - 圖基演算法

庫利 - 圖基演算法是最常見的 FFT 演算法。這一方法以分治法為策略遞迴地將長度為 $N=N_{1}N_{2}$ 的離散傅立葉變換分解為長度為 $N_{1}$ 的 $N_{2}$ 個較短序列的離散傅立葉變換，以及與 $\mathrm {O} (N)$ 個旋轉因子的複數乘法。

這種方法以及 FFT 的基本思路在 1965 年 J. W. Cooley 和 J. W. Tukey 合作發表 _An algorithm for the machine calculation of complex Fourier series_ 之後開始為人所知。但後來發現，實際上這兩位作者只是重新發明了高斯在 1805 年就已經提出的演算法（此演算法在歷史上數次以各種形式被再次提出）。

庫利 - 圖基演算法最有名的應用，是將序列長為 _N_ 的 DFT 分割為兩個長為 _N/2_ 的子序列的 DFT，因此這一應用只適用於序列長度為 2 的冪的 DFT 計算，即基 2-FFT。實際上，如同高斯和 Cooley 與 Tukey 都指出的那樣，Cooley-Tukey 演算法也可以用於序列長度 _N_ 為任意因數分解形式的 DFT，即混合基 FFT，而且還可以應用於其他諸如分裂基 FFT 等變種。儘管 Cooley-Tukey 演算法的基本思路是採用遞迴的方法進行計算，大多數傳統的演算法實現都將顯式的遞迴演算法覆寫為非遞迴的形式。另外，因為 Cooley-Tukey 演算法是將 DFT 分解為較小長度的多個 DFT，因此它可以同任一種其他的 DFT 演算法聯合使用。

### 設計思想

下面，我們用 **N 次單位根** $W_{N}$ 來表示 $e^{-j{\frac {2\pi }{N}}}$ 。

 $W_{N}$ 的性質：

1. **週期性**， $W_{N}$ 具有週期 $N$ ，即 $W_{N}^{k+N}=W_{N}^{k}$
2. **對稱性**： $W_{N}^{k+{\frac {N}{2}}}=-W_{N}^{k}$ 。
3. 若 $m$ 是 $N$ 的因數， $W_{N}^{mkn}=W_{\frac {N}{m}}^{kn}$

為了簡單起見，我們下面設待變換序列長度 $n=2^{r}$ 。根據上面單位根的對稱性，求級數 $y_{k}=\sum _{n=0}^{N-1}W_{N}^{kn}x_{n}$ 時，可以將求和區間分為兩部分：

$$
{\begin{matrix}y_{k}=\sum _{n=2t}W_{N}^{kn}x_{n}+\sum _{n=2t+1}W_{N}^{kn}x_{n}\\=\sum _{t}W_{\frac {N}{2}}^{kt}x_{2t}+W_{N}^{k}\sum _{t}W_{\frac {N}{2}}^{kt}x_{2t+1}\\=F_{even}(k)+W_{N}^{k}F_{odd}(k)&&&&&&(i\in \mathbb {Z} )\end{matrix}}
$$

 $F_{odd}(k)$ 和 $F_{even}(k)$ 是兩個分別關於序列 $\left\{x_{n}\right\}_{0}^{N-1}$ 奇數號和偶數號序列 N/2 點變換。由此式只能計算出 $y_{k}$ 的前 N/2 個點，對於後 N/2 個點，注意 $F_{odd}(k)$ 和 $F_{even}(k)$ 都是週期為 N/2 的函式，由單位根的對稱性，於是有以下變換公式：

 * $y_{k+{\frac {N}{2}}}=F_{even}(k)-W_{N}^{k}F_{odd}(k)$
* $y_{k}=F_{even}(k)+W_{N}^{k}F_{odd}(k)$ 。

這樣，一個 N 點變換就分解成了兩個 N/2 點變換。照這樣可繼續分解下去。這就是**庫利 - 圖基快速傅立葉變換**演算法的基本原理。根據主定理不難分析出此時演算法的時間複雜度為 $\mathrm {O} (N\log N)$

### 演算法實現

* 蝶形結網路和位反轉（Bit Reversal）：

　
: 首先將 $n=2^{N}$ 個輸入點列按二進位進行編號，然後對各個編號按位元倒置並按此重新排序。例如，對於一個 8 點變換，

　
: 001 倒置以後變成 100

　
: 000 → 000

　
: 001 → 100

　
: 010 → 010

　
: 011 → 110

　
: 100 → 001

　
: 101 → 101

　
: 110 → 011

　
: 111 → 111

　
: 倒置後的編號為 {0,4,2,6,1,5,3,7}。

　
: 然後將這 n 個點列作為輸入傳送到蝶形結網路中，注意將因子 $W_{N}^{k}$ 逐層加入到蝶形網路中。

### 演算法複雜度

由於按蝶形結網路計算 n 點變換要進行 log _n_ 層計算，每層計算 n 個點的變換，故演算法的時間複雜度為 $\mathrm {O} (n\log n)$ 。

## 其他演算法

在 Cooley-Tukey 演算法之外還有其他 DFT 的快速演算法。對於長度 $N=N_{1}N_{2}$ 且 $N_{1}$ 與 $N_{2}$ 互質的序列，可以採用基於中國剩餘定理的互質因子演算法將 N 長序列的 DFT 分解為兩個子序列的 DFT。與 Cooley-Tukey 演算法不同的是，互質因子演算法不需要旋轉因子。

Rader-Brenner 演算法是類似於 Cooley-Tukey 演算法，但是採用的旋轉因子都是純虛數，以增加加法運算和降低了數值穩定性為代價減少了乘法運算。這方法之後被 split-radix variant of Cooley-Tukey 所取代，與 Rader-Brenner 演算法相比，有一樣多的乘法量，卻有較少的加法量，且不犧牲數值的準確性。

Bruun 以及 QFT 演算法是不斷的把 DFT 分成許多較小的 DFT 運算。（Rader-Brenner 以及 QFT 演算法是為了 2 的指數所設計的演算法，但依然可以適用在可分解的整數上。Bruun 演算法則可以運用在可被分成偶數個運算的數位）。尤其是 Bruun 演算法，把 FFT 看成是 $z^{N}-1$ ，並把它分解成 $z^{M-1}$ 與 $z^{2M}+az^{M}+1$ 的形式。

另一個從多項式觀點的快速傅立葉變換法是 Winograd 演算法。此演算法把 $z^{N}-1$ 分解成 cyclotomic 多項式，而這些多項式的係數通常為 1，0，-1。這樣只需要很少的乘法量（如果有需要的話），所以 winograd 是可以得到最少乘法量的快速傅立葉演算法，對於較小的數位，可以找出有效率的算方式。更精確地說，winograd 演算法讓 DFT 可以用 $2^{k}$ 點的 DFT 來簡化，但減少乘法量的同時，也增加了非常多的加法量。Winograd 也可以利用剩餘值定理來簡化 DFT。

Rader 演算法提出了利用點數為 N（N 為質數）的 DFT 進行長度為 N-1 的迴旋摺積來表示原本的 DFT，如此就可利用摺積用一對基本的 FFT 來計算 DFT。另一個 prime-size 的 FFT 演算法為 chirp-Z 演算法。此法也是將 DFT 用摺積來表示，此法與 Rader 演算法相比，能運用在更一般的轉換上，其轉換的基礎為 Z 轉換（Rabiner et al., 1969）。

## 實數或對稱資料專用的演算法

在許多的運用當中，要進行 DFT 的資料是純實數，如此一來經過 DFT 的結果會滿足對稱性：

　
: $\mathbf {X} _{N-k}$ = $\mathbf {X} _{k}^{*}$

而有一些演算法是專門為這種情形設計的（e.g. Sorensen, 1987）。另一些則是利用舊有的演算法（e.g. Cooley-Tukey），刪去一些不必要的演算步驟，如此省下了記憶體的使用，也省下了時間。另一方面，也可以把一個偶數長度且純實數的 DFT，用長度為原本一半的複數型態 DFT 來表示（實數項為原本純實數資料的偶數項，虛數項則為奇數項）。

**在 Matlab 上用一次 N 點 FFT 計算兩個 N 點實數訊號 x,y 的 FFT:**

一度人們認為，用離散哈特利轉換（Discrete Hartley Transform）來處理純實數的 DFT 會更有效率，但接著人們發現，對於同樣點數的純實數 DFT，經過設計的 FFT，可以比 DHT 省下更多的運算。Bruun 演算法是第一個試著從減少實數輸入的 DFT 運算量的演算法，但此法並沒有成為人們普遍使用的方法。

對於實數輸入，且輸入為偶對稱或奇對稱的情形，可以更進一步的省下時間以及記憶體，此時 DFT 可以用離散餘弦轉換或離散正弦轉換來代替（Discrete cosine/sine transforms）。由於 DCT/DST 也可以設計出 FFT 的演算法，故在此種情形下，此方法取代了對 DFT 設計的 FFT 演算法。

DFT 可以應用在頻譜分析以及做摺積的運算，而在此處，不同應用可以用不同的演算法來取代，列表如下：

用來做頻譜分析的情況下，DFT 可用下列的演算法代替：

* DCT
* DST
* DHT
* 正交基底的擴充（orthogonal basis expantion）包括正交多項式（orthogonal polynomials）以及 CDMA.
* Walsh（Hadamard）轉換
* Haar 轉換
* 小波（wavelet）轉換
* 時頻分布（time-frequency distribution）

用來做摺積的情況下，DFT 可用下列的演算法代替：

* DCT
* DST
* DHT
* 直接做摺積（direct computing）
* 分段式 DFT 摺積（sectioned DFT convolution）
* 威諾格拉德快速傅立葉變換演算法
* 沃爾什（Walsh、Hadamard）轉換
* 數論轉換

## 複雜度以及運算量的極限

長久以來，人們對於求出快速傅立葉變換的複雜度下限以及需要多少的運算量感到很有興趣，而實際上也還有許多問題需要解決。即使是用較簡單的情形，即 $2^{k}$ 點的 DFT，也還沒能夠嚴謹的證明出 FFT 至少需要 $\Omega (N\log N)$ （不比 $N\log N$ 小）的運算量，目前也沒有發現複雜度更低的演算法。通常數學運算量的多寡會是運算效率好壞最主要的因素，但在現實中，有許多因素也會有很大的影響，如快取記憶體以及 CPU 均有很大的影響。

在 1978 年，Winograd 率先導出一個較嚴謹的 FFT 所需乘法量的下限： $\Theta (N)$ 。當 $N=2^{k}$ 時，DFT 只需要 $4N-2\log _{2}^{2}N-2\log _{2}N-4$ 次無理實數的乘法即可以計算出來。更詳盡，且也能趨近此下限的演算法也一一被提出（Heideman & Burrus, 1986; Duhamel, 1990）。很可惜的是，這些演算法，都需要很大量的加法計算，目前的硬體無法克服這個問題。

對於所需加法量的數目，雖然我們可以在某些受限制的假設下，推得其下限，但目前並沒有一個精確的下限被推導出來。1973 年，Morgenstern 在乘法常數趨近巨大的情形下（對大部分的 FFT 演算法為真，但不是全部）推導出加法量的下限： $\Omega \left(N\log N\right)$ 。Pan（1986）在假設 FFT 演算法的不同步的情形有其極限下證明出加法量的下限 $\Omega (NlogN)$ ，但一般來說，此假設相當的不明確。長度為 $N=2^{k}$ 的情形下，在某些假設下，Papadimitriou（1979）提出使用 Cooley-Tukey 演算法所需的複數加法量 $N\log _{2}N$ 是最少的。到目前為止，在長度為 $N=2^{k}$ 情況，還沒有任何 FFT 的演算法可以讓複數的加法量比 $N\log _{2}N$ 還少。

還有一個問題是如何把乘法量與加法量的總和最小化，有時候稱作 "演算複雜度"（在這裡考慮的是實際的運算量，而不是漸近複雜度）。同樣的，沒有一個嚴謹下限被證明出來。從 1968 年開始， $N=2^{k}$ 點 DFT 而言，split-radix FFT 演算法需要最少的運算量，在 $N>1$ 的情形下，其需要 $4N\log _{2}N-6N+8$ 個乘法運算以及加法運算。最近有人導出更低的運算量： ${\frac {34}{9}}N\log _{2}N$ 。（Johnson and Frigo, 2007; Lundy and Van Buskirk, 2007）

大多數嘗試要降低或者證明 FFT 複雜度下限的人都把焦點放在複數資料輸入的情況，因其為最簡單的情形。但是，複數資料輸入的 FFT 演算法，與實數資料輸入的 FFT 演算法，離散餘弦轉換（DCT），離散哈特列轉換（DHT），以及其他的演算法，均有很大的關連性。故任何一個演算法，在複雜度上有任何的改善的話，其他的演算法複雜度也會馬上獲得改善（Duhamel & Vetterli, 1990）。

## 快速演算法查表

當輸入訊號長度為 N 時:

N = 1~60

N < 1000

N > 1000

## 參閱

* 離散傅立葉變換
* 並列快速傅立葉變換
* 快速數論變換

## 參考資料

## 延伸閱讀