教你理解 FFT 算法及其在音频重采样中的应用
FFT 算法介绍
FFT 算法英文全称是 Fast Fourier Transform，首先我们拆成三个词：Fast（快速的），Fourier（傅里叶，人名），Transform（变换）。傅里叶变换是傅里叶发明的，但是这个变换算法一开始并不是「快速的」。毕竟当年只有数学模型，没有计算机嘛，人们也不会想到去优化，而且也没有人傻到去亲自算加减乘除、泰勒展开等计算。
傅里叶的算法本来是叫「离散傅里叶变换」，Discrete Fourier Transform。离散，意思就是它不是连续的。如何理解离散呢？用程序员的思维去理解的话，你可以这样理解：sin(x) 这是个函数，它作为函数，你输入任意 x，它给你计算正弦值，这是连续的。而反过来，离散是什么意思呢？那就是 float sin_x[360]; 它不是函数，而是数组，你用一个 index 去数组里面取值，你的 index 必须得是整数，对吧？那么它就是不连续的。这就叫离散。
离散傅里叶变换，指的就是把这种「数组」，也就是离散的散点，从时域变为频域（也是离散的散点，而且数量相同），或者反过来从频域再变回时域，这就叫离散傅里叶变换。
离散傅里叶算法代码
// 因为要用到 `size_t` 来表示“大小”，所以包含 `stddef.h`，这是个好习惯来着
// 坏习惯就是使用 `long` 来表示“大小”
#include <stddef.h>
// 要用三角函数就需要用到数学库，所以还需要一个头文件 `math.h`。
// 如果你要用 gcc 编译，你需要链接数学库，编译参数要加上一个 `-lm`
#include <math.h>
// 全程我们使用的浮点数类型是 double，因为 float 的精度太差。
// 圆周率
static const double PI = 3.1415926535897932384626433832795;
// 圆周率×2
static const double PI2 = 6.283185307179586476925286766559;
// 先定义复数
typedef struct Complex_s
{
double re; // 实部
double im; // 虚部
} Complex;
// 离散傅里叶变换算法，输入和输出的长度需要是一致的，必须都是 N。
void DiscreteFourierTransform(size_t N, double *input, Complex *output)
{
for(size_t k = 0; k < N; k ++) // k 是频率
{
Complex result = {0, 0};
for(size_t n = 0; n < N; n++)
{
// 计算相位角 (2πkn/N)
double x = n * k * PI2 / N;
// 累加实部和虚部
result.re += input[n] * cos(x);
result.im += input[n] * sin(x);
}
// 要对输出的结果进行标准化
result.re /= N;
result.im /= N;
output[k] = result;
}
}
// 离散傅里叶变换算法逆变换，输入和输出的长度需要是一致的，必须都是 N。
void InverseDiscreteFourierTransform(size_t N, Complex *input, double *output)
{
for (size_t n = 0; n < N; n++) // 遍历输出序列的每个点
{
double result = 0.0f;
for (size_t k = 0; k < N; k++) // 遍历所有频率分量
{
// 计算相位角 (2πkn/N)
double x = PI2 * k * n / N;
// 计算复指数 e^(j*angle) 的实部和虚部
double cos_val = cos(x);
double sin_val = sin(x);
// 累加复数乘法结果的实部部分
result += input[k].re * cos_val + input[k].im * sin_val;
}
// 存储计算结果（实部）
output[n] = result;
}
}
// 我们需要 `stdio.h` 提供 `printf()` 函数
#include <stdio.h>
// 我们需要 `stdlib.h` 提供 `rand()` 函数用于取得随机数
#include <stdlib.h>
int main()
{
// 离散傅里叶变换并不要求变换长度必须是 2 的 N 次方，随意长度都可，作为例子，我瞎几把定它的长度
#define N 3333
// 我们生成随机数
double data[N];
for(size_t n = 0; n < N; n++)
{
// 直接生成随机数，此时得到的随机数的取值范围是 [0, 1]
data[n] = (double)rand() / RAND_MAX;
// 将随机数的取值范围变为 [-1, 1]，其实也可以不用这么整
data[n] = data[n] * 2 - 1;
}
// 创建复数数组用于接收离散傅里叶算法的结果
Complex dft_result[N];
// 进行变换
DiscreteFourierTransform(N, data, dft_result);
// 创建实数数组用于接收逆变换的结果
double inversed[N];
// 进行逆变换
InverseDiscreteFourierTransform(N, dft_result, inversed);
// 比较逆变换后的结果是不是和原始结果相同（注意浮点数有误差，我们取容差 0.00001）
for (size_t i = 0; i < N; i ++)
{
double diff = fabs(data[i] - inversed[i]);
if (diff > 0.00000000001)
{
printf("在第 %d 个相位点上，算出来的结果的差异值绝对值是 %f\n", (int)i, diff);
}
}
return 0;
}
请在 https://godbolt.org/ 上进行测试，你会发现经过正变换后再经过逆变换，结果无差异。
经过傅里叶变换后，你会得到离散的频域的数组，它是复数。这个复数是干嘛的呢？请这样理解：
- 复数的数组的
index 是正弦函数的频率值（相对于整段波形的长度的频率值）。
- 复数的实部和虚部其实可以被当作一个二维向量的
x 和 y。
- 这个向量与坐标轴的夹角可以理解为正弦函数的相位。
- 这个向量的模可以理解为正弦函数的振幅。
要想将这个频域的复数数组转换为时域的实数数组，只需这么做：
- 以复数的
index 为频率值。
- 以复数的实部和虚部构成二维向量
(x, y)。
- 以二维向量的模为振幅。
- 以二维向量到坐标轴的夹角为相位。
- 有了频率，振幅，有了相位，你就可以弄出对应的正弦波。
- 把所有的正弦波全部叠加到一起，就能还原出之前的
input 里面的波形。
不过实际上没必要傻乎乎地去真的去算什么夹角（用 atan2()，但是这个算法弔慢，即便它是 FPU 支持的功能），看上面代码提供的逆变换就根本没有在使用什么夹角——这正是复数运算的精髓所在：直接在笛卡尔坐标系完成所有向量操作，这也是为啥 DFT、FFT 计算需要用到复数的原因。
离散傅里叶算法的来历
法国人傅里叶大佬 经常生病但是当时医疗条件不足，他就不知道从哪里听来一个说法， 认为「加热人体可以帮助人体治疗一切疾病」。在研究热传导的时候，他于 1805 年提出一个算法，如果不是连续的波形而是这种波形样本的数列，他可以根据相同的数列长度的个数的正弦波或者余弦波以不同的频率、振幅、相位来构成完全相同的数列。他发稿的时候直接被退回，被认为是一派胡言，骂他的人还包括拉普拉斯，说：「傅立叶先生，您对数学一无所知」。
傅里叶努力了 15 年进行抗争，抗争的结果是 1822 年出版了《热的解析理论》。在 1830 年，傅里叶又生病了，他又要进行加热疗法来治疗自己，在一个暑热的夏天，他在家里烧火炉，裹着厚厚的被子，浑身出汗。最终死于中暑。哎。
到了 1965 年，有 Cooley 大佬和 Tukey 大佬研究了傅里叶变换，发明了素因子分解法（其实就是现在的蝶形 FFT 算法的雏形）也还是发稿的时候惨遭退回，原因是蝶形算法只适用于 2ⁿ，被认为局限性太大，没啥卵用。
再到 1973 年，Tom Stockham 和 Charile Rader 画了一个蝶形图，终于大家看懂了这是啥玩意儿了，感叹道：哇！这 FFT 的加速效果怎么这么明显，一个本来时间复杂度弔贵的算法被简化成几十倍的高效率，那是相当强大的发明。终于大家认可了 FFT 变换算法的价值。
快速傅里叶算法，真的只能处理 2ⁿ 大小的数据么？
首先我们要搞清楚快速傅里叶算法的本质：那就是快速的离散傅里叶算法。快速是什么意思？意思就是快。我举个最简单的例子：
void DiscreteFourierTransform(size_t N, double *input, Complex *output)
{
// 莽夫式快速傅里叶算法：用多线程暴力硬怼原始DFT
#pragma omp parallel for
for(size_t k = 0; k < N; k ++) // k 是频率
{
Complex result = {0, 0};
for(size_t n = 0; n < N; n++)
{
// 计算相位角 (2πkn/N)
double x = n * k * PI2 / N;
// 累加实部和虚部
result.re += input[n] * cos(x);
result.im += input[n] * sin(x);
}
// 要对输出的结果进行标准化
result.re /= N;
result.im /= N;
output[k] = result;
}
}
看，我只需要启用 OpenMP，然后在主循环那里加一句 #pragma omp parallel for，使其从单线程计算变为多线程计算，它是不是就是变快了呢？只要在多核的处理器上这么一整，所有的 CPU 都参与进来进行离散傅里叶变换。它是不是会变快？我这是不是就是快速傅里叶变换？你能反驳吗？
然后，当你承认我这个算法是快速傅里叶变换后，回到之前的问题：快速傅里叶算法，真的只能处理 2ⁿ 大小的数据么？显然不是的。其实，程序员有程序员的办法，数学家有数学家的办法。蝶形算法只是快速傅里叶变换的算法的一种（利用了复数相乘实现复平面上的相位旋转），其它的还有：
- 素因子分解算法（Prime Factor Algorithm）
- 通过互质分解避免使用旋转因子，直接组合子DFT结果。
- 布鲁斯坦FFT（Bluestein FFT）
- 对输入的散点，用上卷积乘法。
- 卷积核，其实就是预先算好的三角函数值序列。通过预先计算三角函数来降低对三角函数的使用次数。
- 可以处理任意长度。
- 雷德算法（Rader's FFT）
- 温诺格拉德傅里叶变换算法（Winograd FFT）
- 混合基FFT（Mixed-Radix FFT）
有的算法通过减少重复的三角函数计算来加快 DFT 变换，有的算法通过加法换乘法，有的算法减少复数相乘的次数，八仙过海各显神通。
当然实际上我那个使用 OpenMP 计算的算法并不是合格的 FFT 变换——它并没有降低计算的时间复杂度，而是靠堆砌硬件算力来莽夫式解决问题。
更多的算法这些我就不一一列举了。当你使用开源的成熟的 FFT 库的时候，这些库会根据你需要的处理长度来灵活选择最合适的算法。
使用快速傅里叶变换进行音频重采样
就目前为止我给的例子都是 C 语言写的，但是现在我可能需要提供一些 Rust 的代码——Rust 的 rustfft 库真的很好用，而我们想要研究的重点，是如何利用好 FFT 变换，对音频进行重采样，并避免在重采样过程中造成音频的杂波或音质的降低。我会把我的 Rust 的代码写得尽可能像 C/C++。
音频重采样处理的思路
重采样的本质就是拉长或压缩波形（的数组），这个过程中涉及到插值算法，直接在时域进行插值的话，无论你使用什么插值算法，都会引入杂波。最合理的方式当然就是要先用 FFT 将波形转换为频域，然后在频域里进行拉伸、压缩处理，在频域里进行插值，最后再转换为时域的波形信号。这样你就可以精准控制你想要操作的声波的频率了。
在 FFT 变换后的频域的数据里，我们可以把各种各样的插值算法都拿来用，就包括图像的拉伸和压缩一样，甚至最简单的临近插值都是可以使用的。做好插值后，再转换回时域的信号后，你会发现不会引入杂波。在频域里使用更好的插值算法只是锦上添花，当然使用临近插值也太敷衍了。我使用线性插值，我觉得足够了。
但是，FFT 变换不会改变散点的数量，而我们要做的重采样，是会改变散点数量的。所以不仅要进行 FFT 变换，在变换前和变换后，我们还有一些东西需要处理。
音频重采样第一步：将音频剪切成一段一段的来处理
直接对一整段音频进行 FFT 变换的话，除非你用蝶形算法，否则时间复杂度会变得极高，而使用蝶形算法的话，最差情况下你需要处理几乎两倍于原始音频长度的 FFT 变换，即使时间复杂度比别的 FFT 变换算法的时间复杂度低，但空间的使用率会变得很高。
所以要将音频分段处理。分段处理还有一个好处就是：对于可以被视作无限长度的音频，或者即时录音的音频流的重采样，你也可以处理。而分段的长度需要合理计算。人耳听觉的最低到最高的频率范围是 20 Hz 到 20 kHz。低于 20 Hz 的音频，我们叫它「次声波」。
当我们把音频按 每秒 20 段 的长度进行分割的话，即使我们的重采样算法并非完美，引入了次声波，那 人类也听不到。除非你的算法需要满足犬类、猫类、蟑螂类客户的需求，届时需要重新调整分段的长度。
另外，进行分割时，我们需要确保分段的 长度能够被 2 整除。因为在进行 FFT 变换时，我们的音频数据需要转换为复数。虽然上述我提供的 DFT 算法输入的音频信号是实数数组，但是流行的 FFT 算法库都是 输入和输出都是复数 的。
第二步：取一个合适的 FFT 变换的长度
我们的目标是要在频域范围内进行频率的缩放，当我们 拉长音频 的时候，在频域我们需要 缩小插值 散点，使高频的波形往低频跑；而当我们需要 缩短音频 的时候，在频域内我们需要 放大插值，使低频的波形往高频跑。
FFT 处理的长度至关重要：
- 它影响 FFT 库采用的 FFT 算法的类型。如果你用的长度是 2ⁿ 长度，它就能给你选择时间复杂度最低的蝶形算法；而如果你用的是质数的长度，它可以给你应用雷德算法。不同的算法有不同的时间复杂度和空间复杂度。
- 它影响精度。FFT 处理的点位数越多，精度越高，处理下来的音质越好，时间复杂度越高，但是我们正好可以往 2ⁿ 的长度上去取，这样正好可以利用时间复杂度最低的蝶形算法。
- 它影响你拉长音频（降低频率）的时候，最长你能拉多长。比如我的 FFT 处理大小是
32768，那么不论输入的音频是多少，最长也只能拉长到 32768 个样本。
我推荐的处理长度是你的音频的采样率向上取一个 2ⁿ 长度，比如你的采样率是 44100 Hz，那么 FFT 处理长度就取 65536。与此同时，按照第一步，将音频进行分段的过程，我们每次处理的一段音频的长度是 2205 个样本。这样下来，我们能把音频拉长的最大倍数约等于 14.86，已经非常绰绰有余了。
第三步：将每段音频转换为 FFT 库所需的复数
我们继续按 44100 Hz 的采样率为例，我们需要把 2205 个音频样本放入到长度为 65536 的复数数组里。这里有两个方案：
- 我们把样本按对称共轭的方式填充 FFT 数组的前
1103 和后 1103 个复数里（其实最靠近中间的两个复数只填充了实部，因为音频样本用完了，实际上我们填充的是前后各 1102.5 个复数），同时填充实部和虚部，其余复数全部置零。FFT 变换后，将复数数组的前 2205 个复数的实部和虚部都取出来作为我们的时域样本。这是精度损耗最小的方式。
- 填充方式：
void RealToComplex(float *samples, size_t num_samples, Complex *fft_buffer, size_t fft_len)
{
size_t half = fft_len / 2;
size_t back = fft_len - 1;
size_t j = 0; // j 用来从 samples 里取样本。
// 先把复数数组清零。
memset(fft_buffer, 0, fft_len * sizeof fft_buffer[0]);
for (size_t i = 0; i < half; i ++)
{
if (j + 1 < num_samples)
{
// 前部
fft_buffer[i].re = samples[j + 0];
fft_buffer[i].im = samples[j + 1];
// 后部
fft_buffer[back - i].re = samples[j + 0];
fft_buffer[back - i].im = -samples[j + 1]; // 共轭复数
}
else if (j + 1 == num_samples) // 出现这种情况的时候，样本是奇数个，现在是最后一个样本，填写到实部
{
// 前部
fft_buffer[i].re = samples[j];
fft_buffer[i].im = 0;
// 后部
fft_buffer[back - i].re = samples[j];
fft_buffer[back - i].im = 0;
}
else // 样本全部取完。
{
// 其余复数使其实部虚部为零
break;
}
j += 2;
}
// 如果运行到这里，j 大于 num_samples，说明你给的样本数太多了，多出来的部分全部被丢弃。那是错误。
assert(j <= num_samples);
}
- 逆变换后，取回样本的方式：
void ComplexToReal(float *samples, size_t num_samples, Complex *fft_buffer, size_t fft_len)
{
size_t half = fft_len / 2;
size_t back = fft_len - 1;
size_t j = 0;
// 复数数组后半部分数据全部丢弃。
for (size_t i = 0; i < num_samples; i ++)
{
// 取回我的样本。
if ((i & 1) == 0)
num_samples[i] = fft_buffer[j].re;
else
num_samples[i] = fft_buffer[j++].im;
}
}
- 另一种方式：只把所有实数填写到 FFT 处理的复数的实数部分，虚部留空。这种方式会对音质造成损耗，样本的剪切边界会有咔哒声。
- 填充方式：
void RealToComplex(float *samples, size_t num_samples, Complex *fft_buffer, size_t fft_len)
{
assert(num_samples <= fft_len);
// 先把复数数组清零。
memset(fft_buffer, 0, fft_len * sizeof fft_buffer[0]);
for (size_t i = 0; i < num_samples; i++)
{
// 只填实部
fft_buffer[i].re = samples[i];
}
}
- 逆变换后，取回样本的方式：
void ComplexToReal(float *samples, size_t num_samples, Complex *fft_buffer, size_t fft_len)
{
for (size_t i = 0; i < num_samples; i++)
{
// 丢弃虚部
samples[i] = fft_buffer[i].re;
}
}
第四步：进行 FFT 变换
不论是使用 DFT 变换，还是莽夫式多线程 DFT 再加 SIMD 的变换，还是预计算三角函数然后搞卷积，还是蝴蝶算法，总之结果都一样：我们得到了频域的数据，它是对称的（至于什么叫「对称」，一会儿你就知道了）。
但按照之前我们取用的 FFT 大小，它是 2ⁿ，正好最适合使用蝴蝶算法，时间复杂度最低。
你也许应该注意到了，我没有把所有的音频样本全部都填充到 FFT 处理的缓冲区里面去。因为本来重采样就不是这么使用 FFT 的。
| 实0 |
虚0 |
实1 |
虚1 |
实2 |
虚2 |
实3 |
虚3 |
..略.. |
实3 |
虚3 |
实2 |
虚2 |
实1 |
虚1 |
实0 |
虚0 |
观察上面的表格，每个复数由一个实部和虚部构成，表格里写的数字，其实是这个复数频域点的频率值。
而且，这个复数的数组，它是轴对称的，最左边和最右边的都是最低频的部分，越靠近中间，频率值越高。
我们把音频样本填入复数数组的时候，应当采用第一种方案。我们也是对称填入进去的。两边的复数共为共轭复数。
第五步：频域内进行插值
按照之前说的，如果要拉长波形，那就要压缩频域点；如果要压缩波形，那就要拉长频域点。
根据第四步的结果，你可以看到，FFT 输出的复数是对称的。我们做插值的时候，也是要进行对称处理。
- 如果要压缩频域点，左半边的频域点往左半边压缩；右半边的频域点往右半边压缩。
- 如果要拉长频域点，左半边的频域点要往中间放大；右半边的频域点也要往中间放大。
这里说的「压缩」「拉长」「缩小」「放大」都是线性缩放。每个频域点的缩放比都是相同的。就像你用 Photoshop 拉伸、压缩图像一样。
直接对复数进行线性插值即可。我们使用的 FFT 变换的缓冲区很大，有 65536 个点，变换精度非常高，那么采样算法就可以不用太复杂，线性插值刚刚好。
缩放比 = 目标采样率 ÷ 原始采样率。根据缩放比，对时域点进行缩放。
第六步：进行逆变换，得到时域点（长度等于 FFT 处理的长度），裁剪为我们需要的目标长度
按照上述例子，我们取了 65536 这个 FFT 长度，但是我们给的音频长度很短，就是 20ms 的一段（按照 44100 Hz 采样率，我们给的音频长度是 2205 个时域点）
那假设，我们的重采样是想把 44100 Hz 的音频变为 48000 Hz 的音频，缩放比约为 1.088435 ，实际重采样处理完后，我们得到了 65536 个时域点，我们需要按照缩放比对这么多的时域点（绝大多数都是多余的）进行裁切。
(48000 / 44100) * 2205 = 2400
从这 65536 个时域点里，取 2400 个时域点即可。这 2400 个时域点，就是我们经过重采样后，把 44100 Hz 的音频的 2205 个时域点的波形按照比例拉长后得到的 48000 Hz 频率的 2400 个时域点。就这样，将音频分成一段一段处理，最后再拼接起来，就完成了一整个音频的重采样处理。
这大家也都看到了，原理就是这样的，需要我提供代码吗？
let mut fftbuf: Vec<Complex<f64>> = Self::real_to_complex(samples);
if fftbuf.len() <= self.fft_size {
fftbuf.resize(self.fft_size, Complex{re: 0.0, im: 0.0});
} else {
return Err(ResamplerError::SizeError(format!("The input size {} must not exceed the FFT size {}", fftbuf.len(), self.fft_size)));
}
// 进行 FFT 正向变换
self.fft_forward.process(&mut fftbuf);
// 准备进行插值
let mut fftdst = vec![Complex::<f64>{re: 0.0, im: 0.0}; self.fft_size];
let half = self.fft_size / 2;
let back = self.fft_size - 1;
let scaling = desired_length as f64 / input_size as f64;
if input_size > desired_length {
// 输入大小超过目标大小。意味着音频的压缩。
// 在时域进行时域点的拉伸。
for i in 0..half {
let scaled = i as f64 * scaling;
let i1 = scaled.trunc() as usize;
let i2 = i1 + 1;
let s = scaled.fract();
if INTERPOLATE_DNSCALE { // 做线性插值
fftdst[i] = interpolate(fftbuf[i1], fftbuf[i2], s);
fftdst[back - i] = interpolate(fftbuf[back - i1], fftbuf[back - i2], s);
} else { // 做临近插值
fftdst[i] = fftbuf[i1];
fftdst[back - i] = fftbuf[back - i1];
}
}
} else {
// 输入大小小于目标大小，意味着音频的拉长。
// 在时域进行时域点的压缩。
for i in 0..half {
let i1 = (i as f64 * scaling).trunc() as usize;
let i2 = ((i + 1) as f64 * scaling).trunc() as usize;
if i2 >= half {break;}
let j1 = back - i2;
let j2 = back - i1;
if INTERPOLATE_UPSCALE { // 做平均插值
fftdst[i] = get_average(&fftbuf[i1..i2]);
fftdst[back - i] = get_average(&fftbuf[j1..j2]);
} else { // 做临近插值
fftdst[i] = fftbuf[i1];
fftdst[back - i] = fftbuf[back - i1];
}
}
}
// 进行逆变换
self.fft_inverse.process(&mut fftdst);
// 得到时域点
let mut real_ret = Self::complex_to_real(&fftdst);
// 剪切时域点，使其长度达到我们想要的长度
real_ret.truncate(desired_length);
结论
- 所谓快速傅里叶变换算法其实就是各种各样的针对离散傅里叶变换算法的优化版。
- 快速傅里叶变换算法里，只有蝴蝶算法要求样本点数为 2ⁿ 大小，其它算法则有其它算法的大小要求，也有的算法对处理大小无要求。
- 进行重采样计算的时候，FFT 的处理大小和你的音频的分段大小并不绑定。只是，FFT 的处理大小越大，你的音质损耗越低。就比如，我用
65536 的处理大小去处理音频，但我把音频分割成每 20ms 一段（按照 44100 Hz 采样率，20ms 相当于 2205 个样本）来进行处理（只给 FFT 处理器输入这么长的样本）。处理完后得到的音频长度是 65536，但只有最前面那一段是我需要的。我通过直接裁切它，并最终把这些裁切出来的一段段的音频拼接成完整音频的方式来获取我想要的结果。
参考源码
|Archiver|小黑屋|技术宅的结界
( 滇ICP备16008837号 )|网站地图
发表于 
发表于 
发表于