FFT（快速傅里葉變換）算法詳解

多項式的點值表示(Point Value Representation)

設多項式的系數表示(Coefficient Representation)：

\[\begin{align*} \mathrm P_a(x)&=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \\ &= \sum_{i=0}^{n-1}a_ix^i \end{align*} \]

則我們對上面的式子可以代入不同的 \(n\) 個 \(x\) 的值，構成一個 \(n\) 維向量：

\[\begin{bmatrix} \mathrm P_a(x_0) \\ \mathrm P_a(x_1) \\ \mathrm P_a(x_2) \\ \vdots \\ \mathrm P_a(x_{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{x-1}^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{bmatrix} \]

更簡潔的寫法：

\[\mathbf P_a = \mathbf X \mathbf \alpha \]

對上式觀察后發現，\(\mathbf X\) 是所謂的范德蒙德矩陣(Vandermonde's Matrix)，在 \(n\) 個 \(x\) 的值不同的情況下，其行列式的值為：

\[\det(\mathbf X) = \prod_{0\leqslant i\lt j \leqslant n-1} (x_j - x_i) \]

很明顯，當所有 \(n\) 個 \(x\) 取值不同時，其行列式不為零，因此 \(\mathbf X\) 可逆。
所以我們可以唯一確定多項式系數構成的向量 \(\alpha\)：

\[\mathbf \alpha = \mathbf X^{-1} \mathbf P_a \]

也就是說，多項式 \(\mathrm P_a(x)\) 還可以由 \(n\) 個 \(x\) 代入得到的 \(n\) 個點值來唯一表示：

\[\{ \left [x_0, \mathrm P(x_0) \right ], \left [ x_1, \mathrm P(x_1) \right ], \left [x_2, \mathrm P(x_2) \right ], \cdots,\left [ x_{n-1}, \mathrm P(x_{n-1}) \right ] \} \]

這就是多項式的點值表示。

多項式的點值表示是指，對於 \(n\) 次多項式，可以用 \(n\) 個不同的 \(x\) 和與之對應的多項式的值 \(\mathrm P(x)\) 構成一個長度為 \(n\) 的序列，這個序列唯一確定多項式，並且能夠與系數表示相互轉化。

\(n\) 次單位根

了解了多項式的點值表示，一個很自然的問題是：如何選擇 \(x\) 的值，來防止其指數大小爆炸型增長呢？這里可以借用復數的單位根。
簡單回顧一下，復數有兩種表示方法：迪卡爾積坐標表示和極坐標表示，這里我們用到的是后者：

\[z=re^{i\theta} \]

\(i\) 是虛數單位，\(r\) 表示模長，\(\theta\) 表示相角。
復數的 \(n\) 次單位根 \(\omega\) 需要滿足條件：

\[\omega^n_n=1 \]

了解復數乘法及其幾何意義的同學知道，復數相乘則相角^[1]相加，相當於復數點逆時針轉動；\(n\) 個復數相乘則說明有 \(n\) 個相角相加，\(n\) 次逆時針轉動。因此：

\[\omega_n^n = 1 = e^{i \cdot 2\pi} \]

則 \(n\) 次單位根為：

\[\sqrt[n] {\omega^n} = e^{i \cdot \frac{2\pi}{n}} \]

很容易想到，\(n\) 等分 \(2\pi\) 相當於 \(n\) 等分圓。下圖是 \(n=16\) 的例子：

引入了 \(n\) 次單位根，我們就可以找到任意 \(n\) 個不同的點值 \(x\)，並且不用擔心指數增長帶來的大小爆炸性增長的問題。

離散傅里葉變換(Discrete Fourier Transform)及其反變換

DFT

設長度為 \(n\) 的離散序列 \(\alpha^{\mathrm T}=[a_0, a_1, \cdots,a_{n-1}]\)，構建多項式：

\[\mathrm P_a(x) = \sum_{i=0}^{n-1}a_ix^i = \mathbf x \mathbf \alpha^{\mathrm T} \]

其中，\(\mathbf x^{\mathrm T} = [1, x, x^2, \cdots, x^{n-1}]\)

用 \(n\) 次單位根生成 \(n\) 個不同的值：\(\omega_n^0,\ \omega_n^1,\ \omega_n^2,\ \cdots,\ \omega_n^{n-1}\)
對應的點值表示可以用下面的矩陣運算完成：

\[\begin{bmatrix} \mathrm P_a(\omega_n^0) \\ \mathrm P_a(\omega_n^1) \\ \mathrm P_a(\omega_n^2) \\ \vdots \\ \mathrm P_a(\omega_n^{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{bmatrix} \]

更簡潔的寫法：

\[\mathbf P_a = \Omega \alpha \]

序列 \(\mathbf P_a\) 被稱為 \(\alpha\) 的離散傅里葉變換(DFT)，也稱為頻域序列。
很明顯，DFT 的計算時間復雜度是 \(O(n^2)\)

IDFT

離散傅里葉反變換，就是根據 DFT 得到的頻域序列算出多項式的系數（也稱時域序列）。這可以根據單位根矩陣 \(\Omega\) 的逆矩陣 \(\Omega^{-1}\) 得到

\[\alpha = \Omega^{-1}\mathrm P_a \]

一般來說，計算矩陣的逆的時間復雜度高達 \(O(n^3)\)。所幸，單位根矩陣的逆 \(\Omega^{-1}\) 有個優良的性質可以省去龐大的計算量：

\[\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix}^{-1} = \frac 1 n \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)} \\ 1 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)} \end{bmatrix} \]

也就是說，求 \(\Omega^{-1}\) 就是將里面元素的指數改成其相反數，再對所有元素除以 \(n\)。
有了這個性質，我們就可以得到：

\[\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{bmatrix}=\frac1 n \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)} \\ 1 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} \mathrm P_a(\omega_n^0) \\ \mathrm P_a(\omega_n^1) \\ \mathrm P_a(\omega_n^2) \\ \vdots \\ \mathrm P_a(\omega_n^{n-1}) \end{bmatrix} \]

很明顯，DFT 和 IDFT 的計算時間復雜度都是 \(O(n^2)\)，效率並不高。

快速傅里葉變換(Fast Fourier Transform)及其反變換

FFT

FFT 是用分治法（Divide & Conquer）的思想，用來優化 DFT 計算矩陣相乘的時間復雜度過高這一問題的算法。
設 \(n\) 次多項式 \(\mathrm P(x)\)^[2]：

\[\mathrm P(x)=a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \]

我們把多項式 \(\mathrm P(x)\) 按照系數下標的奇偶性分成兩部分

\[\mathrm P_{e}(x^2) = a_0 + a_2x^2 + a_4\left(x^2\right)^2 + \cdots + a_{n-2}x^{n-2} \\ x\mathrm P_{o}(x^2) = x\left [a_1 + a_3x^2+a_5\left(x^2\right)^2+\cdots+a_{n-1}x^{n-2} \right ] \]

則多項式 \(\mathrm P(x)\) 就是奇偶兩部分之和

\[\mathrm P(x)=\mathrm P_e(x^2)+x\mathrm P_o(x^2) \]

從上式中可以看出，多項式 \(\mathrm P(x)\) 可以由兩個系數個數為 \(n/2\) 的多項式 \(\mathrm P_{e}(x^2)\) 和 \(\mathrm P_o(x^2)\) 計算得到。
對於 \(\mathrm P_e(x^2)\) 和 \(\mathrm P_{o}(x^2)\)，我們令 \(x=x^2\)，就會發現這很明顯是一個遞歸的過程：

\[\mathrm P_{e}(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n-2}{2}} \\ \mathrm P_{ee}(x^2)=a_0+a_4x^2+a_8x^4+\cdots+a_{n-4}x^{\frac{n-2}{2}-1} \\ x\mathrm P_{eo}(x^2)=x\left [ a_2+a_6x^2+a_{10}x^4+\cdots+a_{n-2}x^{\frac{n-2}{2}-1} \right ] \\ \]

就可以求出

\[\mathrm P_e(x)=\mathrm P_{ee}(x^2)+xP_{eo}(x^2) \]

同理

\[\mathrm P_{o}(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n-2}{2}} \\ \mathrm P_{oe}(x^2)=a_1+a_5x^2+a_9x^4+\cdots+a_{n-3}x^{\frac{n-2}{2}-1} \\ x\mathrm P_{oo}(x^2)=x\left [ a_3+a_7x^2+a_{11}x^4+\cdots+a_{n-1}x^{\frac{n-2}{2}-1} \right ] \]

同樣可以求出

\[\mathrm P_o(x)=\mathrm P_{oe}(x^2)+xP_{oo}(x^2) \]

遞歸的終止條件很明顯，就是當遇到多項式中只有一個系數時返回該系數。
因此我們將 \(n\) 個單位根 \(\omega_n^0,\ \omega_n^1,\ \omega_n^2,\ \cdots,\ \omega_n^{n-1}\) 帶入多項式 \(\mathrm P(x)\)：

\[\mathrm P(\omega_n^{k})=\mathrm P_e(\omega_n^{2k})+\omega_n^k\mathrm P_o(\omega_n^{2k}) \quad (k=0, 1, \cdots,n-1) \]

剛剛提到，多項式 \(\mathrm P_{e}(x^2)\) 和 \(\mathrm P_o(x^2)\) 的系數個數為 \(n/2\) ，恰好 \(n\) 個單位根平方后也只剩下 \(n/2\) 個不同的單位根，簡單證明如下：

首先證明：

\[\left( \omega_n^k \right)^2=e^{i\frac{2\pi k}{n}\cdot2}=e^{i\frac{2\pi k}{n/2}}=\omega_{n/2}^k \]

因此相當於 \(n\) 等分圓變成了 \(n/2\) 等分圓。下面證明 \(k=0,\ 1,\ \cdots,\ \frac n 2-1\)：
設 \(k_1=m\), \(k_2=m+n/2\)，\(m=0,1,\cdots,\frac{n}{2}-1\)。則

\[\omega_{n/2}^{k_{1}} = \omega_{n/2}^{m}=e^{i\frac{2\pi m}{n/2}} \\ \omega_{n/2}^{k_2} = \omega_{n/2}^{m+n/2}=e^{i\frac{2\pi (m+n/2)}{n/2}}=e^{i\left(\frac{2\pi m}{n/2}+2\pi{}\right)}=e^{i\frac{2\pi m}{n/2}} \]

即：

\[\omega_{n/2}^{m}=\omega_{n/2}^{m+n/2}\quad (m=0,1,\cdots,\frac n 2-1) \]

由此可以說明，\(\omega_n^k\) 是周期序列，其最小正周期為 \(n/2\)。因此 \(k=0, 1, \cdots, \frac n 2 -1\)

因此多項式 \(\mathrm P(\omega_n^{k})\) 可改寫為

\[\mathrm P(\omega_n^{k})=\mathrm P_e(\omega_{n/2}^{k})+\omega_n^k\mathrm P_o(\omega_{n/2}^{k}) \quad (k=0, 1, \cdots,\frac n 2-1) \]

從上式中可以看出，計算 \(\mathrm P_e(\omega_{n/2}^{k})\) 和 \(\mathrm P_o(\omega_{n/2}^{k})\) 各需要 \(n/2\) 次乘法運算，即一共 \(n\) 次乘法運算；而計算\(\mathrm P_{ee}(\omega_{n/4}^{k})\) 、 \(\mathrm P_{eo}(\omega_{n/4}^{k})\) 、\(\mathrm P_{oe}(\omega_{n/4}^k)\) 和 \(\mathrm P_{oo}(\omega_{n/4}^k)\) 各需要 \(n/4\) 次乘法運算，即一共 \(n\) 次乘法運算……以此類推，這種每層遞歸規模減半的遞歸深度很明顯是 \(\log n\)，因此 FFT 算法的時間復雜度就是 \(O(n\log n)\)

IFFT

快速傅里葉變換反變換同樣是優化 IDFT 計算矩陣相乘的時間復雜度。由於 DFT 和 IDFT 核心操作一樣，都是矩陣相乘，因此 FFT 和 IDFT 的核心操作就是利用分治的思想減少矩陣相乘的運算量。可以想到，FFT 的過程可以直接移植到 IFFT 中來，需要修改的兩個地方是

1. 單位根指數部分改成其相反數：\(\omega_n^{-k}\)
2. 得到的結果均除以 \(n\)

因此，IFFT 的時間復雜度也是\(O(n\log n)\)

---

Refer: COMP3121/9101, CSE UNSW
Written with StackEdit.

---

1. 本文相角一律用弧度表示。 ↩︎

2. 這里為了推導過程的整潔，假設 \(n\) 是 2 的整數次冪。這個假設的合理性在於：即使 \(n\) 不是 2 的整數次冪，我們也可以在多項是后面后面補零達到離 \(n\) 最近的下一個 2 的整數次冪。 ↩︎