快速傅里叶变换

前言

终于下定决心要啃掉这块骨头。大四选修数字图像处理的时候，第一次接触傅里叶变换，然后实现了时间度复杂度为 $O(n^2)$ 的算法，后来才知道还有快速傅里叶变换，时间复杂度是 $O(nlogn)$。

拉格朗日等数学家发现某些周期函数可以由三角函数的和表示，而傅里叶则猜测任意周期函数都可以写成三角函数的和，也就是高等数学中的傅里叶级数。大家当然都不相信啦，但又不能给出有力的论据。直到20年后（1829年）狄利克雷才给出了让人满意的答案：只有满足一定条件时，周期函数才能展开成傅里叶级数。这个条件被称为狄利克雷条件。傅里叶级数用吴文俊先生的话说就是：

把质的困难转化成量的复杂

而最早具有这种想法的人是泰勒，他在 1715 年研究出泰勒公式。泰勒级数把函数表达式转化为多项式函数来近似，是求导函数组成的特征函数和，反应变化剧烈程度。傅里叶级数是频谱叠加的三角函数和，反应变化频率本质属性。

实函数的傅立叶级数和泰勒级数不过是观察复幂级数的两种不同方式

傅里叶级数

傅里叶级数具体构造过程可以参考马同学高等数学的文章。假设 $f(x)$ 的周期为 $T$，并且满足狄利克雷条件，则构造结果如下所示：
$$
f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\Big(a_ncos(\frac{2\pi nx}{T})+b_nsin(\frac{2\pi nx}{T})\Big),a_0\in\mathbb{R}
$$
这就是傅里叶级数，其中 $a_n, b_n$ 为傅里叶系数。利用三角函数的正交性等性质（详细推导过程可参考高数课本），可以推导出：
$$
\begin{cases}
a_n=\frac{2}{T}\int_{x_0}^{x_0+T}f(x)cos(\frac{2\pi nx}{T})dx \\
b_n=\frac{2}{T}\int_{x_0}^{x_0+T}f(x)sin(\frac{2\pi nx}{T})dx
\end{cases}
$$
举个简单的例子，有以周期为 $2\pi$ 的三角波函数，在一个周期内的表达式为：
$$
f(x)=\begin{cases}
1+\frac{x}{\pi}, & -\pi\leq x<0 \\
1-\frac{x}{\pi}, & 0\leq x<\pi
\end{cases}
$$
即 $f(x)=1-2\times |nint(\frac{x}{2\pi})-\frac{x}{2\pi}|$，其中 nint() 为向下取整函数，在 Python 中相当于 round() 函数。其函数图（时域图）如下所示：

由于 $f(x)$ 为偶函数，所以正弦分量的幅值 $b_n=0$。解得:
$$
f(x)=\frac{1}{2}+\frac{4}{\pi^2}\Big(cos(x)+\frac{cos(3x)}{3^2}+\frac{cos(5x)}{5^2}+…\Big)
$$
傅里叶级数实际上相当于是把 $f(x)$ 当做了基的向量。当 $n$ 取 5 时该函数的基为：
$$
\lbrace 1, cos(x), cos(2x), cos(3x), cos(5x)\rbrace
$$
该函数相当于向量：
$$
(1, \frac{4}{\pi^2}, 0, \frac{4}{3^2\pi^2}, 0, \frac{4}{5^2\pi^2})
$$

根据这个向量就能得到函数的幅值频谱图，其中纵坐标为幅值，也就是上面的向量。横坐标为角频率（角速度的标量），也就是对应的基中的 $\omega_n=\frac{2\pi n}{T}$，单位为 $Hz$。

非奇非偶函数既有正弦分量也有余弦分量，所以就需要有两个频域图。当 $n$ 取值越大，则越能逼近原函数：

复数形式

傅里叶级数具有两个频域图，而另一种角度，即复数形式的傅里叶级数就可以结合正弦和余弦分量，详情可见马同学高等数学。复数形式的傅里叶级数如下所示：
$$
f(x)=\sum_{n=-\infty}^{\infty}c_ne^{i\frac{2\pi nx}{T}}
$$
其中
$$
c_n=\frac{1}{T}\int_{x_0}^{x_0+T}f(x)e^{-i\frac{2\pi nx}{T}}dx
$$
由于 $c_n$ 是复数，周期函数的频域图还是不好画，所以上面采取的是三角级数。那复数形式有啥用？欸，先来看看非周期函数和傅里叶变换。

傅里叶变换

傅里叶级数只能用于周期函数，后来的数学家将其扩展到非周期函数上，其思想就是让周期 $T$ 趋于无穷。由傅里叶级数公式可知，两个相邻频域直接的差越小，即频率越密集，详情可见马同学高等数学。因此时域上周期的函数其频域是离散的，在时域上非周期的函数其频域是连续的。对于非周期函数，周期趋于无穷，令 $F(\omega_n)=T\times c_n$，有：
$$
\begin{align*}
f(x)&=\lim_{T\rightarrow \infty}\sum_{n=-\infty}^{\infty}\frac{1}{T}F(\omega_n)e^{i\frac{2\pi nx}{T}} \\
&=\lim_{T\rightarrow \infty}\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}F(\omega_n)e^{i\omega_nx}(\omega_n-\omega_{n-1}) \\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i\omega x}d\omega \\
\end{align*}
$$
其中
$$
F(\omega)=\int_{-\infty}^{\infty}f(x)e^{-i\omega x}dx
$$

总结

$F(\omega)$ 就是傅里叶变换，得到的就是频域曲线，即各频率分量的幅值。$f(x)$ 是傅里叶变换的逆变换，所以时域和频域就是看待同一个数学对象的两种角度，可以来回切换。

离散时间傅里叶变换

上面提及的傅里叶变换 $F(\omega)$ 是连续时间傅里叶变换，即信号在时域上是连续的。但是计算机存储的信号都是采样得到的离散的、有限长的信号。在数学上可以用原信号与脉冲函数（具有筛选性）的乘积表示采样，假设采样周期为 $T$，即每隔时间 $T$ 进行一次采样，则 $f[nT]$ 相当于 $f(x)$ 在时间点 $nT$ 上的采样，即：
$$
\begin{align*}
f_{discrete}(x)&=f(x)\sum_{n=-\infty}^{\infty}\delta(x-nT) \\
&=\sum_{n=-\infty}^{\infty}f[nT]\delta(x-nT) \\
\end{align*}
$$
所以
$$
F(e^{i\omega})=\sum_{n=-\infty}^{\infty}f[nT]e^{-i\omega nT}
$$
这就是离散时间傅里叶变换（DTFT），由于离散时间傅里叶变换可以被看作 Z 变换的特例，因此通常用 $F(e^{i\omega})$ 表示离散时间傅里叶变换，其逆变换为：
$$
f[nT]=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(e^{i\omega})e^{i\omega nT}d\omega
$$
由于 $e^{i\theta}$ 的周期是 $2\pi$，因此对于离散时间傅里叶变换，$n$ 为整数，只要采样周期 $T$ 设置为整数，就有：
$$
\begin{align*}
F(e^{i\omega})&=\sum_{n=-\infty}^{\infty}f[nT]e^{-i\omega nT} \\
&=\sum_{n=-\infty}^{\infty}f[nT]e^{-i(\omega nT+2\pi MnT)} \\
&=\sum_{n=-\infty}^{\infty}f[nT]e^{-i(\omega+2\pi M)nT} \\
&=F(e^{i(\omega+2\pi M)}) \\
\end{align*}
$$
其中 $M$ 为整数，所以离散时间傅里叶变换是频率 $\omega$ 的周期函数，周期为 $2\pi$。而对于连续时间傅里叶变换，时间是连续的时间 $x$，则：
$$
\int_{-\infty}^{\infty}f(x)e^{-i\omega x}dx\neq\int_{-\infty}^{\infty}f(x)e^{-i(\omega +2\pi M)x}dx
$$

总结

时域和频域具有以下关系：

时域	频域
周期	离散
非周期	连续
连续	非周期
离散	周期

即大部分的信号为：连续时间周期信号、连续时间非周期信号、离散时间非周期信号和离散时间周期信号。而计算机存储的信号为时域上采样得到的信号，因此在计算机中我们用不到连续时域信号，以下内容均默认为时域离散，即频域是周期的。

离散傅里叶变换

由于计算机只能处理有限长的信号，对于时域周期的信号，则可以对其一个周期内的信号进行傅里叶变换，然后扩展到整个时域上。对于无限长非周期信号就需要对信号进行截断，然后再当成是周期信号的一个周期来处理

例如在 $[0, L]$ 对离散时域的信号进行截断，设采样周期为 $T$，则时域的采样点数为 $N=\frac{L}{T}$，有：
$$
f_{discrete}(x)=\sum_{n=0}^{N-1}f[nT]\delta(x-nT)
$$
截断后的离散时间傅里叶变换为：
$$
F(e^{i\omega})=\sum_{n=0}^{N-1}f[nT]e^{-i\omega nT}
$$
由于非周期信号的频域是连续的，也无法表示，所以也要在频域上采样，就是相当于以截断的这一段信号当成是周期信号的一个周期，对周期信号的一个周期进行傅里叶变换。根据采样定理，频域的采样间隔应为 $\frac{1}{L}$，频域采样点数为：
$$
\frac{1/T}{1/L}=N
$$
即频域采样的点数和时域采样同为 $N$，频域的采样点为 $\lbrace\omega_n=2\pi k/NT\rbrace_{0\leq k<N}$。这就是离散傅里叶变换（DFT），令 $T=1$ 并且将其归一化，得：
$$
F[k]=\sum_{n=0}^{N-1}f[n]e^{-i\frac{2\pi}{N}nk}, k=0, 1, …, N-1
$$
逆变换为：
$$
f[n]=\frac{1}{N}\sum_{k=0}^{N-1}F[k]e^{i\frac{2\pi}{N}nk}, n=0, 1, …, N-1
$$

总结

综上所述，计算机对现实中的信号的时域采样（离散化）得到了离散时间傅里叶变换（DTFT），这个过程中频域会被周期化；然后对信号进行截断，即只收集一段时间的信号；继而对频域进行采样（离散化）得到离散周期傅里叶级数（DFS），这个过程中时域会被周期化。最后只取一个周期进行分析，即只取一个周期内的 $N$ 个采样点。

快速傅里叶变换

1965年，Cooley 和 Tukey 提出了计算离散傅里叶变换（DFT）的快速算法，将 DFT 的运算量减少了几个数量级。快速傅里叶变换（FFT）中的”快速”等价于快速排序中的“快速”。也就是说，我们需要在计算机上实现排序，有各种各样的排序算法，快速排序这个算法比较好用。快速傅里叶变换也一样，利用了分治策略，减小了算法的时间复杂度。对于 $N$ 点信号进行离散傅里叶变换，令 $W_N=e^{-i\frac{2\pi}{N}}$，有：
$$
F[k]=\sum_{n=0}^{N-1}f[n]W_N^{nk}, k=0, 1, …, N-1
$$
直接计算的话，算法的时间复杂度为 $N^2$。快速傅里叶变换算法就能像快速排序算法一样，将时间复杂度从 $N^2$ 降到 $Nlog_2N$。快速傅里叶变换是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的，它对傅里叶变换的理论并没有新的发现，但是对于在计算机中应用离散傅立叶变换，可以说是进了一大步。

大部分快速傅里叶变换是基于 2 的，即信号的长度为 2 的整数幂（不足则补 0）。快速傅里叶算法每次将长度为 $N$ 的信号分解为两个长度为 $\frac{N}{2}$ 的信号进行处理（按奇数和偶数），即：
$$
\begin{align*}
F[k]&=\sum_{n=0}^{N-1}f[n]W_N^{nk} \\
&=\sum_{n=0}^{\frac{N}{2}-1}f[2n]W_N^{(2n)k}+\sum_{n=0}^{\frac{N}{2}-1}f[2n+1]W_N^{(2n+1)k} \\
\end{align*}
$$
因为
$$
W_N^{(2n)k}=e^{-i\frac{2\pi}{N}(2n)k}=e^{-i\frac{2\pi}{\frac{N}{2}}nk}=W_{\frac{N}{2}}^{nk}
$$

所以
$$
F[k]=E[k]+W_N^kO[k], 0\leq k\leq N-1
$$
其中，偶数样本点信号的离散傅里叶变换为：
$$
E[k]=\sum_{n=0}^{\frac{N}{2}-1}f[2n]W_{\frac{N}{2}}^{nk}
$$
奇数样本点信号的离散傅里叶变换为：
$$
O[k]=\sum_{n=0}^{\frac{N}{2}-1}f[2n+1]W_{\frac{N}{2}}^{nk}
$$
由于 $W_N^{nk}$ 对 $n$ 和 $k$ 都具有周期性，即
$$
W_N^{nk}=W_N^{(n+N)k}=W_N^{n(k+N)}
$$
所以
$$
\begin{cases}
E[k]=E[k+\frac{N}{2}] \\
O[k]=O[k+\frac{N}{2}]
\end{cases}
$$

蝶形图

以 $N=8$ 为例，通过蝶形图可视化快速傅里叶具体计算过程：

如上图所示，有：
$$
\begin{cases}
F[0]=E[0]+W_8^0O[0] \\
… \\
F[4]=E[4]+W_8^4O[4]=E[4-\frac{8}{2}]+W_8^4O[4-\frac{8}{2}]=E[0]+W_8^4O[0] \\
…
\end{cases}
$$
同理，可以继续分解

直到最后

从蝶形图中可以看出，一共有 24 次加法操作。由于 $W^0=1$，所以乘法操作的次数小于加法操作的次数。因此快速傅里叶变换的时间复杂度为 $O(Nlog_2N)$。

测试

import numpy as np
import time

def DFT(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

if __name__ == '__main__':
    x = np.random.random(1024)

    start = time.time()
    dft = DFT(x)
    end = time.time()
    print(end - start)

    start = time.time()
    fft = np.fft.fft(x)
    end = time.time()
    print(end - start)

    ret = np.allclose(dft, fft)
    print(ret)

对于一个长度为 1024 的序列，手动实现的离散傅里叶变换算法需要的时间为 0.2748403549194336 秒，而numpy 的快速傅里叶变换算法只需要 0.0009996891021728516 秒，最终计算结果一致。

总结

花了整整一周的时间终于总结完了快速傅里叶变换相关的内容，虽然复数部分不是很深入，但是目前至少能比较熟练得使用这个 20 世纪最伟大的算法了。以上内容之所以花了这么多时间，还是因为概念不清晰，了解傅里叶变换、离散时间傅里叶变换、离散傅里叶变换和快速傅里叶变换就花了不少时间，如果这次没有记录下来，我相信下次看到还是不会！

参考文献

马同学高等数学. 如何理解傅立叶级数公式？. https://www.matongxue.com/madocs/619.html.
马同学高等数学. 从傅立叶级数到傅立叶变换. https://www.matongxue.com/madocs/712.html.