快速傅立葉轉換 (FFT)

--- 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.