傅里叶变换

傅里叶的原理说明：任何连续测量的时序或信号，都可以用不同频率的正弦波信号的无限叠加来表示。

在数学的角度来看，傅里叶变换算法利用直接测量得到的初始信号，用累加的方式来计算该初始信号中不同正弦波信号的 频率、幅值和相位。

而从物理学的角度来看，傅里叶变换可以帮助我们 将时域的信号转为频域 来分析。有些在时域上很难分析的一些信号，在被转换到频域后，更容易找到其特征。此外，频谱 也可以通过傅里叶变换得出。下图可以更好理解时域到频域的转换：

FFT.png

Figure1 从时域到频域的变化

傅里叶变换的4种类型

  1. 非周期性连续信号傅里叶变换（Fourier Transform）
  2. 周期性连续信号傅里叶级数 (Fourier Series)
  3. 非周期性离散信号离散时域傅里叶变换（Discrete Time Fourier Transform）
  4. 周期性离散信号离散傅里叶变换 (Discrete Fourier Transform, DFT)

FFT4.png

Figure2 四种变换类型

快速傅里叶变换

因为计算机只能处理有限长的离散数值序列，因此适用的是 离散傅里叶变换（DFT）。FFT 是 DFT 的一种快速算法。DFT 运算过程如下：

X(k) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi nk/N} \\ x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) e^{j2\pi nk/N}

其中，

\text{X}(k) — 频域的值
\text{x}(n) — 时域的采样点
\text{n} — 时域采样点的序列索引
\text{k} — 频域值的索引
\text{N} — 采样点数量

对于在毕设中出现的平面曲线的 descriptor，\text{X}(k) 和 \text{Y}(k) 分别代表 \text{x} 轴和 \text{y} 轴的坐标。

S(k) = X(k) + jY(k)
a(u)=\sum_{k=0}^{K-1} s(k)e^{-juk2\pi/K}

其中，

\text{a}(u) — 频域的值
\text{s}(k) — 复数表示的平面坐标点 \text{k}
\text{k} — 点的索引
\text{u} — 频域值的索引
\text{K} — 曲线中的总点数

快在哪？

由公式可以看出，DFT 后的输出包含与输入数量相同的采样点。对于实数信号，FFT 结果具有共轭对称性：

\text{X}(k) = \text{X}^*(\text{N}-\text{k})

因此实际上只有 0 到 \text{N}/2 共 \frac{\text{N}}{2}+1 个点包含有用信息。

FFT 相比 DFT，原本时间复杂度为 O(\text{N}^2) 的运算，变成了 O(\text{N}\log \text{N})。

快速算法实现举例

以下式为例（\text{N}=4）：

X(k) = \sum_{n=0}^{3} x(n) e^{-j2\pi nk/4}

直接计算需要 4 \times 4 = 16 次复数乘法。FFT 采用 Cooley-Tukey 分治策略，将序列按索引的奇偶性分为两组：

X(k) = \underbrace{\sum_{m=0}^{1} x(2m) e^{-j\frac{2\pi}{2}mk}}_{E(k)} + W_4^k \underbrace{\sum_{m=0}^{1} x(2m+1) e^{-j\frac{2\pi}{2}mk}}_{O(k)}

其中 \text{W}_4^\text{k} = e^{-j\frac{2\pi}{4}\text{k}} 为旋转因子。

先计算两个 2 点 DFT：

  • \text{E}(0) = \text{x}(0) + \text{x}(2), \text{E}(1) = \text{x}(0) - \text{x}(2)
  • \text{O}(0) = \text{x}(1) + \text{x}(3), \text{O}(1) = \text{x}(1) - \text{x}(3)

然后通过 蝶形运算 组合得到最终结果：

X(0) = E(0) + W_4^0 \cdot O(0) = E(0) + O(0) \\ X(1) = E(1) + W_4^1 \cdot O(1) = E(1) - j \cdot O(1) \\ X(2) = E(0) + W_4^2 \cdot O(0) = E(0) - O(0) \\ X(3) = E(1) + W_4^3 \cdot O(1) = E(1) + j \cdot O(1)

只需 4 次复数加法 + 2 次复数乘法，远低于直接 DFT 的 16 次复数乘法。当 \text{N} 很大时，这种分治策略将复杂度从 O(\text{N}^2) 降至 O(\text{N}\log \text{N})。

变换前后信号的对应关系

原始信号 幅值 为 \text{A}，FFT 结果的每个点（除第一个直流分量外），幅值就是 \text{A} 的 \frac{\text{N}}{2} 倍。直流（DC）分量，也就是第一个点（0 Hz），它的幅值是直流分量的 \text{N} 倍。

每个点的 相位，就是 该频率下 的信号的相位。

\text{N} 个点均匀分布在 0 到 \text{F}_{\text{s}} 的频率轴上（不含 \text{F}_{\text{s}}），第 \text{k} 个点（从 0 开始计数）的频率为：

f_k = k \cdot \frac{F_s}{N}, \quad k = 0, 1, ..., N-1

在某个点 \text{n} 处所能分辨的 频率 为 \frac{\text{F}_{\text{s}}}{\text{N}}。\text{F}_{\text{s}} 为采样频率，\text{N} 为采样点数。

对于实数信号，由于共轭对称性，有效信息仅在 0 到 \text{F}_{\text{s}}/2（Nyquist 频率） 之间。

FFT3.png

Figure3 FFT 变换后信号的变化

参考文献

[1] FFT的物理意义 available @7/28/2022

[2] 数字图像处理中的傅里叶（DFT/FFT） available @7/28/2022