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

文章目录

• 1. FFT推导
• 2. FFT为什么快?
• 3. 一些加速措施
• • 3.1 查表法计算三角函数
  • 3.2 奇偶分解
• 4. FFT代码

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

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

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

1. FFT推导

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

X [ k w 0 ] = ∑ n = < N > x [ n ] e − j k w 0 n X[kw_0]=\sum_{n=<N>}x[n]e^{-jkw_0n} X[kw0​]=n=<N>∑​x[n]e−jkw0​n

x [ n ] = 1 N ∑ k = < N > X [ k w 0 ] e j k w 0 n x[n]=\frac{1}{N}\sum_{k=<N>}X[kw_0]e^{jkw_0n} x[n]=N1​k=<N>∑​X[kw0​]ejkw0​n

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

为了方便表示，用 X [ k ] X[k] X[k]来表示这是DFT的第几个点（毕竟是离散的），对于DFT可以写成更具体的形式：

X [ k ] = ∑ n = 0 N − 1 x [ n ] e − j k w 0 n X[k]=\sum_{n=0}^{N-1}x[n]e^{-jkw_0n} X[k]=n=0∑N−1​x[n]e−jkw0​n

这种形式是更符合计算机表示的。

将上式按照 n n n的奇偶性分成两部分：

X [ k ] = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k w 0 ( 2 n ) + ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k w 0 ( 2 n + 1 )   = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k w 0 ( 2 n ) + e − j k w 0 ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k w 0 ( 2 n )   = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k 2 π N / 2 n + e − j k w 0 ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k 2 π N / 2 n X[k]=\sum_{n=0}^{N/2-1} x[2n]e^{-jkw_0(2n)}+\sum_{n=0}^{N/2-1}x[2n+1]e^{-jkw_0(2n+1)}\\ \space\\ =\sum_{n=0}^{N/2-1} x[2n]e^{-jkw_0(2n)}+e^{-jkw_0}\sum_{n=0}^{N/2-1}x[2n+1]e^{-jkw_0(2n)}\\ \space\\ =\sum_{n=0}^{N/2-1} x[2n]e^{-jk\frac{2\pi}{N/2}n}+e^{-jkw_0}\sum_{n=0}^{N/2-1}x[2n+1]e^{-jk\frac{2\pi}{N/2}n} X[k]=n=0∑N/2−1​x[2n]e−jkw0​(2n)+n=0∑N/2−1​x[2n+1]e−jkw0​(2n+1) =n=0∑N/2−1​x[2n]e−jkw0​(2n)+e−jkw0​n=0∑N/2−1​x[2n+1]e−jkw0​(2n) =n=0∑N/2−1​x[2n]e−jkN/22π​n+e−jkw0​n=0∑N/2−1​x[2n+1]e−jkN/22π​n

注意 k k k的取值范围仍然是 [ 0 , N − 1 ] [0, N-1] [0,N−1]，但是由于

e − j ( k + N ) 2 π N n = e − j k 2 π N n ⋅ e − j 2 n π = e − j k 2 π N n 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)N2π​n=e−jkN2π​n⋅e−j2nπ=e−jkN2π​n

同理

e − j ( k + N / 2 ) 2 π N / 2 n = e − j k 2 π N / 2 n ⋅ e − j 2 n π = e − j k 2 π N / 2 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} e−j(k+N/2)N/22π​n=e−jkN/22π​n⋅e−j2nπ=e−jkN/22π​n

所以不妨令 k 1 k_1 k1​的取值范围为 [ 0 , N / 2 − 1 ] [0, N/2-1] [0,N/2−1]。令

X 1 [ k 1 ] = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k 1 2 π N / 2 n X_1[k_1]=\sum_{n=0}^{N/2-1} x[2n]e^{-jk_1\frac{2\pi}{N/2}n} X1​[k1​]=n=0∑N/2−1​x[2n]e−jk1​N/22π​n

X 2 [ k 1 ] = ∑ n = 0 N / 2 − 1 x [ 2 n + 1 ] e − j k 1 2 π N / 2 n X_2[k_1]=\sum_{n=0}^{N/2-1}x[2n+1]e^{-jk_1\frac{2\pi}{N/2}n} X2​[k1​]=n=0∑N/2−1​x[2n+1]e−jk1​N/22π​n

则

X [ k ] = X 1 [ k 1 ] + e − j k w 0 X 2 [ k 1 ] X[k]=X_1[k_1]+e^{-jkw_0}X_2[k_1] X[k]=X1​[k1​]+e−jkw0​X2​[k1​]

其中， k 1 = k % ( N / 2 ) k_1=k \% (N/2) k1​=k%(N/2)，百分号是取余数。

很显然， X 1 [ k 1 ] X_1[k_1] X1​[k1​]就是 x [ n ] x[n] x[n]的偶数项的DFT，而 X 2 [ k 1 ] X_2[k_1] X2​[k1​]就是 x [ n ] x[n] x[n]的奇数项的DFT，只是基波频率有所变化（毕竟样点数变成了原来的一半）。

如果将上式继续按照这种方式分解。

X 1 [ k 1 ] = ∑ n = 0 N / 2 − 1 x [ 2 n ] e − j k 1 2 π N / 2 n   = ∑ n = 0 N / 4 − 1 x [ 4 n ] e − j k 1 2 π N / 2 ( 2 n ) + ∑ n = 0 N / 4 − 1 x [ 4 n + 2 ] e − j k 1 2 π N / 2 ( 2 n + 1 )   = ∑ n = 0 N / 4 − 1 x [ 4 n ] e − j k 1 2 π N / 4 n + e − j k 1 2 π N / 2 ∑ n = 0 N / 4 − 1 x [ 4 n + 2 ] e − j k 1 2 π N / 4 n X_1[k_1]=\sum_{n=0}^{N/2-1} x[2n]e^{-jk_1\frac{2\pi}{N/2}n}\\ \space\\ =\sum_{n=0}^{N/4-1}x[4n]e^{-jk_1\frac{2\pi}{N/2}(2n)}+\sum_{n=0}^{N/4-1}x[4n+2]e^{-jk_1\frac{2\pi}{N/2}(2n+1)}\\ \space\\ =\sum_{n=0}^{N/4-1}x[4n]e^{-jk_1\frac{2\pi}{N/4}n}+e^{-jk_1\frac{2\pi}{N/2}}\sum_{n=0}^{N/4-1}x[4n+2]e^{-jk_1\frac{2\pi}{N/4}n} X1​[k1​]=n=0∑N/2−1​x[2n]e−jk1​N/22π​n =n=0∑N/4−1​x[4n]e−jk1​N/22π​(2n)+n=0∑N/4−1​x[4n+2]e−jk1​N/22π​(2n+1) =n=0∑N/4−1​x[4n]e−jk1​N/42π​n+e−jk1​N/22π​n=0∑N/4−1​x[4n+2]e−jk1​N/42π​n

同样，可以令 k 2 = k 1 % ( N / 4 ) k_2=k_1\%(N/4) k2​=k1​%(N/4)，即 k 2 k_2 k2​的取值范围为 [ 0 , N / 4 − 1 ] [0, N/4-1] [0,N/4−1]。于是可得：

X 1 [ k 1 ] = X 3 [ k 2 ] + e − j k 1 w 0 2 X 4 [ k 2 ] X_1[k_1]=X_3[k_2]+e^{-jk_1w_02}X_4[k_2] X1​[k1​]=X3​[k2​]+e−jk1​w0​2X4​[k2​]

这样就很显然了，整个公式就是把长的DFT转化为短的DFT，每次把一个长度为N的DFT按照序号的奇偶分为两个长度为 N / 2 N/2 N/2的DFT。一直到最后分成一个数的DFT，一个数的DFT就是这个数本身。

X [ 0 ] = ∑ n = 0 0 x [ n ] e − j 0 = x [ 0 ] X[0]=\sum_{n=0}^{0}x[n]e^{-j0}=x[0] X[0]=n=0∑0​x[n]e−j0=x[0]

当然上面是以计算系统序号从0开始的情况来说的（绝大多数计算机语言的序号都是从0开始的），其实如果计算系统序号从1开始也是一样的：

X [ 1 ] = ∑ n = 1 1 x [ n ] e − j w 0 = x [ 1 ] X[1]=\sum_{n=1}^{1}x[n]e^{-jw_0}=x[1] X[1]=n=1∑1​x[n]e−jw0​=x[1]

因为是一个数的DFT，因此 N = 1 N=1 N=1， w 0 = 2 π / N = 2 π w_0=2\pi/N=2\pi w0​=2π/N=2π。由欧拉公式 e − j w 0 = 1 e^{-jw_0}=1 e−jw0​=1。

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

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

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

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

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

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

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

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

然后看那个系数问题。

根据之前的分解过程，第一次分解后，系数是

e − j k w 0 e^{-jkw_0} e−jkw0​

第二次是

e − j k 1 w 0 2 e^{-jk_1w_02} e−jk1​w0​2

第三次是

e − j k 2 w 0 4 e^{-jk_2w_04} e−jk2​w0​4

对于取值范围， k n k_{n} kn​总是 k n − 1 k_{n-1} kn−1​的一半。

第n次分解，就是

e − j k n w 0 2 n − 1 e^{-jk_nw_02^{n-1}} e−jkn​w0​2n−1

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

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

2. FFT为什么快?

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

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

而使用FFT时，我们需要计算 N + N / 2 + N / 4 + ⋯ + 4 N + N/2 + N/4 + \cdots + 4 N+N/2+N/4+⋯+4（等比数列）次复数幂计算，总共有 l o g 2 N − 1 log_2N-1 log2​N−1项，也就是 2 N ( 1 − ( 1 2 ) l o g 2 N − 1 ) 2N(1-(\frac{1}{2})^{log_2N-1}) 2N(1−(21​)log2​N−1)。复数乘法以及加法都要进行 N l o g 2 N Nlog_2N Nlog2​N次。

3. 一些加速措施

3.1 查表法计算三角函数

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

3.2 奇偶分解

首先看欧拉公式，对于DFT公式中e的幂次方部分，：

e − j k 2 π N n = c o s ( k 2 π N n ) − j s i n ( k 2 π N n ) e^{-jk\frac{2\pi}{N}n}=cos(k\frac{2\pi}{N}n)-jsin(k\frac{2\pi}{N}n) e−jkN2π​n=cos(kN2π​n)−jsin(kN2π​n)

我们首先看一下不同频率的复信号。由于是只看频率，因此可剔除 n n n，只看 k k k。（这里不理解的同学可以回顾一下上一篇讲DFT的文章中通俗理解的那部分）。

e − j k 2 π N = c o s ( k 2 π N ) − j s i n ( k 2 π N ) e^{-jk\frac{2\pi}{N}}=cos(k\frac{2\pi}{N})-jsin(k\frac{2\pi}{N}) e−jkN2π​=cos(kN2π​)−jsin(kN2π​)

从这两个三角函数中可以看到， c o s ( k 2 π N ) cos(k\frac{2\pi}{N}) cos(kN2π​)是关于 k = a N 2 k=a\frac{N}{2} k=a2N​偶对称的，而 s i n ( k 2 π N ) sin(k\frac{2\pi}{N}) sin(kN2π​)是关于 k = a N 2 k=a\frac{N}{2} k=a2N​奇对称的。其中a是整数。

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

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

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

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

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

x E [ n ] = x [ n ] + x [ N − n ] 2   x O [ n ] = x [ n ] − x [ N − n ] 2 x_E[n]=\frac{x[n] + x[N - n]}{2}\\ \space\\ x_O[n]=\frac{x[n] - x[N - n]}{2} xE​[n]=2x[n]+x[N−n]​ xO​[n]=2x[n]−x[N−n]​

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

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

X [ k ] = X r [ k ] + X i [ k ] X[k]=X_r[k]+X_i[k] X[k]=Xr​[k]+Xi​[k]

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

设 X e X_{e} Xe​是 x [ n ] x[n] x[n]的偶数项的FFT结果， X o X_{o} Xo​是奇数项的FFT结果（注意此时奇数项被我们视作纯虚数序列）。 X e r X_{er} Xer​是 X e X_{e} Xe​的实部序列， X e i X_{ei} Xei​是 X e X_{e} Xe​的虚部序列。 X o r X_{or} Xor​是 X o X_{o} Xo​的实部序列， X o i X_{oi} Xoi​是 X o X_{o} Xo​的虚部序列。

X e [ k ] = X e r [ k ] + X e i [ k ] X o [ k ] = X o r [ k ] + X o i [ k ] X_{e}[k] = X_{er}[k] + X_{ei}[k]\\ X_{o}[k] = X_{or}[k] + X_{oi}[k] Xe​[k]=Xer​[k]+Xei​[k]Xo​[k]=Xor​[k]+Xoi​[k]

那么

X r [ k ] = X e r [ k ] + X o r [ k ] X_r[k]=X_{er}[k] + X_{or}[k] Xr​[k]=Xer​[k]+Xor​[k]

由前面的结论， X e r X_{er} Xer​是偶对称的， X o r X_{or} Xor​是奇对称的。故

X e r [ k ] = X r [ k ] + X r [ N 2 − k ] 2   X o r [ k ] = X r [ k ] − X r [ N 2 − k ] 2 X_{er}[k]=\frac{X_r[k] + X_r[\frac{N}{2}-k]}{2}\\ \space\\ X_{or}[k]=\frac{X_r[k] - X_r[\frac{N}{2}-k]}{2} Xer​[k]=2Xr​[k]+Xr​[2N​−k]​ Xor​[k]=2Xr​[k]−Xr​[2N​−k]​

同理

X e i [ k ] = X i [ k ] − X i [ N 2 − k ] 2   X o i [ k ] = X i [ k ] + X i [ N 2 − k ] 2 X_{ei}[k]=\frac{X_i[k] - X_i[\frac{N}{2}-k]}{2}\\ \space\\ X_{oi}[k]=\frac{X_i[k] + X_i[\frac{N}{2}-k]}{2} Xei​[k]=2Xi​[k]−Xi​[2N​−k]​ Xoi​[k]=2Xi​[k]+Xi​[2N​−k]​

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

4. FFT代码

代码请访问我的github：FFT

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