快速傅里叶变换

// 来源: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;
}
https://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B
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)
输出结果。正逆变换一轮之后 ifft 接近但不等于输入的数据,原因是浮点运算产生的误差

注意到测试用例的输入采样数据是一个正弦曲线,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++ 里了,特别是各种处理数组的算法。但区别于其他容器,它的元素类型必须满足数值类型的具名要求,如整数、浮点数、复数等等。