数字信号处理3: 快速傅里叶变换（FFT）（含代码）

在之前的文章《傅里叶变换》中，我们已经推导了连续傅里叶变换和离散傅里叶变换。由于计算机的发展，离散傅里叶变换（DFT）可谓是信号处理的杀手锏。但是离散傅里叶变换计算量巨大，通常在实时信号处理时是无法使用的，直到快速傅里叶变换（FFT）算法被发现。

与DFT不同，FFT是一种算法而非理论，因此它无法只靠一两个公式就描述出来。接下来我们从DFT出发进行FFT的推导。

当然，大多数人可能来看这个博客就只是为了寻找可用的FFT代码，但是肯定还有一些同学对这个过程会感到好奇，因此我还是先讲推导过程，代码在最后一部分。

1. FFT推导

首先大家应该还记得DFT的公式

$$X[k{w}_0]=\sum _{n=<N>}x[n]e^{-jk{w}_{0}n}$$

$$x[n]=\frac{1}{N}\sum _{k=<N>}X[k{w}_0]e^{jk{w}_{0}n}$$

其中，$N$为DFT长度，也就是多少个数据点，$w_0=2\pi /N$为基波频率。首先从第一个公式开始。

为了方便表示，用$X[k]$来表示这是DFT的第几个点（毕竟是离散的），对于DFT可以写成更具体的形式：
$$X[k]=\sum _{n=0}^{N-1}x[n]e^{-jk{w}_{0}n}$$
这种形式是更符合计算机表示的。
将上式按照$n$的奇偶性分成两部分：
$$\begin{gathered} X[k]=\sum _{n=0}^{N/2-1}x[2n]e^{-jk{w}_0(2n)}+\sum _{n=0}^{N/2-1}x[2n+1]e^{-jk{w}_0(2n+1)} \\ =\sum _{n=0}^{N/2-1}x[2n]e^{-jk{w}_0(2n)}+e^{-jk{w}_0}\sum _{n=0}^{N/2-1}x[2n+1]e^{-jk{w}_0(2n)} \\ =\sum _{n=0}^{N/2-1}x[2n]e^{-jk\frac{2\pi }{N/2}n}+e^{-jk{w}_0}\sum _{n=0}^{N/2-1}x[2n+1]e^{-jk\frac{2\pi }{N/2}n} \end{gathered}$$
注意$k$的取值范围仍然是$[0,N-1]$，但是由于
$$e^{-j(k+N)\frac{2\pi }{N}n}=e^{-jk\frac{2\pi }{N}n}\cdot e^{-j2n\pi }=e^{-jk\frac{2\pi }{N}n}$$
同理
$$e^{-j(k+N/2)\frac{2\pi }{N/2}n}=e^{-jk\frac{2\pi }{N/2}n}\cdot e^{-j2n\pi }=e^{-jk\frac{2\pi }{N/2}n}$$
所以不妨令$k_1$的取值范围为$[0,N/2-1]$。令
$$X_1[k_1]=\sum _{n=0}^{N/2-1}x[2n]e^{-j{k}_1\frac{2\pi }{N/2}n}$$
$$X_2[k_1]=\sum _{n=0}^{N/2-1}x[2n+1]e^{-j{k}_1\frac{2\pi }{N/2}n}$$
则
$$X[k]=X_1[k_1]+e^{-jk{w}_0}{X}_2[k_1]$$
其中，$k_1=k\%(N/2)$，百分号是取余数。
很显然，$X_1[k_1]$就是$x[n]$的偶数项的DFT，而$X_2[k_1]$就是$x[n]$的奇数项的DFT，只是基波频率有所变化（毕竟样点数变成了原来的一半）。

如果将上式继续按照这种方式分解。
$$\begin{gathered} X_1[k_1]=\sum _{n=0}^{N/2-1}x[2n]e^{-j{k}_1\frac{2\pi }{N/2}n} \\ =\sum _{n=0}^{N/4-1}x[4n]e^{-j{k}_1\frac{2\pi }{N/2}(2n)}+\sum _{n=0}^{N/4-1}x[4n+2]e^{-j{k}_1\frac{2\pi }{N/2}(2n+1)} \\ =\sum _{n=0}^{N/4-1}x[4n]e^{-j{k}_1\frac{2\pi }{N/4}n}+e^{-j{k}_1\frac{2\pi }{N/2}}\sum _{n=0}^{N/4-1}x[4n+2]e^{-j{k}_1\frac{2\pi }{N/4}n} \end{gathered}$$
同样，可以令$k_2=k_1\%(N/4)$，即$k_2$的取值范围为$[0,N/4-1]$。于是可得：
$$X_1[k_1]=X_3[k_2]+e^{-j{k}_1{w}_{0}2}{X}_4[k_2]$$

这样就很显然了，整个公式就是把长的DFT转化为短的DFT，每次把一个长度为N的DFT按照序号的奇偶分为两个长度为$N/2$的DFT。一直到最后分成一个数的DFT，一个数的DFT就是这个数本身。
$$X[0]=\sum _{n=0}^{0}x[n]e^{-j0}=x[0]$$
当然上面是以计算系统序号从0开始的情况来说的（绝大多数计算机语言的序号都是从0开始的），其实如果计算系统序号从1开始也是一样的：
$$X[1]=\sum _{n=1}^{1}x[n]e^{-j{w}_0}=x[1]$$
因为是一个数的DFT，因此$N=1$，$w_0=2\pi /N=2\pi$。由欧拉公式$e^{-j{w}_0}=1$。

稍微有些跑题，接着讲。我们计算出来了一个数的DFT之后，根据上面分解的那些公式显然可以用它们来组合出两个数的DFT，进而得出四个数的DFT，一直到$N$。其实这整个过程也隐含了一个很重要的条件，就是FFT的长度必须为2的幂次方。

好了，现在我们知道，长的DFT可以通过短的DFT组合计算。但是还有两个头疼的问题：

• 在我们分解的过程中，$x[n]$的顺序已经被打乱了，那么这个顺序要怎么排？
• 每次分解后的奇数项求和前面有个系数，这个系数如何确定？

接下来用一个具体的例子。

假设现在要求一个8个点的DFT。

表中绿色表示每次分解后的偶数项，黄色表示每次分解的奇数项。分解3次后到底。可以预见，对于长度为$N$的序列(N必须为2的幂次方)，分解需要进行$lo{g}_{2}N$步。

乍一看似乎看不出有什么规律，其实这种序号变化可以使用一种位反转的方法来计算。如下表:

需要注意，位反转并非二进制取反，而是“低位与高位交换”，并且反转点需要根据具体的长度来定。比如如果$N=16$时，原序号为1，反转后就是8。

然后看那个系数问题。
根据之前的分解过程，第一次分解后，系数是
$$e^{-jk{w}_0}$$
第二次是
$$e^{-j{k}_1{w}_{0}2}$$
第三次是
$$e^{-j{k}_2{w}_{0}4}$$
对于取值范围，$k_n$总是$k_{n-1}$的一半。
第n次分解，就是
$$e^{-j{k}_n{w}_0{2}^{n-1}}$$

OK，以上基本就是所有FFT的推导了，结合上述过程，我们就可以得到FFT计算的蝶形图：

图中$W_M^k=e^{-jk\frac{2\pi }{M}}$，$k$的取值范围为$[0,M-1]。$这里的$M$并非原序列$x[n]$的长度，而是每次被分解的那个序列的长度。比如在第二次分解时，我们得到$e^{-j{k}_1{w}_{0}2}=e^{-j{k}_1\frac{2\pi }{N/2}}$，$M=N/2=4$。

2. FFT为什么快?

在未使用任何程序上的加速技巧时：

对于一个长度为N的序列，使用传统DFT方法，我们需要计算$N^2$次复数乘法，$N^2$次求和，$N^2$次复数幂的计算（当然你也可以使用欧拉公式将其转化为三角函数的计算，但这仍然是不小的开支）。

而使用FFT时，我们需要计算$N+N/2+N/4+\cdots +4$（等比数列）次复数幂计算，总共有$lo{g}_{2}N-1$项，也就是$2N(1-(\frac{1}{2}{)}^{lo{g}_{2}N-1})$。复数乘法以及加法都要进行$Nlo{g}_{2}N$次。

3. 一些加速措施

3.1 查表法计算三角函数

对于计算机来说，无论是计算e的复数幂，或者求三角函数，都是非常耗时的。在arm加速库中，对于三角函数是采用查表法，预先生成一个周期的sin值保存在数组里，需要精度高就保存数组长一些，使用的时候根据角度在数组中查找对应值。当然，如果你不嫌烦的话，仅生成1/4个周期的sin值即可，不过查找起来比较麻烦。

3.2 奇偶分解

首先看欧拉公式，对于DFT公式中e的幂次方部分，：
$$e^{-jk\frac{2\pi }{N}n}=cos(k\frac{2\pi }{N}n)-jsin(k\frac{2\pi }{N}n)$$
我们首先看一下不同频率的复信号。由于是只看频率，因此可剔除$n$，只看$k$。（这里不理解的同学可以回顾一下上一篇讲DFT的文章中通俗理解的那部分）。
$$e^{-jk\frac{2\pi }{N}}=cos(k\frac{2\pi }{N})-jsin(k\frac{2\pi }{N})$$
从这两个三角函数中可以看到，$cos(k\frac{2\pi }{N})$是关于$k=a\frac{N}{2}$偶对称的，而$sin(k\frac{2\pi }{N})$是关于$k=a\frac{N}{2}$奇对称的。其中a是整数。

因此我们可以得到一个重要结论：

对于纯实信号，DFT后$X[k]$的实部是关于$\frac{N}{2}$偶对称的，$X[k]$的虚部是关于$\frac{N}{2}$奇对称的。
对于纯虚信号，DFT后$X[k]$的实部是关于$\frac{N}{2}$奇对称的，$X[k]$的虚部是关于$\frac{N}{2}$偶对称的。

由于我们接触的基本都是实信号，如此一来如果按照传统FFT，就会有相当多的计算是浪费的。因此，我们不妨把纯实序列$x[n]$当作复数序列，它的偶数项是实部，奇数项是虚部。这样一来我们仅需要做计算$N/2$长度的FFT即可。

待上一步计算完成后，我们需要把实部和虚部的FFT结果分离出来。这个时候就要用到之前对称性的结论，可以知道$X[k]$的实部和虚部都是由一个偶对称序列和一个奇对称序列叠加而成。如何把偶对称序列和奇对称序列分离出来？这就需要使用奇偶分解：

$$\begin{gathered} x_E[n]=\frac{x[n]+x[N-n]}{2} \\ {x}_O[n]=\frac{x[n]-x[N-n]}{2} \end{gathered}$$

N为$x[n]$的长度。同时我们令$x[N]=x[0]$。这对于$X[k]$也是适用的，因为如果不限制k的取值范围，$X[k]$是周期为N的周期函数。这一点可以从DFT公式得出。

接下来我们就可以分解出偶对称序列和奇对称序列。$X[k]$由实部序列和虚部序列组合而成：

$$X[k]=X_r[k]+X_i[k]$$

由于我们把长度为$N$的实数FFT转化为长度为$N/2$的复数FFT，因此上式中k的范围是$[0,N/2-1]$。$X_r[k]$是实部序列，$X_i[k]$是虚部序列。

设$X_e$是$x[n]$的偶数项的FFT结果，$X_o$是奇数项的FFT结果（注意此时奇数项被我们视作纯虚数序列）。$X_{er}$是$X_e$的实部序列，$X_{ei}$是$X_e$的虚部序列。$X_{or}$是$X_o$的实部序列，$X_{oi}$是$X_o$的虚部序列。

$$\begin{gathered} X_e[k]=X_{er}[k]+X_{ei}[k] \\ {X}_o[k]=X_{or}[k]+X_{oi}[k] \end{gathered}$$
那么
$$X_r[k]=X_{er}[k]+X_{or}[k]$$
由前面的结论，$X_{er}$是偶对称的，$X_{or}$是奇对称的。故
$$\begin{gathered} X_{er}[k]=\frac{X_r[k]+X_r[\frac{N}{2}-k]}{2} \\ {X}_{or}[k]=\frac{X_r[k]-X_r[\frac{N}{2}-k]}{2} \end{gathered}$$

同理

$$\begin{gathered} X_{ei}[k]=\frac{X_i[k]-X_i[\frac{N}{2}-k]}{2} \\ {X}_{oi}[k]=\frac{X_i[k]+X_i[\frac{N}{2}-k]}{2} \end{gathered}$$

这样一来，我们就得到了$X_e[k]$和$X_o[k]$。回顾我们将实数FFT转化为复数DFT的那步，我们将偶数项和奇数项分开做FFT，这和蝶形运算非常相似，不过在蝶形运算中，我们仍然将奇数项看作实数，因此，我们将$X_o[k]$每一项乘$-j$，即可得到奇数项的FFT结果，此时它是正常的实数FFT结果。接下来，我们就可以进行蝶形运算的最后一步合并，得到的就是最终的FFT结果。

4. FFT代码

代码请访问我的github：FFT

在代码中，我使用了查表法。对于纯实数的FFT，我将其转换为了复数FFT。在代码中并没有使用复数，而是使用float数组代替复数数组，格式为{real, imag, real ,imag …}。