---
tags: signal
---
# 快速傅立葉轉換 (FFT)
[TOC]
## 動機:多項式乘法
### Coefficient Representation
給定一個$n-1$階多項式(或稱***degree bound***為$n$)$A(x) = \sum_{j = 0}^{n - 1} a_j x^j$,我們可以用$n$維向量$(a_0, a_1, \ldots , a_{n-1})$來唯一表示這個多項式,其中分量$a_j$代表的是$x^j$對應的係數,此方式稱之為***coefficient representation***。在這個表示之下,我們可以輕易的在$O(n)$的時間內將兩個多項式$A(x)$(其coefficient representation為$a = (a_0, a_1, \ldots, a_{n-1})$)和多項式$B(x)$(其coefficient representation為$b= (b_0, b_1, \ldots, b_{n-1})$)相加所得的新多項式求出來,其coefficient representation即為$c = (c_0 , c_1, \ldots, c_{n-1})$。
然而,當我們想要得到兩個多項式相乘的結果之coefficient representation時(次方為$2n$的多項式),就必須計算$c = a \otimes b$,用最單純的定義計算必須花費$O(n^2)$的計算成本。
### Point-value Representation
另一種表示$n-1$階多項式的方式如下:我們可以取不同的$n$個點$x_0, \ldots, x_{n-1}$,並且計算該多項式在這$n$個點上的取值$y_0, \ldots, y_{n-1}$,如此一來我們可以用$n$個pair來表示該多項式$A(x)$:
$$
\begin{align}
\{(x_0, y_0), (x_1, y_1), \ldots, (x_n, y_{n-1})\} \\
\end{align}
$$
其中
$$
\begin{align}
y_k = A(x_k)
\end{align}
$$
當給定一個多項式$A(x)$的coefficient representation,$y_k$可以用***Horner's rule***在$O(n)$的複雜度下計算出來:
$$
A(x) = a_0 + x (a_1 + x(a_2 + \cdots + x(a_{n-1})) \cdots)
$$
而由於總共有$n$個函數值要計算,在這種直觀的計算方法下,由coefficient representation轉換為point-value representation的時間複雜度為$O(n^2)$。
相反地,當我們想由point-value representation推回coefficient representation時,此過程稱為插值(interpolation),以下的定理證明了插值是一個一對一的映射,也就是說,給定特定的$n$階多項式,其point-value representation和coefficient representation皆為唯一:
#### Theorem 1
For any set $\{(x_0, y_0), \ldots (x_{n-1}, y_{n-1})\}$ of $n$ point-value pairs such that all the $x_k$ values are distinct, there is a **unique** polynomial $A(x)$ of degree-bound $n$ such that $y_k = A(x_k)$ for $k = 0, 1, \ldots, n-1$.
**proof**
From the definition of point-value representation and coefficient representation, we have
$$
\begin{align}
\begin{pmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\
1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1}
\end{pmatrix}
\begin{pmatrix}
a_0 \\ a_1 \\ \vdots \\ a_{n-1}
\end{pmatrix}
= \begin{pmatrix}
y_0 \\ y_1 \\ \vdots \\ y_{n-1}
\end{pmatrix}
\end{align}
$$
The matrix on the left is known as a [Vander-monde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) and has determinant $\prod_{0 \le j < k \le n-1} (x_k - x_j)$. Thus, it's invertible if $x_k$ are distinct.
**remark**
如果由上面的矩陣乘法來由point-value representation求得coefficient representation,則透過LU分解演算法奇複雜度為$O(n^3)$,而透過[Lagranges's formula](https://en.wikipedia.org/wiki/Lagrange_polynomial)的複雜度則縮減為$O(n^2)$。
#### Point-Value Representation的優勢
如果給定兩個$n-1$階多項式$A(x)$和$B(x)$的point-value representation 如果我們想要知道$C(x) = A(x) B(x)$的point-value representation,那其實就是兩個函數在同一個點上的函數值相乘的結果,而由於$C(x)$會是一個degree bound為$2n$的多項式,也就是說我們可以先將$A(x)$和$B(x)$在$2n - 1$個點上取值,所以我們會得到兩個$2n-1$維的point-value representation:
$$
\{ (x_0, y_0), \ldots , (x_{2n-1}, y_{2n - 1}) \}
$$
和
$$
\{ (x_0, y'_0), \ldots , (x_{2n-1}, y'_{2n - 1}) \}
$$
則$C(x)$的point-value representation:
$$
\{ (x_0, y_0 y'_0), \ldots , (x_{2n-1}, y_{2n - 1} y'_{2n-1}) \}
$$
此方式的複雜度為$O(n)$,比起coefficient representation的$O(n^2)$還要有效率了許多。
從上方的發現在coefficient representation中的多項式加法和point-value representation中的多項式乘法的複雜度都只有$O(n)$。因此,我們只要找到兩種表示法之間有效率的轉換演算法,就能夠較有效率的計算上述兩種運算。而此想法的關鍵核心就在於$x_k$的選取,而其實,當我們將$x_k$選為$\omega_k = e^{-2\pi ki/n}$的時候,此轉換問題就會等價於離散傅立葉轉換的計算,也因此在接下來的小節中,我們探討的是如何將離散傅立葉轉換的時間複雜度下降至$O(n \log n)$。
## 一維快速傅立葉轉換
### 相關定義及定理
#### 離散傅立葉轉換
對於週期為$N$的離散訊號$x[j]$,其傅立葉級數的係數求法為
$$
\begin{align}
X[k] = \frac{1}{n}\sum_{j = 0}^{n - 1}
x[j] e^{- 2\pi k j/n}
\end{align}
$$
其中$i = \sqrt{-1}$,此亦稱為$x[i]$的離散傅立葉轉換(Discrete Fourier Transform, DFT),而在上一節中我們想要計算多項式
$$
A(x) = \sum_{j = 0}^{n-1} a_j x^j
$$
在點$x_k$的取值
$$
y_k = A(x_k) = \sum_{j = 0}^{n-1} a_j x_k^j
$$
當我們將$x_k$設為$e^{- 2\pi k i/n}$時,上面的式子就變成了傅立葉轉換的樣子,其中$a_j = x[j]$。另一種選擇$x_k$的方式是將取樣的點$x_k$設為$\omega^{2\pi ki / n} := \omega_k$(與之對應的轉換稱為離散傅立葉逆轉換,IDFT),本文接下來會會採取後者,不過在接下來的幾個定理和性質中,無論是上面哪種取法皆會成立,也因此兩種方式都能套用所介紹快速傅立葉轉換背後的演算法概念。
#### Proposition 1
Let $\omega_n = e^{2\pi i/n}$, then we have
$$
\omega_n^{nk} = 1, \forall k \in \mathbb Z_n
$$
That is, $\omega^k_n$ is the $n$th roots of unity. Moreover, we call $w_n$ the ***principal $\mathbf{n}$th root of unity***.
**proof**
Using Euler's formula $e^{iu} = \cos(u) + i \sin (u)$ completes the proof.
#### Lemma 1
$\forall n, k, d\in \mathbb Z,$
$$
w_{dn}^{dk} = w_n^k
$$
This implies that $w_k^2 = w_{n/2}$.
**proof**
$$
\begin{align}
w_{dn}^{dk} &= (e^{2 \pi i / dn})^{dk} \\
&= (e^{2 \pi i / dn})^k \\ &= w_n^k
\end{align}
$$
#### Lemma 2
$\forall n = 2k, k \in \mathbb Z$,
$$
w_n^{n/2} = w_2 = -1
$$
**proof**
Proved by definition and Euler formula.
#### Lemma 3 (Halving lemma)
If $n > 0$ is even, then the squares of the $n$ complex $n$th roots of unity are the $n/2$ complex $(n/2)$th roots of unity.
**proof**
By lemma 1, we have $(w_n^k)^2 = w_{n /2 }^k$ for $k \in \mathbb Z$.
這裡有一個直得注意的點是$(w_n^{k + n/2})^2 = (w_n^k)^2$。此引理是快速傅立葉轉換的一個關鍵性質。
#### Lemma 4 (Summation lemma)
$\forall n \ge 1$ and nonzero integer $k$ not divisible by $n$,
$$
\sum_{j = 0}^{n -1} (w_n^k)^j = 0
$$
**proof**
$$
\begin{align}
\sum_{j = 0}^{n - 1} (w_n^k)^j = \frac{w_n^{kn} - 1}{w_n^k - 1} = \frac{1 - 1}{w_n^k - 1} = 0
\end{align}
$$
### 快速傅立葉轉換(FFT)演算法
#### 原理
如同我們在本節一開始所說,當我們將$x_k$設為$w_n^k = e^{2\pi k / n}$時,會得到
$$
\begin{align}
y_k = A(x_k) = A(w_n^k) = \sum_{j=0}^{n-1} a_j w_n^{kj} \tag{1} \label{1}
\end{align}
$$
現在我們假設$n$是2的整數次方。下面我們將用divide-and-conquer的策略將問題拆解成兩個子問題。
我們可以觀察到
$$
\begin{align}
A(x) =& \ a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \cdots + a_{n-2} x^{n-2} + a_{n-1} x^{n-1} \\
=& \ a_0 + a_2 x^2 + a_4 x^4 + \cdots + a_{n-2} x^{n-2} \\
+& \ x(a_1 + a_3 x^2 + a_5 x^4 + \cdots + a_{n-1} x^{n-2})
\end{align}
$$
當我們令
$$
\begin{align}
A^{[0]}(x) = a_0 + a_2 x + a_4 x^2 + \cdots a_{n-2} x^{n / 2 -1}\\
A^{[1]}(x) = a_1 + a_3 x + a_5 x^2 + \cdots a_{n-1} x^{n / 2 -1}\\
\end{align}
$$
則上面的式子可以進一步表示成
$$
\begin{align}
A(x) = A^{[0]} (x^2) + x A^{[1]} (x^2) \tag{2} \label{2}
\end{align}
$$
也就是說,計算$A(x)$在$w_n^0, w_n^1, \ldots, w_n^{n-1}$上之值的問題可化簡為
1. 分別計算兩個degree-bound為n/2的多項式$A^{[0]}$和$A^{[1]}$在$(w_n^0)^2, (w_n^1)^2, \ldots, (w_n^{n-1})^2$上的值。
2. 利用式$(\ref{2})$將1.的兩個結果結合
由[Lemma 3](#Lemma-3-(Halving-lemma))可知上面的$n/2$個點$(w_n^0)^2, (w_n^1)^2, \ldots, (w_n^{n-1})^2$其實就是$w_{n/2}^0, w_{n/2}^1, \ldots, w_{n/2}^{n-1}$,也就是說就算$A^{[0]}(x^2)$和$A^{[1]}(x^2)$的值變成了兩個個和原問題$(\ref{1})$相同形式的子問題!(不同之處在於degree-bound減半為$n/2$)。也就是說,我們可以不斷的如此遞迴下去直到遇見degree-bound為$1$的最簡單子問題(此時輸入只有一個係數,只要將它回傳即可)。
#### 步驟
根據上面的想法,我們可以寫出以下的P最後演算法的pseudocode來計算$a = (a_0, a_1, \ldots, a_{n-1})$的傅立葉轉換$y$(維度為$n$的向量):
$$
\begin{align}
&\text{RECURSIVE-FFT}(a)&\quad & \quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad \\
1 & \quad n = a.length \\
2 & \quad \textbf{if } n == 1\\
3 & \quad \quad \textbf{return } a\\
4 & \quad w_n = e^{2\pi i/n}\\
5 & \quad w = 1\\
6 & \quad a^{[0]} = (a_0, a_2, \ldots, a_{n-2})\\
7 & \quad a^{[1]} = (a_1, a_3, \ldots, a_{n-1})\\
8 & \quad y^{[0]} = \text{RECURSIVE-FFT}(a^{[0]}) \quad \quad // \text{ array of size } n/2\\
9 & \quad y^{[1]} = \text{RECURSIVE-FFT}(a^{[1]})\\
10 & \quad \textbf{for } k = 0 \textbf { to } n/2 -1\\
11 & \quad \quad y_k = y_k^{[0]} + w y_k ^{[1]} \quad \quad // \ \text{array of size } {n}/{2}\\
12 & \quad \quad y_{k + (n/2)} = y_k^{[0]} - w y_k^{[1]}\\
13 & \quad \quad w = w w_n \quad \quad // \ \text{compute } w_n^{k+1} \text{for next iteration}\\
14 & \quad \textbf{return } y \quad \quad \quad // \ y \text{ is assumed to be a column vector} \\
\end{align}
$$
**說明**
上面的8,9行就是透過divide-and-conquer的方式分別計算兩個大小為$n/2$的傅立葉轉換,其中
$$
y_k^{[0]} = A^{[0]} (w_{n/2}^k), \\
y_k^{[1]} = A^{[1]} (w_{n/2}^k)
$$
且$k = 0, 1, \ldots, n/2 - 1$。而11,12行就是將兩個子問題的結果結合起來計算原本DFT的結果,第11行透過$(\ref{2})$計算了$y_0, y_1, \ldots, y_{n/2 - 1}$的值,而因為
$$
\begin{align}
y_{k + (n/2)} &= A(w_n^{k + (n/2)}) \\
&\overset{(\ref{2})}= A^{[0]}(w_n^{2k + n}) + w_n^{k + (n/2)} A^{[1]} (w_n^{2k + n}) \\
&= A^{[0]}(w_n^{2k}) + w_n^{k + (n/2)} A^{[1]} (w_n^{2k})
& (w_n^n = 1)\\
&= A^{[0]}(w_n^{k}) - w_n^{k} A^{[1]} (w_n^{2k}) & \text{(Lemma 2)} \\
&= y_k^{[0]} - w_n^k y_k^{[1]}
\end{align}
$$
也因此,$y_{k + n/2}$的計算可以跟$y_k$的計算合併在同一個iteration中,也就是第12行在做的事。而13行做的事情就是利用第$k$次iteration的$w_n^k$來計算出下一次iteration所需的$w_n^{k + 1}$,省去每次iteration都要重新計算$w_n^{k}$的計算成本。
#### 時間複雜度
根據上面的演算法,在for迴圈之中我們需要$\Theta(n)$次複數加法和乘法,也因此當我們令input vector的長度為$T(n)$時,該遞迴的時間複雜度為:
$$
T(n) = 2 T(n/2) + \Theta(n) = \Theta(n \log n)
$$
### 更有效率的FFT實作方式
透過一些FFT特殊的性質,我們可以將上一節的$\text{RECURSIVE-FFT}$改寫成非遞迴的形式,並且剩下一些計算及空間上的成本(演算法本身複雜度不變)。
首先,在$\text{RECURSIVE-FFT}$的for迴圈之中,我們重複計算了兩次$w_n^k y_k^{[1]}$,也因此我們可以透過暫存該運算式到變數$t$來減少一次乘法的運算,因此迴圈可改寫如下:
$$
\begin{align}
& \textbf{for } k = 0 \textbf{ to } n/2 - 1 &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad \\
& \quad t = w y_k^{[1]} \\
& \quad y_k = y_k^{[0]} + t \\
& \quad y_{k + (n/2)} = y_k^{[0]} - t \\
& \quad w = w w_n
\end{align}
$$
在這個迴圈之中我們將$w = w_n^k$乘上$y_k^{[1]}$,並將結果存進變數$t$,接著對$y_k^{[0]}$分別加上合減去$t$,這個運算被稱為***butterfly operation***,圖示如下:
![](https://i.imgur.com/QLqnASG.png =300x)
現在我們的目標是將FFT演算法從遞迴的形式轉換為迭代的形式,首先我們先觀察原先遞迴形式的呼叫結構,假設$n=8$,我們可以畫出下面的樹狀圖:
![](https://i.imgur.com/FRjI6YY.png =600x)
在此圖中,每個節點都代表一次$\text{RECURSIVE-FFT}$函數的呼叫,而節點中的陣列則是該次呼叫的輸入參數,由於該函式會向下遞迴呼叫兩次函數,圖中節點的left child代表的是第一次呼叫(第8行),right child代表的是第二次呼叫(第9行)。可以發現,只要我們適當的編排最一開始的輸入陣列$a$的元素使得元素的順序和圖中葉節點從左至右的順序相同(令此編排後的陣列為$A[0, \ldots, n-1]$),則我們可以**由下往上**運行FFT,而非由上往下遞迴呼叫。姑且我們先忽略這個順序是如何透過副函式$\text{BIT-REVERSE-COPY}$生成的,當我們要執行FFT時,是從最底部的節點先開始,而節點的合併就是執行上面描述的butterfly operatoin所獲得的兩倍長陣列,當底部第一層的節點都合併為$n/2$個第二層節點後,再繼續下一層的合併,以此類推直到最後合併出來的結果就是FFT的輸出。
而要實作出上面描述的演算法,我們可以引入變數$s$來表示目前迴圈所在的層數,則上述的演算法大致迭代結構如下:
$$
\begin{align}
& 1 \quad \textbf{for } s = 1 \textbf{ to } \log n \\
& 2 \quad \quad \textbf{for } k = 0 \textbf{ to } n - 1 \text{ by } 2^s \\
& 3 \quad \quad \quad \text{combine the two } 2^{s-1} \text{-element DFTs in} \\
& \ \quad \quad \quad \quad A[k, \ldots, k + 2^{s-1} -1] \text{ and }
A[k + 2^{s-1}, \ldots, k + 2^s -1] \\
& \ \quad \quad \quad \quad \text{into one } 2^s \text{-element DFT in} A[k, \ldots, k+ 2^s - 1]
\end{align}
$$
因此,最後演算法的pseudocode如下:
$$ \begin{align} & \text{ITERATIVE-FFT}(a) &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad \\
1 & \quad \text{BIT-REVERSE-COPY}(a,A) \\
2 & \quad n = a.length \quad \quad \quad \quad \quad \quad \quad \quad // \ n \text{ is a power of 2} \\
3 & \quad \textbf{for } s=1 \textbf{ to } \log n \\
4 & \quad \quad m = 2^s \\
5 & \quad \quad w_m = e^{2 \pi i/m} \\
6 & \quad \quad \textbf{for } k = 0 \textbf{ to } n-1 \textbf{ by } m \\
7 & \quad \quad \quad w = 1 \\
8 & \quad \quad \quad \textbf{for } j = 0 \textbf{ to } m/2 - 1 \\
9 & \quad \quad \quad \quad t = w A[k + j + m/2]\\
10 & \quad \quad \quad \quad u = A[k+j]\\
11 & \quad \quad \quad \quad A[k+j] = u + t \\
12 & \quad \quad \quad \quad A[k+j+m/2] = u -t \\
13 & \quad \quad \quad \quad w = w w_m\\
14 & \quad \textbf{return } A\\
\end{align}
$$
最後的重點就是$\text{BIT-REVERSE-COPY}(a, A)$是如何決定順序的,這點我們可以從上面的樹狀圖中觀察,當根節點分裂為兩個子節點時,在$a$中的元素是根據其index的奇偶性來分群的,也就是說,當我們將$a$中元素的index以二進制來表示時$(000, 001, \ldots)$,我們看的是最右邊的位元值,若值為$0$,則元素會被分到left child,反之若值為$1$,則元素會被分類至right child,而到下一層的分群時,則是看index的第二個位元,並且套用相同的規則,以此類推到最底層。因此當我們想判斷兩個葉節點的相對位置時,要先比較的是原先元素index的第一個位元大小,較小的元素會在另一個元素的左邊,而如果該位元值相同,則繼續朝下一個位元做比較,以此類推。也就是說,其實我們只要把元素的index位元**翻轉過來**,並且比較反轉後之二進制數的大小,較小者會在較大者的左邊,用這個方式我們就可以決定元素在底層葉節點的順序了:那就是將元素原本在$a$的index表示成二進制並進行翻轉,翻轉的結果就是該元素在$A$的index。寫成pseudocode的描述方式如下:
$$
\begin{align}
& \text{BIT-REVERSE-COPY}(a, A) &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad &\quad \\
1 & \quad n = a.length \\
2 & \quad \textbf{for } k = 0 \textbf{to } n-1 \\
3 & \quad \quad A[\text{rev}(k)] = a_k
\end{align}
$$
## 多維傅立葉轉換
一維傅立葉轉換可以被推廣到多維的情況,假設一個$d$維陣列$A = (a_{j_1, j_2, \ldots, j_d})$的維度是$n_1, n_2, \ldots, n_d$,且$n_1 n_2 \cdots n_d = n$。$n$維的傅立葉轉換被定義為
$$
\begin{align}
y_{k_1, k_2, \ldots, k_d} =
\sum_{j_1=0}^{n_1-1} \sum_{j_2=0}^{n_2-1} \cdots \sum_{j_d=0}^{n_d-1}
a_{j_1, j_2, \ldots, j_d} w_{n_1}^{j_1 k_1} w_{n_2}^{j_2 k_2} \cdots w_{n_d}^{j_d k_d}
\end{align}
$$
當我們想計算上面的式子時,可以先將$A$作為輸入計算陣列中沿著第一個維度方向上的一維傅立葉轉換(共有$n/n_1$個長度為$n_1$的一維向量),接著以其結果(維度和$A$相同)做輸入,沿著第二個維度的方向做一維傅立葉轉換(共有$n/n_2$個長度為$n_2$),以此類推,最後得到的結果就是整個$n$維傅立葉轉換的結果。可以發現,如果我們使用快速傅立葉轉換的話,這個方式的複雜度總共是
$$
\begin{align}
O(\frac{n}{n_1} \cdot n_1\log {n_1}) +
O(\frac{n}{n_2} \cdot n_2\log {n_2}) + \cdots
O(\frac{n}{n_d} \cdot n_3\log {n_d}) &=
O(n(\log n_1 + \log n_2 + \cdots + \log n_d)) \\ &= O(n \log n)
\end{align}
$$
也就是說,$n$維快速傅立葉轉換的複雜度只與元素的個素有關,而跟維度的數量$d$無關。
## 參考資料
[1] Introduction to Algorithms, 3rd Edition (The MIT Press) (3rd ed.). (2009). MIT Press.