快速傅里叶变换
在 Rosetta Code 上看到的快速傅里叶变换算法，其中用到了标准库中的 complex 和 valarray。我第一次注意到有 valarray 这个库，代码看起来很有函数式处理数组的味道，刚好很适合用在蝶形结递归里。代码简洁，并且顺便代入一个正弦波采样试了试。
// 来源：https://rosettacode.org/wiki/Fast_Fourier_transform
// 有少量修改
#include <complex>
#include <valarray>
#include <iostream>
const double PI = std::acos(-1);
typedef std::complex<double> complex_t;
typedef std::valarray<complex_t> complex_array;
void fft(complex_array& x) {
const size_t N = x.size();
if (N <= 1) return;
complex_array even = x[std::slice(0, N/2, 2)];
complex_array odd = x[std::slice(1, N/2, 2)];
fft(even); fft(odd);
for (size_t k = 0; k < N/2; ++k) {
complex_t t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
void ifft(complex_array& x) {
x = x.apply(std::conj);
fft( x );
x = x.apply(std::conj);
x /= x.size();
}
int main() {
const complex_t test[] = {
1.0, 1.7, 2.0, 1.7, 1.0, 0.3, 0.0, 0.3,
1.0, 1.7, 2.0, 1.7, 1.0, 0.3, 0.0, 0.3,
};
complex_array data(test, 16);
fft(data);
std::cout << "fft" << std::endl;
for (auto i: data)
std::cout << i << std::endl;
ifft(data);
std::cout << std::endl << "ifft" << std::endl;
for (auto i: data)
std::cout << i << std::endl;
}
fft
(16,0)
(0,0)
(6.89019e-16,-7.9598)
(0,0)
(0,0)
(0,0)
(1.9916e-16,0.040202)
(0,0)
(0,0)
(0,0)
(-1.9916e-16,-0.040202)
(0,0)
(0,0)
(0,0)
(-6.89019e-16,7.9598)
(0,0)
ifft
(1,-0)
(1.7,-5.55112e-17)
(2,1.2326e-32)
(1.7,-0)
(1,-0)
(0.3,5.55112e-17)
(0,-1.2326e-32)
(0.3,-0)
(1,-0)
(1.7,-5.55112e-17)
(2,1.2326e-32)
(1.7,-0)
(1,0)
(0.3,5.55112e-17)
(0,-1.2326e-32)
(0.3,-0)注意到测试用例的输入采样数据是一个正弦曲线，N=16，
可以看到整个采样内出现了恰好两个周期。不妨视这段采样恰为一段单位时间，那么上图可以将其视为频率为 2Hz 的信号和一个 0Hz 直流信号的叠加。如果用我个人觉得比较直观的记法，也就是说，$$f=1 \sin(t, 0\mathrm{Hz}) + 1 \sin(t, 2\mathrm{Hz}).$$
接下来试着看着 FFT 的输出结果反推这个频率构成。根据傅里叶变换的对称性，原函数如果全为实数（“电平”输入的16个点显然是实数），那么变换后的函数恰好是左右对称的。所以只需要看前一半8个点。需要注意最开始的第零个点被重复叠加了一倍。
回过头看上面的输出结果，在这里 FFT 输出的前8个点（为了简洁起见略去了那些极小的量）可以视为：$$16, 0, -7.9598i, 0, 0, 0, 0,0$$
第0个点代表了 0Hz 上的信号（的两倍），而第2个点代表了 0Hz 上的信号。取这些复数的模即可得到该频率上的振幅，取相角即可得到该频率的相位。但取模之后应该除以 N/2 才能得到真实幅值：
$$
\begin{align}
0\mathrm{Hz}&: & 16\div 2&=8,\\
& & 8\div \frac{N}{2} &=1.\\
2\mathrm{Hz}&: & 7.9598\div \frac{N}{2}&=1.
\end{align}
$$
代码里出现的标准库以及 cppreference 链接：
- complex：复数库
- valarray：数值数组库，函数式数组操作框架
数值数组库便于递归地处理蝶形结网络，弥补了类 C 数组里难以处理不连续切片的缺陷，在这个案例里 Cooley-Tukey FFT 算法里的蝶形结网络的分治是奇偶交错的。gcc 默认命令不会开启优化，开启优化编译以达到较好的性能。
第一次见 valarray 这个库文件。感到有些惊喜，这么一来好用的函数式写法可以塞进 C++ 里了，特别是各种处理数组的算法。但区别于其他容器，它的元素类型必须满足数值类型的具名要求，如整数、浮点数、复数等等。