傅里叶变换

波的基础知识¶

波是一种通过介质传递能量的现象，例如水中的涟漪或声音。我们可以使用随空间（$x$）和时间（$t$）变化的正弦函数来为一个简单的波（或正弦波）建模。本章介绍傅里叶变换，这是一个强大的工具，用于将复杂的波分解为简单正弦波的组合。在深入研究变换本身之前，让我们先建立一种更强大的方法来描述这些基本的振荡。

我们可以认为正弦波可以在时间和空间上都发生变化。如果我们绘制出不同位置的变化，每个时间快照都将是一个随位置变化的正弦波。请看下图，其中一个固定点在$x=2.5$处显示为一个红点。当然，你也可以看到特定位置随时间的变化，你可以自己绘制这个图。

从正弦波到复指数¶

为了理解傅里叶变换，我们首先需要一种比标准正弦或余弦函数更强大的方法来描述振荡。这通过使用复数并将其与圆周运动联系起来实现。

复数简介¶

一个复数 $y$ 是一个可以表示为以下形式的数：

其中 $a$ 和 $b$ 是实数，$i$ 是虚数单位，满足方程 $i^2 = -1$。

• $a$ 是实部 (real part)。

• $b$ 是虚部 (imaginary part)。

我们可以在一个称为复平面的二维平面上可视化复数，其中水平轴代表实部，垂直轴代表虚部。

$y = a + ib$ 的共轭复数 (complex conjugate) 记作 $\bar{y}$ 或 $y^*$，定义为 $\bar{y} = a - ib$。在几何上，它表示复数关于实轴的镜像。

通过圆周运动理解波的术语¶

虽然波通常用正弦函数描述，但它们的属性与匀速圆周运动密切相关。下图有助于在这种背景下定义关键的波的术语：

• 振幅 (Amplitude, A): 振动体或波上的点与其平衡位置之间的最大位移或距离。在圆周运动的背景下，它就是圆的半径。

• 频率 (Frequency, f): 单位时间内绕圆周运动的圈数（例如，每秒的圈数，以赫兹为单位）。

• 周期 (Period, T): 完成一次完整旋转所需的时间 ($T = 1/f$)。

• 角频率 (Angular Frequency, $\omega$): 旋转速率，以弧度/秒为单位 ($\omega = 2\pi f$)。

• 相位 (Phase, $\phi$): 描述了波形周期上某一点在时间上的位置。它指示了相对于参考点已过去的波周期的分数，通常用度或弧度表示。在此背景下，它是指在时间 $t=0$ 时点在圆上的起始角度。

这些术语共同构成了我们熟悉的正弦波方程：

然而，一种更强大的表示方法来自复平面。

欧拉公式及其与波的联系¶

圆周运动和波之间的联系通过欧拉公式得以形式化，该公式描述了复平面上单位圆上的一个点：

这个方程表明，复数 $e^{i\theta}$ 代表单位圆上角度为 $\theta$ 的一个点。其实部是 $\cos(\theta)$，虚部是 $\sin(\theta)$。

欧拉公式通过复数将指数函数与三角函数联系起来。来源: Britannica

为了表示具有特定振幅 $A$ 的波，我们对单位圆的表示进行缩放。半径为 $A$ 的圆上的一个点可以描述为一个复数 $y$：

使用欧拉公式，可以将其简化为紧凑的指数形式：

如果我们让角度随时间变化，使得 $\theta = \omega t + \phi$，那么 $y(t) = Ae^{i(\omega t + \phi)}$ 这一项就代表一个以角频率 $\omega$、振幅 $A$ 和初始相位 $\phi$ 连续旋转的点。

• 这个旋转点在实轴上的投影描绘出一条余弦波：$A\cos(\omega t + \phi)$。

• 在虚轴上的投影则描绘出一条正弦波：$A\sin(\omega t + \phi)$。

下面的动画直观地展示了这种基本联系。

我们甚至可以在三维空间中将这种关系可视化，其中复指数 $e^{i\theta}$ 随着时间描绘出一条螺旋线。这条螺旋线在实-时平面和虚-时平面上的投影分别揭示了余弦波和正弦波。

这种复指数表示法是傅里叶变换的基础。它使我们能够将任何简单的波（正弦波）表示为一个单一的旋转向量（“相量”），它优雅地结合了振幅和相位信息。这通过将困难的三角学转化为更易于管理的代数来简化复杂信号的分析。

生成两个正弦波，时间在0到1秒之间，频率分别为5Hz和10Hz，都以100Hz的频率采样。绘制这两个波并观察差异。计算1秒内有多少个周期。

生成两个正弦波，时间在0到1秒之间。两个波的频率都为5Hz，采样率为100Hz，但相位分别为0和10。此外，两个波的振幅分别为5和10。绘制这两个波并观察差异。

离散傅里叶变换 (DFT)¶

从上一节中，我们学习了如何用周期/频率、振幅、相位来轻松地描述一个波。但这些对于简单的周期性信号，如正弦波或余弦波，是很容易的。对于复杂的波，像那样描述就不容易了。例如，下面是一个相对更复杂的波，很难说出这个波的频率、振幅是多少，对吗？

现实世界中有更复杂的情况，如果我们有一种方法可以用来分析波的特性，那就太好了。傅里叶变换可以用于此目的，它将任何信号分解为简单的正弦和余弦波的总和，我们可以轻松地测量其频率、振幅和相位。傅里叶变换可以应用于连续或离散的波，在本章中，我们只讨论离散傅里叶变换（DFT）。

使用DFT，我们可以将上述信号分解为一系列正弦曲线，它们各自具有不同的频率。下面的动画演示了如何通过叠加正弦波（本轮）来重建复杂形状，这是傅里叶级数的核心思想。

下面的三维图显示了DFT背后的思想，即我们看到的上述信号实际上是3个不同正弦波相加的结果。时域信号，也就是我们看到的上述信号，可以转换成频域中的一个图，称为DFT振幅谱，其中信号频率显示为垂直条。归一化后，条的高度是时域中信号的振幅。你可以看到3个垂直条对应于正弦波的3个频率，这些频率也绘制在图中。

在本节中，我们将学习如何使用DFT来计算和绘制DFT振幅谱。

正向DFT¶

DFT可以将一系列均匀间隔的信号转换为关于所有正弦波频率的信息，这些正弦波需要相加以构成时域信号。其定义如下：

其中

• N = 采样点数

• n = 当前采样点

• k = 当前频率, 其中 $k\in [0,N-1]$

• $x_n$ = 采样点n处的正弦值

• $X_k$ = DFT，包含振幅和相位信息

此外，上述方程中的最后一个表达式源自欧拉公式，它将三角函数与复指数函数联系起来：$e^{i\cdot x} = cosx+i\cdot sinx$

由于变换的性质，$X_0 = \sum_{n=0}^{N-1}x_n$。如果$N$是奇数，元素$X_1, X_2, ..., X_{(N-1)/2}$包含正频率项，而元素$X_{(N+1)/2}, ..., X_{N-1}$包含负频率项，按负频率递减的顺序排列。如果$N$是偶数，元素$X_1, X_2, ..., X_{N/2-1}$包含正频率项，而元素$X_{N/2},...,X_{N-1}$包含负频率项，按负频率递减的顺序排列。在我们的输入信号$x$是实值序列的情况下，正频率的DFT输出$X_n$是负频率值$X_n$的共轭，频谱将是对称的。因此，通常我们只绘制对应于正频率的DFT。

请注意，$X_k$是一个复数，它编码了函数$x_n$的复正弦分量$e^{i\cdot 2\pi kn/N}$的振幅和相位信息。信号的振幅和相位可以计算如下：

其中$Im(X_k)$和$Re(X_k)$是复数的虚部和实部，$atan2$是$arctan$函数的双参数形式。

如果我们用采样点数对DFT返回的振幅进行归一化，它们就等于输入到DFT的信号的振幅。请注意，这样做会将功率分配在正负两侧，如果输入信号是如上所述的实值序列，正负频率的频谱将是对称的，因此，我们只看DFT结果的一侧，并且我们用$N/2$而不是$N$来除，以获得对应于时域信号的振幅。

现在我们对DFT有了基本的了解，让我们看看如何使用它。

生成3个频率为1Hz、4Hz和7Hz，振幅为3、1和0.5，相位均为零的正弦波。以100Hz的采样率将这3个正弦波相加，你会看到它与我们在本节开头展示的信号相同。

编写一个函数DFT(x)，它接受一个参数，x - 输入的一维实值信号。该函数将计算信号的DFT并返回DFT值。将此函数应用于我们上面生成的信号并绘制结果。

我们可以从这里看到，DFT的输出在采样率的一半处是对称的（你可以尝试不同的采样率来测试）。采样率的这一半被称为奈奎斯特频率或折叠频率，它是以电子工程师哈里·奈奎斯特的名字命名的。他和克劳德·香农提出了奈奎斯特-香农采样定理，该定理指出，如果一个信号只包含低于采样频率一半的频率分量，那么以该速率采样的信号可以被完全重建，因此DFT输出的最高频率是采样率的一半。

通过绘制DFT结果的前半部分，我们可以看到在频率1Hz、4Hz和7Hz处有3个清晰的峰值，振幅分别为3、1、0.5，符合预期。这就是我们如何使用DFT通过将任意信号分解为简单的正弦波来分析它。

逆DFT¶

当然，我们可以轻松地进行DFT的逆变换。

我们将把它作为一个练习留给你来编写一个函数。

DFT的局限性¶

上述DFT实现的主要问题是，如果我们有一个包含许多数据点的信号，它的效率就不高。如果信号很大，计算DFT可能需要很长时间。

试一试 编写一个函数来生成一个具有不同采样率的简单信号，并通过改变采样率来看计算时间的差异。

我们可以看到，随着数据点数量的增加，使用这种DFT会花费大量的计算时间。幸运的是，Cooley和Tukey在他们1965年的论文中推广了快速傅里叶变换（FFT），有效地解决了这个问题，这将是下一节的主题。

快速傅里叶变换 (FFT)¶

**快速傅里叶变换（FFT）**是一种计算序列DFT的高效算法。它最早在Cooley和Tukey于1965年发表的经典论文中被描述，但这个想法实际上可以追溯到高斯在1805年未发表的工作。它是一种分治算法，通过递归地将DFT分解为更小的DFT来减少计算量。因此，它成功地将DFT的复杂度从$O(n^2)$降低到$O(nlogn)$，其中$n$是数据的大小。这种计算时间的减少是显著的，特别是对于具有大$N$的数据，因此使得FFT在工程、科学和数学中被广泛使用。FFT算法被《科学与工程计算》杂志评为20世纪十大算法之一。

在本节中，我们将向您介绍FFT如何减少计算时间。本节内容主要基于Jake VanderPlas编写的这篇优秀的教程。

DFT中的对称性¶

FFT如何加速DFT计算的答案在于利用了DFT中的对称性。让我们来看一下DFT中的对称性。根据DFT方程的定义

我们可以计算出

注意，$e^{-i2\pi{n}} = 1$，因此，我们有

稍作扩展，我们可以得到

这意味着在DFT内部，我们显然有一些可以用来减少计算的对称性。

FFT中的技巧¶

既然我们知道DFT中存在对称性，我们可以考虑利用它来减少计算，因为如果我们需要同时计算$X_k$和$X_{k+N}$，我们只需要计算一次。这正是FFT背后的思想。Cooley和Tukey表明，如果我们继续将问题分解成更小的问题，我们可以更有效地计算DFT。让我们首先将整个序列分成两部分，即偶数部分和奇数部分：

我们可以看到，上述方程中两个较小的项（大小只有$\frac{N}{2}$）是两个较小的DFT。对于每一项，$0\leq m \le \frac{N}{2}$，但$0\leq k \le N$，因此，我们可以看到由于我们上面描述的对称性，一半的值将是相同的。因此，我们只需要计算每一项中一半的字段。当然，我们不必止步于此，我们可以继续将每一项的偶数和奇数值分成一半，直到它达到最后两个数字，那时的计算将非常简单。

这就是FFT算法的本质：递归地将DFT分解为更小的DFT，直到达到基本情况。下面的伪代码说明了这种分治策略：

Cooley-Tukey FFT算法伪代码：

function FFT(x):
    N = length(x)
    
    // 基本情况：如果只有一个元素，直接返回
    if N == 1:
        return x
    
    // 分治：分离成偶数索引和奇数索引元素
    x_even = x[0], x[2], x[4], ..., x[N-2]
    x_odd  = x[1], x[3], x[5], ..., x[N-1]
    
    // 递归：对每一半递归计算FFT
    X_even = FFT(x_even)
    X_odd  = FFT(x_odd)
    
    // 合并：使用旋转因子合并结果
    for k = 0 to N-1:
        factor_k = exp(-2πi * k / N)
        X[k] = X_even[k mod N/2] + factor_k * X_odd[k mod N/2]
    
    return X

关键要点：

• 基本情况：当N=1时，单个点的DFT就是该点本身

• 分治步骤：将输入分成偶数和奇数索引元素

• 递归步骤：递归计算每一半的FFT（大小为N/2）

• 合并步骤：使用复指数"旋转因子" $e^{-2\pi i k/N}$ 合并结果

• 对称性：由于周期性，$X_{even}[k]$ 和 $X_{odd}[k]$ 对所有k值循环使用

现在让我们看一个这个递归FFT算法的Python实现。重要提示：该实现要求输入信号长度必须是2的幂（即N = 2, 4, 8, 16, 32, ...）。如果您的信号长度不是2的幂，应该用零填充到下一个2的幂。例如，长度为100的信号应该填充到128。

使用FFT函数计算上述信号的傅里叶变换。绘制双边和单边频率的振幅谱。

生成一个长度为2048的简单信号，并计时运行FFT所需的时间，并与DFT的速度进行比较。

我们可以看到，对于长度为2048（约2000）的信号，这种FFT的实现使用了16.9毫秒，而不是DFT的120毫秒。请注意，还有很多方法可以优化FFT的实现，使其更快。

我们可以看到，随着数据点数量的增加，使用这种DFT会花费大量的计算时间。幸运的是，Cooley和Tukey在他们1965年的论文中推广了快速傅里叶变换（FFT）。

二维傅里叶变换¶

（二维）傅里叶变换是图像处理中一个非常经典的工具。 它是一维信号傅里叶变换的扩展， 将信号分解为复振荡（实际上是复指数）的和。 在图像处理中，傅里叶变换将图像分解为具有不同频率、相位和方向的振荡的和。 注意，如果像素是实数值，则振荡不是复指数。

因此，傅里叶变换提供了关于图像频率内容的信息， 即图像中的强度如何通过不同频率分布。 Figure 6 显示了将合成图像分解为振荡。 在这个玩具示例中，图像足够简单，仅使用三个振荡就可以分解。 我们将进一步看到，通常的图像需要更多的振荡。

正向傅里叶变换¶

大小为 $M \times N$ 的图像 $f$ 的离散傅里叶变换（DFT）是一个同样大小的图像 $F$，定义为：

接下来，我们用 $\mathcal{F}$ 表示 DFT，因此 $\mathcal{F}[f] = F$。

请注意，傅里叶变换的定义使用了一个复指数。 因此，图像的 DFT 可能是复数的，因此不能用单个图像显示。 这就是为什么我们将分别显示 DFT 的振幅（模）和相位（参数），如 Figure 7 所示。

计算振幅和相位：

由于 $F(u,v)$ 是一个复数，我们可以按以下方式提取其振幅和相位：

其中：

• $Re(F(u,v))$ 是复数DFT系数的实部

• $Im(F(u,v))$ 是复数DFT系数的虚部

• $|F(u,v)|$ 表示振幅（或幅度）谱

• $\angle F(u,v)$ 表示相位谱

• $\text{atan2}$ 是双参数反正切函数，它能正确处理所有四个象限

振幅谱 $|F(u,v)|$ 告诉我们图像中存在的每个频率分量的强度，而相位谱 $\angle F(u,v)$ 包含关于这些频率分量空间位置的信息。两者对于从傅里叶变换重建原始图像都是必不可少的。

振幅和相位表示频率平面中能量的分布。 低频的强度位于图像的中心，高频的强度位于边界附近。

在上图中，松鼠后面的灰色背景主要“包含”低频，因为像素的强度从一个像素到另一个像素缓慢演变。 相反，尾部需要高频来显示毛发和背景之间的快速交替。

逆傅里叶变换¶

逆离散傅里叶变换从其傅里叶变换计算原始图像：

下文中用 $\mathcal{F}^{-1}$ 表示。

示例：实践中的二维傅里叶变换¶

让我们看一个使用 Python 和 NumPy/SciPy 库将二维傅里叶变换应用于图像的完整示例。

示例要点¶

此示例演示了几个重要概念：

1. 计算二维DFT：使用 SciPy 的 fft2() 计算二维图像的傅里叶变换。

2. 频率移位：fftshift() 函数将零频率分量（DC分量）移动到频谱的中心，使可视化更加直观。

3. 振幅和相位提取：

   • 振幅：np.abs(F) 计算 $|F(u,v)| = \sqrt{Re(F)^2 + Im(F)^2}$

   • 相位：np.angle(F) 计算 $\angle F(u,v) = \text{atan2}(Im(F), Re(F))$

4. 对数刻度：振幅谱具有非常大的动态范围。应用 log(1 + amplitude) 有助于可视化大小值。

5. 频域滤波：该示例展示了一个简单的低通滤波器：

   • 在频域创建圆形掩码

   • 仅保留低频（频谱中心）

   • 去除高频（边缘和细节）

   • 产生模糊图像

6. 逆变换：使用 ifft2() 和 ifftshift() 从频域转换回空间域。

这个工作流程（正向FFT → 频率操作 → 逆FFT）是许多图像处理应用的基础，包括去噪、锐化和压缩。

属性¶

• DFT是线性的：

  $$\mathcal{F}[af + bg] = aF + bG \qquad\text{其中}\; a,b\in\mathbb{C}.$$

  (19)

• 两个图像的卷积等同于图像 DFT 的乘积， 前提是卷积是循环的（边缘上的环绕假设）：

  $$f * g = \mathcal{F}^{-1}[F \times G]$$

  (20)

• DFT是可分离的：它可以通过在行上计算一维 DFT，然后在列上计算一维 DFT 来获得。

• DFT在周期 $M$ 和 $N$ 上是周期的（$k, l \in \mathbb{Z}$）：

  $$F(u,v) = F(u+kM,v) = F(u,v+lN) = F(u+kM,v+lN).$$

  (21)

• 图像上的平移意味着 DFT 上的相移：

  $$\mathcal{F}[f(m-m_0,n-n_0)] = F(u,v) \exp\left(-j2\pi\left(\frac{um_0}{M}+\frac{vn_0}{N}\right)\right)$$

  (22)

• 图像上的旋转意味着 DFT 上的相同旋转。