並行譜數值方法/一維離散傅立葉變換

離散傅立葉變換（DFT）採用在有限個點處取樣的函式，並找到最接近該函式的三角多項式線性組合的係數；使用的三角多項式的數量等於樣本點的數量。假設我們有一個在區間 $a\leq x\leq b$ 上定義的函式 $f(x)$ 。由於記憶體限制，計算機只能儲存有限個樣本點的值，即 $a\leq x_{0}<x_{1}<...<x_{n}\leq b$ 。出於我們的目的，這些點將是等間距的，例如 $x_{1}-x_{0}=x_{3}-x_{2}$ ，因此我們可以寫出

$x_{j}=a+jh,\qquad j=0,1,2,...,n$

其中 $x_{j}$ 是樣本點， $n$ 是樣本點的數量，並且

$h={\frac {b-a}{n}}.$

使用標準區間很方便，對於該區間 $0\leq x\leq 2\pi$ 。用標準區間重寫 $x$ 得到

$x_{0}=0,x_{1}={\frac {2\pi }{n}},x_{2}={\frac {4\pi }{n}},x_{j}={\frac {2j\pi }{n}},...,x_{n-1}={\frac {2(n-1)\pi }{n}}$

注意如何省略 $x_{n}=2\pi$ ；週期性意味著函式在 $2\pi$ 處的函式值與函式在 $0$ 處的函式值相同，因此無需包含。我們將使用線性代數的語言介紹 DFT。這種形式主義的很多內容適用於正在被逼近的連續函式。它也使理解演算法的計算機實現變得更容易。許多計算機包和程式都針對透過矩陣運算執行計算進行了最佳化，因此這種形式主義在實際計算變換時也很有用。我們將對 $f(x)$ 在樣本點處的逼近寫成一個有限維向量

${f}=(f_{0},f_{1},...,f_{n-1})^{T}=(f(x_{0}),f(x_{1}),...,f(x_{n-1}))$

其中

$f_{j}=f(x_{j})=f\left({\frac {2j\pi }{n}}\right).$

DFT 將取樣函式 $f(x)$ 分解成復指數的線性組合， $\exp(ikx)$ ，其中 $k$ 是一個索引。由於

$\exp(ikx)=\cos(kx)+i\sin(kx),$

我們還獲得了三角函式的展開，這在微積分和微分方程課程中可能更為熟悉。由於該函式在 $n$ 個點處取樣，因此可以分辨出的最高振盪頻率將具有 $n$ 個振盪。原始函式中任何高於 $n$ 的頻率都無法得到充分的分辨，並導致混疊誤差（例如，參見 Boyd^[1] 或 Uecker^[2] 以瞭解更多資訊）。透過在更多點處取樣，可以減少這種誤差，從而也可以增加逼近指數函式的數量。在提高模擬精度和模擬完成所需時間之間存在權衡。對於許多科學和實際關注的情況，最多包含數千個網格點的模擬可以在相對較短的時間內計算出來。下面我們將解釋如何透過插值三角多項式 $p(x)$ 來逼近函式 $f(x)$ ，使得

$f(x)\approx p(x)=c_{0}+c_{1}e^{2ix}+c_{2}e^{2ix}+...+c_{n-1}e^{(n-1)ix}=\sum \limits _{k=0}^{n-1}c_{k}e^{ikx}$

符號 $\approx$ 表示 $f(x)$ 和 $p(x)$ 在每個樣本點上都一致，即對於每個 $j=0,1,...n-1$ ，都有 $f(x_{j})=p(x_{j})$ ，但插值多項式 $p(x)$ 只是真實解 $f(x)$ 在樣本點以外的近似值。 $c_{n}$ 稱為離散的 傅立葉係數，是我們需要求解的目標。 $p(x)$ 代表了度數為 $\leq n-1$ 的插值三角多項式的值，所以如果我們有這些係數的值，那麼我們就有一個可以使用的函式。由於我們在有限維向量空間中工作，一個有用的方法是將離散傅立葉級數重寫為一個向量。令

$\omega _{k}=(e^{ikx_{0}},e^{ikx_{1}},e^{ikx_{2}},...,e^{ikx_{n}})^{T}$

$=(1,e^{2k\pi i/n},e^{4k\pi i/n},...,e^{2(n-1)k\pi i/n})^{T},$

其中 $k=0,1,...,n-1$ 。插值條件， $f(x_{j})=p(x_{j})$ ，也可以用向量形式重寫

f=c_{0}{\omega _{0}}+c_{1}{\omega _{1}}+...+c_{n-1}{\omega _{n-1}}.

( 1)

這裡 $f$ 是在取樣點處評估的向量，它被分解成向量 ${\omega }_{k}$ ，就像三維空間中的向量可以分解成 $x$ ， $y$ 和 $z$ 方向上的分量。離散傅立葉變換允許我們計算係數 $c_{i}$ ，只要給定函式在取樣點的值。這可能最初看起來沒有動機，但在許多應用中，例如求解微分方程，操作三角多項式的線性組合， ${\omega _{0}},...,{\omega _{n-1}}$ ，比處理原始函式更容易。為了求解 $c_{k}$ ，我們利用基元素 ${\omega _{0}},...,{\omega _{n-1}}$ 的正交性。我們現在解釋如何做到這一點（有關更詳細的解釋，請參見 Olver 和 Shakiban^[3]）。

定義 $\xi _{n}=e^{2\pi i/n}$ 。我們觀察到

$\left(\xi _{n}\right)^{n}=\exp \left({\frac {2\pi in}{n}}\right)=\cos(2\pi )+i\sin(2\pi )=1$

因此， $\xi _{n}$ 被稱為單位根的本原 $n^{\text{th}}$ 次方。還要注意的是，對於 $0\leq k<n$ ，我們有 $(\xi _{n}^{k})^{n}=1$ ，所以所有其他單位根，當取到 $n$ 次方時，可以從單位根的本原 $n^{\text{th}}$ 次方得到。我們將使用它來執行 DFT 演算法，計算公式 1 中的係數 $c_{0},...,c_{k-1}$ 。DFT 演算法背後的主要思想是利用向量 $\omega _{k}$ 的正交性。為了證明向量 $\omega _{k}$ 和 $\omega _{l}$ 之間的正交性，我們用 $\omega _{l}^{*}$ 表示 $\omega _{l}$ 的複共軛，然後取 $\omega _{k}$ 和 $\omega _{l}$ 的內積，發現

$<\omega _{k},\omega _{l}>$

$={\frac {1}{n}}\sum _{m=0}^{n-1}\exp \left({\frac {2\pi ikm}{n}}\right)\left[\exp({\frac {2\pi ilm}{n}})\right]^{*}$

$={\frac {1}{n}}\sum _{m=0}^{n-1}\exp \left({\frac {2\pi i(k-l)m}{n}}\right)$

$={\frac {1}{n}}\sum _{m=0}^{n-1}\cos \left({\frac {\pi (k-l)m}{n}}\right)+i\sin \left({\frac {\pi (k-l)m}{n}}\right)$

$=1$ 如果 $k=l$ 並且 $0$ 在其他情況下。

為了推匯出最後部分，如果 $k=l$ 那麼 $\exp(0)=1$ ，並且如果 $k\neq l$ ，那麼我們正在對正弦和餘弦函式進行取樣，這些函式在整個積分波長上的等間隔點。由於這些函式具有相等的幅度正負部分，因此它們之和為零，就像正弦或餘弦在積分波長上的積分等於零一樣。這意味著我們可以透過取內積來計算離散傅立葉和中的傅立葉係數

$c_{k}=<{f},{\omega _{k}}>={\frac {1}{n}}\sum \limits _{m=0}^{n-1}\xi _{n}^{-mk}f_{j}.$

我們注意到連續設定和離散設定之間的密切聯絡，其中積分被求和所取代。

快速傅立葉變換

使用 DFT 從定義中計算傅立葉係數 $c_{0},...,c_{n-1}$ 當 $n$ 很大時，速度可能會非常慢。計算傅立葉係數 $c_{0},...c_{n-1}$ 需要 $n^{2}-n$ 次複數乘法和 $n^{2}-n$ 次複數加法。1960 年，Cooley 和 Tukey ^[4] 重新發現了一種計算 DFT 的更高效方法，他們開發了一種稱為快速傅立葉變換 (FFT) 的演算法。FFT 將算術運算的次數降低到 $O(n\log n)$ 。對於 $n$ 的較大值，與標準 DFT 相比，這可以使計算時間產生巨大的差異。FFT 如此重要的原因是它被廣泛用於譜方法。Cooley 和 Tukey ^[5] 使用的基本 FFT 演算法在許多地方都有很好的記錄，但是，還有其他演算法實現，要使用的最佳演算法版本在很大程度上取決於計算機體系結構。因此，我們這裡不再做進一步的描述。