🎈前言

傅立叶变换（Fourier Transform，FT）可以将时域上的信号转变为频域上的信号，使用者从而获得感兴趣的变量，在音视频、图像、神经影像等信号处理领域均有广泛的应用。而快速傅立叶变换（Fast Fourier Transform，FFT）是一种离散信号的傅立叶变换的的快速计算方法。

在上篇博客我复习了一下傅立叶变换的数学推导，这里就回顾一下快速傅立叶变换的实现。

✏算法原理

这里记原信号为 $x[n]$，傅立叶变换的结果为 $x[k]$，信号的周期为 $T$，在一个周期内采样 $N$ 次，离散傅立叶变换（DFT）为

$$X[k]=\sum_{k=0}^{N -1}x[n]e^{-2\pi in k/N}$$

直接计算 DFT 需要 $N^2$ 次乘法。FFT 利用 $e^{i\theta}$ 的周期性和分治的思想减少了需要的计算。

当采样信号个数 $N$ 为偶数时，有

$$X[k]=\sum_{n=0}^{N -1}x[n]e^{-2\pi in k/N}\\ =\sum_{n=0}^{N/2 -1}x[2n]e^{-2\pi i\frac{2n k}{N}}+\sum_{n=0}^{N/2 -1}x[2n+1]e^{-2\pi i\frac{ (2n+1) k}{N}}\\ =\sum_{n=0}^{N/2 -1}x[2n]e^{-2\pi i\frac{2n k}{N}}+e^{-2\pi i\frac{ k}{N}}\sum_{n=0}^{N/2 -1}x[2n+1]e^{-2\pi i\frac{2nk}{N}}$$

记

$$E[k]=\sum_{n=0}^{N/2 -1}x[2n]e^{-2\pi i\frac{2n k}{N}}\\O[k]=\sum_{n=0}^{N/2 -1}x[2n+1]e^{-2\pi i\frac{2nk}{N}}\\W[k]=e^{-2\pi i\frac{k}{N}}\\ X[k]=E[k]+W[k]O[k]$$

注意到 $e^{i\theta}=e^{i(\theta+2\pi)}$，$e^{i\theta}=-e^{i(\theta+\pi)}$，有

$$E[k] = E[k+N/2]\\ O[k] = O[k+N/2]\\ W[k] = -W[k+N/2]\\ X[k+N/2]=E[k]-W[k]O[k]$$

这样子，这里只需计算一半的 $X[k]$。进一步地，可以发现 $E[k]$ 和 $O[k]$ 中也可以利用 $e^{i\theta}$ 的周期性减少计算，只要 $N$ 是 2 的正整数次幂（不是的话后面补 0），就可以把 $E[k]$ 和 $O[k]$ 的计算细分为更小的子问题，这就是所谓的蝶形变换。

$$X_0[k]=X[k],k=0,1,...N-1\\$$

$$X_1[k]=E[k]\\X_1[k+N/2]=O[k]\\ X_0[k] = \left\{ \begin{array}{lr} X_1[k]+W[k]X_1[k+N/2],k<N/2\\ X_1[k-N/2]+W[k]X_1[k],else \end{array} \right. \\ k=0,1,...N/2-1\\ ...$$

同样地，$X_1$（如果有的话）也可以根据 $e^{i\theta}$ 周期性分解为更小的子问题：

$$X_{1}[k]=E_{1}[k]+W^2[k]O_{1}[k]\\ X_{1}[k+N/4]=E_{1}[k]-W^2[k]O_{1}[k]\\$$

$N=4$ 时，如下图所示，

（开始胡言乱语…）以此类推，可以注意到以下规律：

$$X_d[k] = \left\{ \begin{array}{lr} X_{d+1}[k] +W^{d+1}[k]X_{d+1}[k+(N/2^d)],\\ \text{ if }!(1\&[2^{d+1}k/N])\\ \text{ }\\ X_{d+1}[k-(N/2^d)] +W^{d+1}[k]X_{d+1}[k],\\\text{ else} \end{array} \right.\\ d=0,1,2...log_2N$$

此时，计算 FFT 只需要 $O(n\log{n})$ 的时间复杂度。

再观察最左边的步骤，$X_{log_2N}$ 所需计算的 $x$ 的下标为 $reverse(N)$，$reverse$ 表示位数为 $log_2N$ 时，1 到 N 的二进制翻转。

对于 $N=4$，$reverse(N)$ 是[0,1,2,3]的二进制表示[00,01,10,11]的左右翻转[00,10,01,11]，即为[0,2,1,3]。

现在已经知道初始状态和状态转移方程，在分治算法中直接向上归并即可。

🎁Code

这里直接给出代码，（虽然刚学的 rust，写起来一点都不优雅）

rust
use num_complex::Complex;
use std::f64::{self, consts::{E, PI}};
fn binary_reversed(n: usize) -> Vec<usize> {
    (0..n).map(|t: usize| {
        let mut res: usize = 0;
        let mut x = t;
        for _i in 0..(n as f64).log(2.0) as usize {
            res = (res << 1) | (x & 1);
            x >>= 1;
        }
        res
    }).collect()
}
fn w(n: f64, k: f64, count: f64) -> Complex<f64> {
    Complex::new(E, 0.0).powc(Complex::new(-2.0_f64.powf(n) * PI * k / count, 0.0) * Complex::I)
}
fn fft(mut x: Vec<Complex<f64>>) -> Vec<Complex<f64>> {
    let len_of_x: usize = x.len();
    let depth = (len_of_x as f64).log2().ceil();
    let expected_len = (2.0_f64).powf(depth);
    for _n in 0..(expected_len - len_of_x as f64) as usize {
        x.push(Complex::new(0.0, 0.0))
    }
    let expected_len_usize = expected_len as usize;
    let mut pre: Vec<Complex<f64>>;
    let mut ans: Vec<Complex<f64>> = vec![];
    let inverted_indexex: Vec<usize> = binary_reversed(expected_len_usize);
    for i in 0..expected_len_usize {
        ans.push(x[inverted_indexex[i]]);
    }
    pre = ans.clone();
    let mut temp = 2.0_f64.powf(depth);
    for d in 0..depth as i32 {
        for k in 0..expected_len_usize {
            let mut _o: Complex<f64>;
            let mut _e: Complex<f64>;
            let _w = w(depth - d as f64, k as f64, expected_len);
            if 1 & (temp as usize * k / expected_len_usize) == 1 {
                _e = pre[k - (expected_len / temp) as usize];
                _o = pre[k];
            } else {
                _o = pre[k + (expected_len / temp) as usize];
                _e = pre[k];
            }
            ans[k] = _e + _o * _w;
        }
        pre = ans.clone();
        temp = temp / 2.0; 
    }
    ans
}
fn main() {
    let data: Vec<_> = vec![
        Complex::new(4.0, 0.0),
        Complex::new(3.0, 0.0),
        Complex::new(2.0, 0.0),
        Complex::new(1.0, 0.0),
    ];
    let ans = fft(data);
    for complex_num in ans.iter() {
        println!("{:?}", complex_num);
    };
}
// Complex { re: 10.0, im: 2.4492127076447545e-16 }
// Complex { re: 1.9999999999999998, im: -2.0000000000000004 }
// Complex { re: 2.0, im: -7.347638122934264e-16 }
// Complex { re: 2.0, im: 1.9999999999999998 }

🗑结语

终于写完这个博客了。这个博客简析了 FFT 的数学推导，并且实现了一个简单的 FFT 程序，虽然感觉自己单纯在写着玩，对工作什么的完全没有任何帮助就是了…