快速傅里叶变换（FFT），真的很细

一、前言

在电赛中，使用FFT算法进行信号频谱分析极其常用，为了给大家科普FFT，本博客将从傅里叶级数到傅里叶变换，再到离散时间傅里叶变换、离散傅里叶变换，之后再简单讲解FFT中简单的按时间抽选的基 2–FFT 算法，主要目的在于为小白科普，如有不严谨或不正确的地方，还请大佬指出。本次博客大量参考其余大佬的博客，笔者只是拾人牙慧的小屁孩。

二、傅里叶变换的前世今生

傅里叶级数展开告诉我们，任何一个函数 $f(t)$ ，都可以展开成如下形式的傅里叶级数:
$$\begin{gathered} f(t)=\frac{1}{2}{a}_0+\sum _{n=1}^{+\infty }[a_{n}cos(nwt)+b_{n}sin(nwt)] \\ w=\frac{2\pi }{T} \\ {a}_n=\frac{2}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)cos(nwt)dt,n=1,2,\mathrm{3...} \\ {b}_n=\frac{2}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)sin(nwt)dt,n=1,2,\mathrm{3...} \end{gathered}$$
关于傅里叶级数待定系数的求解是利用的三角函数的正交性，详情可见这篇知乎。
我们利用欧拉公式得到（傅里叶变换里面有个$j$，要考虑到这个复数和三角函数的关系）：
$$\begin{gathered} cos(nwt)=\frac{e^{jnwt}+e^{-jnwt}}{2} \\ sin(nwt)=\frac{e^{jnwt}-e^{-jnwt}}{2j} \end{gathered}$$
代入傅里叶级数之中：
$$\begin{gathered} f(t)=\frac{1}{2}{a}_0+\sum _{n=1}^{+\infty }[a_{n}cos(nwt)+b_{n}sin(nwt)] \\ f(t)=\frac{1}{2}{a}_0+\sum _{n=1}^{+\infty }(a_n\frac{e^{jnwt}+e^{-jnwt}}{2}+b_n\frac{e^{jnwt}-e^{-jnwt}}{2j}) \\ f(t)=\frac{1}{2}{a}_0+\sum _{n=1}^{+\infty }(\frac{a_n-j{b}_n}{2}{e}^{jnwt}+\frac{a_n+j{b}_n}{2}{e}^{-jnwt}) \end{gathered}$$
当$n=0$时，
$$\frac{a_n-j{b}_n}{2}=\frac{a_0-j{b}_0}{2}=\frac{1}{2}(\frac{2}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)dt-0)=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)dt$$
令$c_0=\frac{a_0-j{b}_0}{2}=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)dt$，
猜想$c_n=\frac{a_n-j{b}_n}{2}=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-jnwt}dt$
假设我们猜想成立，
我们计算$c_{-n}$：
$$\begin{array}{rlr}{c}_{-n}& =\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{jnwt}dt& \\ & =\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)[cos(nwt)+jsin(nwt)]dt& \\ & =\frac{1}{2}[\frac{2}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)cos(nwt)dt+\frac{2j}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)sin(nwt)dt]& \\ & =\frac{a_n+j{b}_n}{2}& \end{array}$$
综上，
$$c_{-n}=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{jnwt}dt,n=1,2,\mathrm{3...}$$
延拓为正整数 $n$，
$$c_n=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-jnwt}dt,n=1,2,\mathrm{3...}$$
之后我们再将傅里叶级数的式子写成：
$$\begin{array}{rl}f(t)& =\frac{1}{2}(a_0+0)+\sum _{n=1}^{+\infty }(\frac{a_n-j{b}_n}{2}{e}^{jnwt}+\frac{a_n+j{b}_n}{2}{e}^{-jnwt})\\ f(t)& =c_0+\sum _{n=1}^{+\infty }(c_n{e}^{jnwt}+c_{-n}{e}^{-jnwt})\\ f(t)& =\sum _{n=-\infty }^{+\infty }{c}_n{e}^{jnwt}\end{array}$$
又因为$c_n=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-jnwt}dt$
得到
$$\begin{array}{rl}f(t)& =\sum _{n=-\infty }^{+\infty }[\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(\tau )e^{-jnw\tau }d\tau ]e^{jnwt}\\ f(t)& =\frac{1}{T}\sum _{n=-\infty }^{+\infty }[{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(\tau )e^{-jnw\tau }d\tau ]e^{jnwt}\end{array}$$
请注意，这里的$n$依然是一个整数，尽管它的范围为$-\infty <n<\infty$，所以接下来我们需要将其连续化。
傅里叶积分定理——当$T\to +\infty$，任何一个非周期函数$f(t)$都可以看成某个周期函数$f_T(t)$转换而来，数学表示即：
$$\underset{T\to +\infty }{\lim }{f}_T(t)=f(t)$$
故，当我们取极限时，
$$\begin{gathered} \underset{T\to +\infty }{\lim }\Delta (nw)=\underset{T\to +\infty }{\lim }[nw-(n-1)w]=\underset{T\to +\infty }{\lim }w \\ =\underset{T\to +\infty }{\lim }\frac{2\pi }{T}=0 \\ i.e.\underset{T\to +\infty }{\lim }\Delta (nw)=0=\underset{T\to +\infty }{\lim }d(nw) \end{gathered}$$
这样，只要$T\to +\infty$，$\Delta (nw)\to d(nw)\to 0$，这样的话，岂不是$nw$的间距就变得非常小了，离散的$\Delta (nw)$就变成了连续的$d(nw)$，求和号就变成了积分号。
我们刚刚的式子也可以写成，
$$\begin{array}{rl}f(t)& =\frac{1}{T}\sum _{n=-\infty }^{+\infty }[{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(\tau )e^{-jnw\tau }d\tau ]e^{jnwt}\\ f(t)& =\underset{T\to +\infty }{\lim }\frac{1}{\frac{2\pi }{dw}}{\int }_{-\infty }^{+\infty }[{\int }_{-\infty }^{+\infty }f(\tau )e^{-jw\tau }d\tau ]e^{jwt}\\ f(t)& =\underset{T\to +\infty }{\lim }\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }[{\int }_{-\infty }^{+\infty }f(\tau )e^{-jw\tau }d\tau ]e^{jwt}dw\end{array}$$
如果我们把${\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(\tau )e^{-jw\tau }d\tau$设为是傅里叶变换，那它可以得到一个含$w$的式子，用于外层积分。于是我们把这个一个函数写成两步：
$$\begin{gathered} F(w)={\int }_{-\infty }^{+\infty }f(t)e^{-jwt}dt \\ f(t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }F(w)e^{jwt}dw \end{gathered}$$
这便是傅里叶变化和逆傅里叶变换。

三、DTFT和DTF

有了傅里叶变换
$$F(w)={\int }_{-\infty }^{+\infty }f(t)e^{-jwt}dt$$
我们把时间离散化，便有了离散时间傅里叶变换（DTFT）：
$$F(w)={\int }_{-\infty }^{+\infty }f(t)e^{-jwn}dn,n\in Z$$
不过，在这里的$n$是从$-\infty$到$\infty$的整数，而且$w$是连续的。该函数在时域上离散，在频域上则是周期的，它一般用来对离散时间信号进行频谱分析。

为了在科学计算和数字信号处理等领域使用计算机进行傅里叶变换，必须将函数定义在离散点上而非连续域内，且须满足有限性或周期性条件，即离散傅里叶变换（DFT）。

把$0$-$N$个点的波形，周期填充到$-\infty$到$\infty$，即为DFT的核心思想。

$wt=2\pi ft=2\pi \frac{t}{T}$，这代表着$t$从$0$到$T$中有一个周期的波形，
而
$$\begin{gathered} 0\le \frac{n}{N}<1 \\ 0\le 2\pi \frac{n}{N}<2\pi \text{ 这代表着0-N中有一个周期的波形} \\ 0\le 2\pi \frac{kn}{N}<2k\pi \text{ 这代表着0-N中有k个周期的波形} \end{gathered}$$
于是我们把DTFT的$wn$改写成$\frac{2\pi kn}{N}$，它们都具有相同的含义。
于是我们感性的得到了DFT（本博客重点不是讲明如何严谨证明，而是通俗易懂解释）：
$$F(k)=\sum _{n=0}^{N-1}{x}_n{e}^{-j2\pi kn/N}$$

四、FFT的蝶形变换

DFT 运算量从$N^2$，因为复数乘法的运算量为：$(a+jb)(c+jd)=(ab-bd)+j(bc+ad)$。
在多点数的DFT运算中，常常会有很长的执行时间。
于是，在 1965 年Cooley 和 Tukey 在《计算机科学 》发表著名的《机器计算傅立叶级数的一种算法》论文后，FFT 开始大规模应用。
FFT也就是Fast Fourier Transform。FFT 充分利用了 DFT 运算中的对称性和周期性,从而将 DFT 运算量从$N^2$ 减少到${\log }_{2}N$。

（一）对称性、周期性和可约性

在式子中，
$$F(k)=\sum _{n=0}^{N-1}f(n)e^{-j2\pi kn/N}$$
我们一般将$e^{-j2\pi kn/N}$写成$W_N^{kn}$
DFT也可以写成
$$F(k)=\sum _{n=0}^{N-1}f(n)W_N^{kn}$$
为什么要这么做，因为$W_N^{kn}$具有很多性质，可以改善DFT 的运算效率，故为了我们方便研究$W_N^{kn}$。

对称性
$(W_N^{kn}{)}^{\ast }=(\cos (-2\pi kn/N)+j\sin (-2\pi kn/N){)}^{\ast }=(\cos (2\pi kn/N)-j\sin (2\pi kn/N){)}^{\ast }=\cos (2\pi kn/N)-j\sin (2\pi kn/N)=W_N^{-kn}$
$W_N^{kn}=W_N^{kn}\times 1=W_N^{kn}\times \cos (-2\pi kN/N)+j\sin (-2\pi kN/N)=W_N^{kn}\times W_N^{kN}=W_N^{k(n+N)}$
$W_N^{kn}=W_N^{kn}\times 1=W_N^{kn}\times \cos (-2\pi Nn/N)+j\sin (-2\pi Nn/N)=W_N^{kn}\times W_N^{Nn}=W_N^{n(k+N)}$

周期性
$W_N^{kn}=W_N^{(n+N)k}$

可约性
$W_N^{kn}=W_{mN}^{mkn}$

特殊点
$W_N^0=1$
$W_N^{N/2}=-1$
$W_N^{k+N/2}=-W_N^k$

（二）FFT的核心思想

FFT 算法的基本思想：
➢ 利用 DFT 系数的特性，合并 DFT 运算中的某些项。
➢ 把长序列 DFT→短序列 DFT，从而减少运算量。
将一个大点数 N 的 DFT 能分解为若干小点数 DFT 的组合，则显然可以达到减少运算工作量的效果。

（三）按时间抽选的基 2–FFT 算法

设输入序列长度为 N = $2^M$（M 为正整数），将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基 2 按时间抽取的 FFT 算法。也称为 Coolkey-Tukey 算法。
其中基 2 表示：N = $2^M$，M 为整数。若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到 N = $2^M$。
下面我们来细说，
$$\begin{array}{rl}F(k)=& DFT[f(n)]\\ =& \sum _{n=0}^{N-1}f(n)W_N^{kn}\\ =& \sum _{r=0}^{N/2-1}f(2r)W_N^{k2r}+\sum _{r=0}^{N/2-1}f(2r+1)W_N^{k(2r+1)}\end{array}$$
把式子分成奇数项和偶数项，接着开始变换，
$$\begin{array}{rl}=& \sum _{r=0}^{N/2-1}f(2r)W_N^{k2r}+\sum _{r=0}^{N/2-1}f(2r+1)W_N^{k(2r+1)}\\ =& \sum _{r=0}^{N/2-1}f(2r)W_N^{k2r}+W_N^k\sum _{r=0}^{N/2-1}f(2r+1)W_N^{k2r}\\ =& \sum _{r=0}^{N/2-1}f(2r)W_{N/2}^{kr}+W_N^k\sum _{r=0}^{N/2-1}f(2r+1)W_{N/2}^{kr}\\ =& F_1(k)+W_N^k{F}_2(k)\end{array}$$
我们便将一个DFT拆分成了两个DFT运算的相加。

我们再利用对称性$W_N^{kn}=W_N^{k(n+N)}$，可以得到，
$$\begin{array}{rl}=& \sum _{r=0}^{N/2-1}f(2r)W_{N/2}^{kr}+W_N^k\sum _{r=0}^{N/2-1}f(2r+1)W_{N/2}^{kr}\\ =& \sum _{r=0}^{N/2-1}f(2r)W_{N/2}^{k(r+N/2)}+W_N^k\sum _{r=0}^{N/2-1}f(2r+1)W_{N/2}^{k(r+N/2)}\\ =& F_1(k+N/2)+W_N^k{F}_2(k+N/2)\end{array}$$
做一些变换，
$$\begin{array}{rl}F(k)=& F_1(k+N/2)+W_N^k{F}_2(k+N/2)\\ F(k+N/2)=& F_1(k+N)+W_N^{k+N/2}{F}_2(k+N)\\ F(k+N/2)=& F_1(k)+W_N^k\times W_N^{N/2}{F}_2(k)\text{ 周期性}\\ F(k+N/2)=& F_1(k)-W_N^k{F}_2(k)\text{ 特殊点}\end{array}$$
于是我们有了两个计算结果，分别结算FFT的前半部分和后半部分
$$\begin{gathered} F(k)=F_1(k)+W_N^k{F}_2(k) \\ F(k+N/2)=F_1(k)-W_N^k{F}_2(k) \end{gathered}$$
有了上面的计算结果后，我们可以得到如下的蝶形运算流图符号：

1. 1 个蝶形运算需要 1 次复乘，2 次复加。(加减都叫复加，这里的复乘和复加意为复数运算)
2. 左边两路为输入，右边两路为输出。
3. 中间以一个小圆表示加减运算（右上路为相加输出，右下路为相减输出）。

分解后的运算量减少了近一半。

我们举个例子：
N=$2^3$=8，
我们如果只分两组，

确实能得到答案，不过还可以继续分下去，提高性能。
再继续分下去的时候，要注意新的奇偶次序，

理清次序之后，我们分四组，进行2点DFT：

下图是按 8 点抽取的 FFT 运算流图：

和我们绘制的本质一致，我们只需增加第一级蝶形即可。
至此，我们成功完成了按时间抽选的基 2–FFT 算法的8点FFT。

四、FFT变换的应用

(一)获取信号的频率幅值相位

一个模拟信号，经过 ADC 采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍（要满足奈奎斯特采样定律）。
现在ADC的采样率为Fs，对峰值为A的信号进行采样，采样点个数为N，我们得到的数字信号就是一个N点数数组。我们对其进行FFT变换，获得N点数FFT结果。（为了方便FFT 运算，通常 N 取 2 的整数次方）
现在我们必须知道FFT 变换结果的物理意义：

1. FFT 的结果的每个点的模值就是 A 的 N/2 倍，第一个点就是直流分量，它的模值就是直流分量的 N 倍。
2. 每个点的相位呢，就是在该频率下的信号的相位——计算的话，非常简单直接将该点对应的实数和虚数求出角度即可。
3. 第一个点表示直流分量（即 0Hz），而最后一个点 N-1的再下一个点N（这个点在数组上不存在的）表示采样频率 Fs，故某点n的频率可以表示为$f_n=\frac{(n-1)\times F_s}{N}$。（含义为0-Fs被分成了N份）
4. 按照3所说，分辨率则为$f_n=\frac{F_s}{N}$。（例如采样频率 Fs为 1024Hz，采样点数为1024 点，则可以分辨到 1Hz。）

可以得出结论，采样时间越长，点数越多，分辨率越高。
但请注意要满足奈奎斯特采样定律（采样频率要大于信号
频率的两倍），当然实际应用中，往往高于两倍。

(二)频谱泄漏

还记得我们计算信号频率时，是将0-Fs被分成了N份，而每个数组下标对应的频率是一个Fs/N的整数倍。
但是实际中没有这么巧合的情况，我们在进行离散化处理的时候，难免会出现截断。

对于频率为 fs 的正弦序列，它的频谱应该只是在 fs 处有离散谱。但是，在利用 DFT 求它的频谱做了截短，结果使信号的频谱不只是在 fs 处有离散谱，而是在以 fs 为中心的频带范围内都有谱线出现，它们可以理解为是从 fs 频率上“泄露”出去的，这种现象称 为频谱“泄露"。

所以在刚刚的图中，正弦的频谱应该是一条直线，但笔者却绘制了一个粗粗的刺，这就是频谱泄露的体现之一。
在实际问题中遇到的离散时间序列f(n)通常是无限长序列（前文提到的DTFT），因而处理这个序列的时候需要将它截断（计算器无法处理无限长序列）。截断相当于将序列乘以窗函数 w(n)。（时域中 f(n)和 w(n)相乘对应于频域中它们的离散傅立叶变换 F(jw)和 W(jw)的卷积）
为了减小频谱“泄露”的影响，往往在 FFT 处理中采用加窗技术，典型的加窗序列有 Hamming、Blackman、Gaussian 等窗序列。此外，增加窗序列的长度（增加采样个数）也可以减少频谱“泄露”。
基础知识补充链接——什么是窗函数?
时域上乘上窗函数，相当于频域进行卷积。长度为无穷长的常数窗函数，频域为 delta 函数，卷积后的结果和原来一样。窗的频谱，越像 delta 函数（主瓣越窄，旁瓣越小），频谱的还原度越高。加窗就不可避免频谱泄漏，典型的加权序列有 Hamming、Blackman、Gaussian 等窗序列主要是为了降低降低旁瓣，对于降低频谱泄漏效果远不如增加窗序列的长度明显。

参考资料

1. 从傅里叶级数到傅里叶变换
2. 傅里叶系列（一）傅里叶级数的推导
3. 安富莱_STM32-V7开发板_第2版DSP数字信号处理教程（V2.7）
4. 什么是窗函数?