并行谱数值方法/一维离散傅里叶变换

离散傅里叶变换（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 算法在许多地方都有很好的记录，但是，还有其他算法实现，要使用的最佳算法版本在很大程度上取决于计算机体系结构。因此，我们这里不再做进一步的描述。