数值计算领域不会梦到快速傅里叶算法

-迪蜡熊如是说(雾)

---

傅里叶变换回顾·

无论是在信号分析还是谱分析, 傅里叶变换总是一个强大的工具, 把时域转换到频域上进行分析. 而在数值计算(计算物理)领域我们则需要把傅里叶变换到计算机上去实现.

对于周期函数 $f(x)$, 其周期为 $2l$, 则其可以写成傅里叶级数

$$\begin{equation*} f(x) = a_0 + \sum_{k}^{\infty} \left( a_{k} \cos \frac{k\pi}{l} x + b_{k} \sin \frac{k\pi}{l} x \right) \end{equation*}$$

其中各个系数为

$$\begin{equation*} a_0 = \frac{1}{2l}\int_{-l}^{l} f(\xi) \mathrm{d}\xi \qquad a_{k} = \frac{1}{l}\int_{-l}^{l} f(\xi) \cos \frac{k\pi\xi}{l} \mathrm{d}\xi \qquad b_{k} = \frac{1}{l}\int_{-l}^{l} f(\xi) \sin \frac{k\pi\xi}{l} \mathrm{d}\xi \end{equation*}$$

而对于非周期函数, 无法采用傅里叶级数, 则使用 $l \to \infty$ 的版本,也就是傅里叶积分.

$$\begin{equation*} f(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} g(\omega) e^{-i \omega x} \mathrm{d}\omega \end{equation*}$$

其中傅里叶系数

$$\begin{equation*} g(\omega) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} f(x) e^{i \omega x} \mathrm{d}x \end{equation*}$$

在计算机实现中, 我们会利用欧拉公式把三角函数全部转化成 $e$ 指数. 并记 $w = e^{- 2\pi i /N}$.

离散傅里叶变换(DFT)·

计算机上的傅里叶变换是离散的, 即只能处理有限个点的函数. 此时的离散傅里叶变换其实就是一个插值问题, 用三角函数作为基函数展开函数.

离散傅里叶变换是对一些点组成的点集进行傅里叶变换. 比如等距分布在 $[0,l]$ 中的 $N$ 个点, 并假设在 $[0,l]$ 外这些点是周期的. 此时这个点集的DFT可以写成

$$f_k = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} g_n e^{i 2\pi k n / N} = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} g_n w^{-kn}$$

而傅里叶系数由

$$g_n = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} f_k e^{-i 2\pi k n / N} = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} f_k w^{kn}$$

给出. 注意到 $g_{n}$ 可能为复数, 而 $g_{k}$ 和 $g_{n-k}$ 互为复共轭, 于是 $f(x)$ 能始终保证为实数.

为了确定一个点集的傅里叶系数, 我们需要求解线性方程组, 比如对于一个 $N=4$ 的点集, 我们有

$$T \begin{bmatrix} g_0 \\ g_1 \\ g_2 \\ g_3\end{bmatrix} = \frac{1}{N} \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & w^{-1} & w^{-2} & w^{-3} \\ 1 & w^{-2} & w^{-4} & w^{-6} \\ 1 & w^{-3} & w^{-6} & w^{-9} \\\end{bmatrix} \begin{bmatrix} g_0 \\ g_1 \\ g_2 \\ g_3 \end{bmatrix} = \begin{bmatrix} f_0 \\ f_1 \\ f_2 \\ f_3 \end{bmatrix}$$

求一个矩阵的逆是 $O(n^{3})$ 的, 但是这个矩阵的逆是很显然的

$$T \cdot \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & w^{1} & w^{2} & w^{3} \\ 1 & w^{2} & w^{4} & w^{6} \\ 1 & w^{3} & w^{6} & w^{9} \\\end{bmatrix} = I$$

于是原本需要矩阵求逆的操作变成了矩阵乘法, 时间复杂度从 $O(n^{3})$ 降到 $O(n^{2})$.

这里提供一个简单的DFT实现，利用的是julia语言：

function dft(f)
    N = length(f)
    w = exp(-2im * pi / N)
    g = zeros(Complex{Float64}, N)
    for n = 1 : N
        for k = 1 : N
            g[n] += f[k] * w^((n-1) * (k-1))
        end
    end
    return g
end

快速傅里叶变换(FFT)·

快速傅里叶变换的思路很简单: 分治. 因为可以观察到DFT的展开式可以被分为奇项和偶项:

$$g_{j} = \sum_{k=0}^{N-1} f_{k} w^{-jk} = \sum_{k=0}^{N/2-1} f_{2k} e^{-i 2\pi j (2k) / N} + \sum_{k=0}^{N/2-1} f_{2k+1} e^{-i 2\pi j (2k+1) / N} = x_{j} + y_{j} w^{j}$$

其中

$$x_{j} = \sum_{k=0}^{N/2-1} f_{2k} ~e^{-i 2\pi j k / (N/2)} \qquad y_{j} = \sum_{k=0}^{N/2-1} f_{2k+1}~ e^{-i 2\pi j k / (N/2)}$$

这里我们忽略了因子 $1/\sqrt{N}$, 这个因子在之后手动加上即可. 于是我们可以把一个 $N$ 阶的DFT分解为两个 $N/2$ 阶的DFT. 这样我们就可以递归地分解到 $N=1$ 的情况, 这时候 $g_{0} = f_{0}$.

同时, $g_{j}$ 和 $g_{j + N /2}$ 之间还存在一重对称性(当 $N$ 是偶数时), 即

$$g_{j} = x_{j} + y_{j} w^{j} \qquad g_{j + N/2} = x_{j} - y_{j} w^{j}$$

这样可以减小我们一半的计算量.

现在来分析FFT算法的复杂度: 假设 $N = 2^{M}$, 在每一级需要的操作是 $O(N)$ 的, 则可以得到递归式

$$T(N) = 2T(N / 2) + O(N)$$

利用主定理(master theorem), 可以直接得到 $T(N) = O(N \log N)$.

---

下面是 $N = 4$ 的一个简单例子:

在普通的DFT中, 我们需要计算:

$$\left\{ \begin{aligned} &g_0 = f_0 w^{0} + f_1 w^{0} + f_2 w^{0} + f_3 w^{0} \\ &g_1 = f_0 w^{0} + f_1 w^{1} + f_2 w^{2} + f_3 w^{3} \\ &g_2 = f_0 w^{0} + f_1 w^{2} + f_2 w^{4} + f_3 w^{6} \\ &g_3 = f_0 w^{0} + f_1 w^{3} + f_2 w^{6} + f_3 w^{9} \\ \end{aligned}\right.$$

而利用了FFT的思想后, 我们的计算变成

$$\left\{ \begin{aligned} g_0 = (f_0 + w^{0}f_2) + w^{0}(f_1 + w^{0}f_3) \\ g_1 = (f_0 - w^{0}f_2) + w^{1}(f_1 - w^{0}f_3) \\ g_2 = (f_0 + w^{0}f_2) + w^{2}(f_1 + w^{0}f_3) \\ g_3 = (f_0 - w^{0}f_2) + w^{3}(f_1 - w^{0}f_3) \\ \end{aligned} \right.$$

在普通的DFT中, 我们需要 $4 \times 4 = 16$ 次乘法和 $(4-1) \times 4 = 12$ 次加法, 而利用FFT的思想后, 我们只需要6次乘法和8次加法.

---

下面是julia语言的简单FFT实现：

function fft(f, w)
    N = length(f)
    if N == 1
        return f
    end

    #一个julia特别的语法: 数组索引从1开始
    f_even = fft(f[1 : 2 : end], w ^ 2) #计算偶数项
    f_odd = fft(f[2 : 2 : end], w ^ 2)  #计算奇数项
    g = zeros(Complex{Float64}, N)
    for k = 1 : N / 2
        g[k] = f_even[k] + w^(k-1) * f_odd[k]
        g[k + N / 2] = f_even[k - (N / 2)] - w^(k-1) * f_odd[k - (N / 2)]
    end
    return g
end

这里有一个问题: 上面的代码仅在 $N$ 是2的幂时才能正确运行. 为了解决这个问题, 我们可以在输入的数组后面补零, 使其长度为 $2^{M}$. 这种方法被称为零填充(zero padding). 零填充使得频域中的频率槽(frequency bin)变得更加密集, 使得计算出的频域振幅图像显得更加光滑, 它不会改变FT之后频域的图像.

---

关于零填充方法几个可能可以带给你帮助的地方:

https://dsp.stackexchange.com/questions/741/why-should-i-zero-pad-a-signal-before-taking-the-discrete-fourier-transform

https://dsp.stackexchange.com/questions/10293/zero-padding-in-fft