概要：

FFT(Fast Fourier transform)：快速傅里叶变换，是DFT的工程化实现方法。
DFT直接求解太过于复杂，FFT方法根据DFT求解过程中旋转因子的性质并引入分治算法思想，大大简化计算过程，被广泛应用在频谱分析的工程实践中，如matlab，C，C++，CUDA等底层实现

一，DFT简介

频谱分析是信号处理中的重要环节，从傅里叶变换FT，到拉普拉斯变换LT，离散时间傅里叶变换DTFT，Z变换ZT，到我们所讲的离散傅里叶变换DFT（他们之间的联系和区别见我的其他博客）。
相比于其他变换，DFT被广泛应用的原因是其输入的时域信号是离散的，输出的频域结果也是离散的。这就极大方便了我们进行基于计算机的频谱计算，存储和分析，没办法数字信号处理是大趋势。
​ DFT变换的公式为：
X [ k ] = ∑ n = 0 N − 1 x [ n ] e − j k 2 π K n , k = 0 , 1 , 2 , . . . , K − 1 X\left[ k \right]=\sum\limits_{n=0}^{N-1}{x\left[ n \right]{{e}^{-jk\frac{2\pi }{K}n}},k=0,1,2,...,K-1} X[k]=n=0∑N−1​x[n]e−jkK2π​n,k=0,1,2,...,K−1
这里不同于一般课本上的是， k k k的取值不再与输入信号 n n n的长度 0 ∼ N − 1 0\sim N-1 0∼N−1相同，而是自己设置。这是为了突出 k k k的设置本质上是为了对以上 2pi 为周期的连续频谱离散化（DFT是DTFT连续频域结果离散化处理后的结果），也即频谱采样。

但为了分析方便，在FFT的计算过程中，我们依然使用 k = 0 ∼ N − 1 k = 0 \sim N-1 k=0∼N−1的选取策略。也即，如下：

X [ k ] = ∑ n = 0 N − 1 x [ n ] e − j k 2 π N n , k = 0 , 1 , 2 , . . . , N − 1 X\left[ k \right]=\sum\limits_{n=0}^{N-1}{x\left[ n \right]{{e}^{-jk\frac{2\pi }{N}n}},k=0,1,2,...,N-1} X[k]=n=0∑N−1​x[n]e−jkN2π​n,k=0,1,2,...,N−1

上式直接求解当然可以，但是需要 N 2 N^2 N2次复数乘法和 N ( N − 1 ) N(N-1) N(N−1)次复数加法：
X [ k ] = ∑ n = 0 N − 1 x [ n ] e − j k 2 π N n = ∑ n = 0 N − 1 x [ n ] W N k n X\left[ k \right]=\sum\limits_{n=0}^{N-1}{x\left[ n \right]{{e}^{-jk\frac{2\pi }{N}n}}\text{=}}\sum\limits_{n=0}^{N-1}{x\left[ n \right]W_{N}^{kn}} X[k]=n=0∑N−1​x[n]e−jkN2π​n=n=0∑N−1​x[n]WNkn​
其中， W W W是需要替换的旋转因子： W N k n = e − j k 2 π N n , W = e − j 2 π N W_{N}^{kn}={{e}^{-jk\frac{2\pi }{N}n}},W={{e}^{-j\frac{2\pi }{N}}} WNkn​=e−jkN2π​n,W=e−jN2π​

​ 使用FFT算法简化DFT计算过程就是依赖旋转因子 W W W的一些性质，简化计算过程。

二、旋转因子 W W W的性质

• 周期性： W N a + N = W N a W_{N}^{a+N}=W_{N}^{a} WNa+N​=WNa​
• 对称性： W N a + N / 2 = − W N a W_{N}^{a+N/2}=-W_{N}^{a} WNa+N/2​=−WNa​
• 缩放性： W N 2 = W N / 2 1 W_{N}^{2}\text{=}W_{N/2}^{1} WN2​=WN/21​

证明方法就是按旋转因子定义，直接拆开就行，就是代数变换。以上性质意味着值相同的就不用了再多计算一遍了，这就能简化DFT的计算过程。

三、FFT蝶形计算证明

​ FFT的计算过程运用了“分治算法”思想，并结合了旋转因子 W W W的性质。具体证明过程如下：

首先，我们把输入的时域信号 x [ n ] , n = 0 , 1 , . . . , N − 1 x[n],n=0,1,...,N-1 x[n],n=0,1,...,N−1根据索引分为奇偶两部分：

f e v e n [ n ] = x [ 2 n ] {{f}_{even}}\left[ n \right]=x\left[ 2n \right] feven​[n]=x[2n]

f o d d [ n ] = x [ 2 n + 1 ] {{f}_{odd}}\left[ n \right]=x\left[ 2n+1 \right] fodd​[n]=x[2n+1]

此时，索引范围为： n = 0 , 1 , 2 , . . . , N / 2 − 1 n=0,1,2,...,N/2-1 n=0,1,2,...,N/2−1，

​ 2. 对DFT公式（3）进行化简：
X [ k ] = ∑ n = 0 N − 1 x [ n ] W N k n , k = 0 , 1 , . . . , N − 1 X\left[ k \right]=\sum\limits_{n=0}^{N-1}{x\left[ n \right]W_{N}^{kn}},k=0,1,...,N-1 X[k]=n=0∑N−1​x[n]WNkn​,k=0,1,...,N−1
得到：
X [ k ] = ∑ n 为 偶 数 x [ n ] W N k n + ∑ n 为 奇 数 x [ n ] W N k n k = 0 , 1 , . . . , N − 1 X\left[ k \right]=\sum\limits_{n为偶数}{x\left[ n \right]W_{N}^{kn}}+\sum\limits_{n为奇数}{x\left[ n \right]W_{N}^{kn}}k=0,1,...,N-1 X[k]=n为偶数∑​x[n]WNkn​+n为奇数∑​x[n]WNkn​k=0,1,...,N−1
这个过程很简单，如下图，就是按奇偶索引把求和分成两部分

详细点写也即：
X [ k ] = ∑ m = 0 ( N / 2 ) − 1 x [ 2 m ] W N k ⋅ 2 m + ∑ m = 0 ( N / 2 ) − 1 x [ 2 m + 1 ] W N k ⋅ ( 2 m + 1 ) , k = 0 , 1 , . . . , N − 1 X\left[ k \right]=\sum\limits_{m=0}^{(N/2)-1}{x\left[ 2m \right]W_{N}^{k\centerdot 2m}}+\sum\limits_{m=0}^{(N/2)-1}{x\left[ 2m+1 \right]W_{N}^{k\centerdot \left( 2m+1 \right)}},k=0,1,...,N-1 X[k]=m=0∑(N/2)−1​x[2m]WNk⋅2m​+m=0∑(N/2)−1​x[2m+1]WNk⋅(2m+1)​,k=0,1,...,N−1
根据旋转因子的缩放性，可以进一步换算：
X [ k ] = ∑ m = 0 ( N / 2 ) − 1 f e v e n [ m ] W N / 2 k ⋅ m + W N k ∑ m = 0 ( N / 2 ) − 1 f o d d [ m ] W N / 2 k ⋅ m , k = 0 , 1 , . . . , N − 1 X\left[ k \right]=\sum\limits_{m=0}^{(N/2)-1}{{{f}_{even}}\left[ m \right]W_{N/2}^{k\centerdot m}}+W_{N}^{k}\sum\limits_{m=0}^{(N/2)-1}{{{f}_{odd}}\left[ m \right]W_{N/2}^{k\centerdot m}},k=0,1,...,N-1 X[k]=m=0∑(N/2)−1​feven​[m]WN/2k⋅m​+WNk​m=0∑(N/2)−1​fodd​[m]WN/2k⋅m​,k=0,1,...,N−1
也即：
X [ k ] = F e v e n [ k ] + W N k F o d d [ k ] , k = 0 , 1 , . . . , N − 1 X\left[ k \right]={{F}_{even}}\left[ k \right]+W_{N}^{k}{{F}_{odd}}\left[ k \right],k=0,1,...,N-1 X[k]=Feven​[k]+WNk​Fodd​[k],k=0,1,...,N−1
其中， F e v e n [ k ] {{F}_{even}}\left[ k \right] Feven​[k]为偶数索引输入 f e v e n [ n ] {{f}_{even}}\left[ n \right] feven​[n]的DFT结果， F o d d [ k ] {{F}_{odd}}\left[ k \right] Fodd​[k]为奇数索引输入 f o d d [ n ] {{f}_{odd}}\left[ n \right] fodd​[n]的DFT结果。

分析 F e v e n [ k ] {{F}_{even}}\left[ k \right] Feven​[k]和 F o d d [ k ] {{F}_{odd}}\left[ k \right] Fodd​[k]以便进一步转换：

无论 F e v e n [ k ] {{F}_{even}}\left[ k \right] Feven​[k]和 F o d d [ k ] {{F}_{odd}}\left[ k \right] Fodd​[k]两个哪一个，他们的时域输入长度都为 ( N / 2 ) − 1 (N/2)-1 (N/2)−1，但此时的$ k=0,1,…,N-1 是 输 入 信 号 长 度 的 2 倍 。 这 就 说 明 是输入信号长度的2倍。这就说明 是输入信号长度的2倍。这就说明{{F}{even}}\left[ k \right] 和 和 和{{F}{odd}}\left[ k \right]$都是周期性的（可以理解为，N个点里包含了2个$2\pi $周期的频谱采样），也即:

F e v e n [ k + N / 2 ] = F e v e n [ k ] {{F}_{even}}\left[ k+N/2 \right] = {{F}_{even}}\left[ k \right] Feven​[k+N/2]=Feven​[k]

F o d d [ k + N / 2 ] = F o d d [ k ] {{F}_{odd}}\left[ k+N/2 \right] = {{F}_{odd}}\left[ k \right] Fodd​[k+N/2]=Fodd​[k]

再次简化：

根据式（13）可以知道，最终结果的前半部分，可以直接被得到：
X [ k ] = F e v e n [ k ] + W N k F o d d [ k ] , k = 0 , 1 , . . . , N / 2 − 1 X\left[ k \right]={{F}_{even}}\left[ k \right]+W_{N}^{k}{{F}_{odd}}\left[ k \right],k=0,1,...,N/2-1 X[k]=Feven​[k]+WNk​Fodd​[k],k=0,1,...,N/2−1
又通过式（14,15），可以将（16）进一步转化：
X [ k + N / 2 ] = F e v e n [ k ] + W N k F o d d [ k ] , k = 0 , 1 , . . . , N / 2 − 1 X\left[ k+N/2 \right]={{F}_{even}}\left[ k \right]+W_{N}^{k}{{F}_{odd}}\left[ k \right],k=0,1,...,N/2-1 X[k+N/2]=Feven​[k]+WNk​Fodd​[k],k=0,1,...,N/2−1
通过，（16,17）可以看出：一个N点的DFT结果，可以被两个奇偶输入的DFT结果计算得到。举个例子也即，8点的DFT，可以被偶4点DFT结果和奇4点DFT结果计算得到，同理奇/偶4点DFT又可以被2点DFT结果计算得到，以此类推，分治求解。

四、FFT计算过程

​ 步骤1：通过二进制镜像的方法，对时域信号的索引进行二进制编号，如下表最右列，从右向左反推输入计算序列，结果可以对应上图。

二进制：对应计算序列转换原始索引：对应二进制

000:0	←	0:000
100:4	←	1:001
010:2	←	2:010
110:6	←	3:011
001:1	←	4:100
101:5	←	5:101
011:3	←	6:110
111:7	←	7:11

​ 步骤2：奇偶项逐渐合并计算

五、其他说明

• 逆DFT过程也可以使用以上方法计算
• 以上方法为基2的方法，还有基4的方法，具体参考《数字信号处理——原理、算法与应用（第四版）》P380
• CSDN对Typora的md语法支持并不好，有些公式显示不出来，也没有标号，我已把PDF文档上传到了CSDN。