DFT, FFT与音频处理

鸡哥

计算机图形学话题下的优秀答主

本文原来只是想把知乎当我的云笔记本记录一下FFT的，备着万一以后要搞下图像处理或者图形学里的基于FFT的海洋模拟。但是又想搞点差异化，于是脑子发热想写点音频处理，没想到坑也不浅。且因为我因为本科坑爹的课程安排导致我甚至还没上过信号处理相关的课，所以本文很可能有错误的地方，欢迎轻怼orz，图侵删

傅里叶变换

我们可能或多或少听说过，一个在长度为 P 的定义域里的可积(integrable)实数函数 f(x) 可以展开成无限个三角函数的和，也就是傅里叶级数 [1](Fourier series)：

f(x)= \frac{a_0}{2}+\sum^N_{n=1} \left( a_n\cos \left(\frac{2\pi nx}{P} \right) + b_n \sin \left(\frac{2\pi nx}{P} \right) \right)\\

其中序列 a_n, b_n 就是各个三角函数分量的权重。我们把这一大坨三角函数加权求和就可以逼近我们的原函数 f(x) ，就像泰勒级数那样。
另外呢，傅里叶级数有其他的形式，就是”振幅-相位”(amplitude-phase)形式：

f(x) = \frac{A_0}{2} + \sum^N_{N=1}A_n \cos \left(\frac{2\pi nx}{P} - \phi_n \right)\\

其中振幅序列 A_n= \sqrt{a^2_n +b^2_n} ，相位 \phi_n=\text{atan2}(b_n,a_n) 。这种形式可能是在做数字信号处理里面，我们比较关心的一种。
那么好了，现在我们有了一个已知的函数 f(x) ，要怎么求得这些三角函数的“权重”序列呢？在 N \rightarrow \infty 的时候，我们就用傅里叶变换就好啦：

F(\xi)=\int^\infty_{-\infty}f(x)e^{-2\pi ix\xi} dx\\

不过还要提一下，前面的级数形式看起来还很正常，突然傅里叶变换就在 e 的指数上加了虚数看起来这么恐怖？？其实可以直接用欧拉公式展开，并把傅里叶变换的这个形式和前文提到的其他傅里叶级数的形式联系起来：

e^{ix}=cosx+isinx\\

这些三角函数的“权重”，就是频域里面的系数。我们每确定一次频率参数 \xi 并计算一次傅里叶变换的积分，就可以得到一个“权重”。

离散傅里叶变换(DFT)

DFT与频谱

下一个问题，在计算机的世界里，那必然不可能用无限个权重，输入的函数(信号)也在绝大部分时候不是连续的。于是我们有了 离散傅里叶变换 (DFT, Discrete Fourier Transform)，给定一个离散的实数序列 x(n) ，我们可以用DFT得到一个离散的频谱 (spectrum) X(k) ，其中频谱第 k 个点计算公式为：

X(k)=\sum^{N-1}_{n=0} x(n) e^{-j \frac{2\pi}{N} k}\\

虽然理论上，你输入的是一堆复数(Complex)，输出也是一堆复数。但在做音频处理这种数字信号处理的场景下， 输入是实数信号 ，你要对每个输入信号采样先加上一个 0i ，再拿去做DFT。
那这个DFT得到的频谱要怎么理解/解析/使用呢？音频处理的话，我们肯定是倾向于用 振幅-相位形式 的傅里叶级数来翻译啦。假设我们的频谱的复数值 X(k)=a_k+b_kj ，那么当前频谱点对应的谐波分量的振幅和相位为：

\begin{aligned} \\ A(k)&=\sqrt{a_k^2+b_k^2} \\ \phi(k)&= \text{atan2}(b_k,a_k) \\ \end{aligned} \\

有了上面的公式，我们可以把离散信号DFT的输出结果解析成 振幅谱 和 相位谱 ，分别给出了每个(余弦)三角函数的振幅(权重)及相位。

DFT及频谱的一些性质

1. 实数序列DFT频谱的对称性 ：因为一般情况下我们的输入序列 x(n) 是实数序列，则其频谱 X(k) 满足共轭对称性，即 X(k)=X^*(N-k) 。这带来的一个结果是，振幅谱 A(k) 关于 k=N/2 圆周偶对称，相位谱 \phi(k) 关于 k=N/2 圆周奇对称。这意味着，我们可以做分析的时候会 把频谱砍掉一半 ，一般情况下抛弃后半段。如果砍半之前的振幅谱叫 双边振幅谱 。

2. 折叠频率 ：假设输入信号的采样频率为 f_s ，则这个 k=N/2 处对应的频率为 折叠频率 。其实折叠频率就是 Nyquist Frequency ，其值为 f_s/2 。

3. 频谱点对应的谐波频率 ：振幅谱/相位谱每个点都有其对应的频率， N/2 个点均匀间隔地分布在 [0,f_s/2] 里面，相当于每个点都掌管一小段频率范围，每小段范围大小是 f_s/N 。每个频率小区域我们称为frequency bin。第k个frequency bin的频率 f_k=kf_s/N 。

4. 输入信号补0 ：假设原始输入序列的点数为 N ，在实际操作的时候我们可以通过在后面疯狂补0，把序列扩充到长度为 2N,3N,4N ….，这样子DFT输出的频谱点数会增加， 提高视觉上的分辨率，但遗憾的是频谱分辨能力并不会因此增加 ，补0前区分不开的频率分量，补0后还是区分不开。补0后只能让振幅谱和相位谱平滑连续一点，不需要自己做插值了。

频谱泄漏

理想的连续Fourier Transform的频谱采样精度是无限，在计算机的世界里是没有这等好事的，DFT/FFT由于不能达到无穷的精度，有时每个频谱点对应的频率，不一定就是对应着输入信号的频率。例如我用1024 point FFT分析采样率20480Hz的信号，其源信号Nyquist Frequency是10240Hz，"有效频谱点"数为1024/2=512，则频谱每个点对应20Hz, 40 Hz, 60Hz, … 512 * 20Hz。如果我的输入信号的频率是20Hz的倍数还好，还能完美捕捉。但假设我有个输入信号不是20的倍数赫兹，那么这个信号不能被FFT完美的捕捉到，会有很多能量泄露并散布在整个频谱中，造成 频谱泄漏(spectral leakage) 。
下面是我写的FFT的一些test cases，输入信号为200Hz, 250Hz, 440Hz的正弦波叠加在一起，FFT window =1024 point，输出超过一定强度阈值的频谱点，应该可以体现Spectral Leakage：

图1的frequency bin真是完美避开了源信号的分量频率，导致还挺严重的频谱泄漏。而图2的bin就真的可以全部精确覆盖到输入谐波的频率，就没有严重的泄漏了（但还是有点浮点数误差的）。这些会在主频率隔壁产生很多小山坡，降低频谱分辨率

快速傅里叶变换(FFT)

DIT Base2-FFT

其实按照公式，DFT的实现就两个for几行代码的事（当然你可能得先实现一个复数类 class Complex)：

//baseline O(n^2) DFT
static void DFT(const std::vector<Complex>& inSeq, std::vector<Complex>& outSeq)
    const size_t N = inSeq.size();
    outSeq.reserve(N);
    // for every frequency domain point
    for (uint32_t k = 0; k < N; ++k)
        Complex X_k; // k-th point in freq-domain
        for (uint32_t i = 0; i < N; ++i)