详解FFT算法：从DFT到快速傅里叶变换
详解FFT算法：从DFT到快速傅里叶变换
傅里叶变换是将时域信号转换为频域信号的重要信号处理方法。本文将详细介绍DFT（离散傅里叶变换）、旋转因子特性以及FFT（快速傅里叶变换）算法，并通过具体实例进行验证。
DFT（离散傅里叶变换）
离散傅里叶变换（Discrete Fourier Transform）是数字信号处理中的重要工具。以固定采样频率Fs对输入信号进行等间隔采样N次后，获得N个离散采样点，其中采样间隔为Fs / N（亦可称为采样分辨率，单位Hz）。采集到的N个采样点为输入信号在时域上的振幅信息。DFT对时域上的N个离散点进行处理后会得到N个复数，每个复数表示了以Fs / N为间隔的频率成分的幅度和相位信息，从而完成了离散信号由时域到频域的转换。
注： 根据奈奎斯特采样定理，为了保证信号的采样不失真，当采样频率为Fs时，可采集的最大输入信号频率为Fs / 2。
公式解析
对于一个长度为N的时域离散信号X(n)，其DFT变换结果为X(k)，则：
X(k) = DFT{X(n)} = \sum_{n=0}^{N-1} X[n]e^{-j2\pi kn/N}, k=0,1,...,N-1
其中：
- X[k] 是频域上的第 k 个频率分量，
- X[n] 是时域上的第 n 个样本，
- j 是虚数单位，
- N 是信号的长度，
- e^{-j2\pi kn/N} 是旋转因子，表示频率 k 的复指数函数
注： 欧拉公式 e^{jx} = cos(x) + jsin(x)
DFT举例验证
Eg1: 假设输入信号为F(t)=cos(2Π3t)，以采样频率8Hz采集8个点，得到以下波形数据：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)
0.0 0.71 -1 0.71 0.0 -0.71 1 -0.71
整个DFT的表达式可如下所示，等式左侧即为经过傅里叶变换后的值：
下一步，进行矩阵的分解运算以便于理解其中的运算过程：
X(0) = (0.0 * (cos(0)+sin(0)j)) + (0.71 * (cos(0)+sin(0)j)) + ... + (-0.71 * (cos(0)+sin(0)j))
= 0+0j
X(1) = (0.0 * (cos(0)+sin(0)j)) + (0.71 * (cos(Π/4)-sin(Π/4)j))+ ... + (-0.71 * (cos(7Π/4) + sin(7Π/4)j))
= 0+0j
X(2) = 0 + 0j
X(3) = 0 - 4j
X(4) = 0 + 0j
X(5) = 0 + 4j
X(6) = 0 + 0j
X(7) = 0 + 0j
DFT序列 复数 幅度 相位
x(0) 0 + 0j 0 0
x(1) 0 + 0j 0 0
x(2) 0 + 0j 0 0
x(3) 0 - 4j 4 0
x(4) 0 + 0j 0 0
x(5) 0 + 4j 4 0
x(6) 0 + 0j 0 0
x(7) 0 + 0j 0 0
注： 复数：a + bj，幅度 = \sqrt{a^2 + b^2}，相位 = arctan(b / a)。
由此可以看出x(3)与x(5)存在共轭属性，根据以上运算过程可知，需进行N²次的复数的乘法运算，在进行N(N-1)次的复数的加法运算。
旋转因子特性
n,k = 0,1,2,...,N-1
Eg2: 假设N等于8，我们可以得到以下旋转因子分布图：
周期性：
∵ e^{-j2\pi kn/N} = e^{-j2\pi (k+N)n/N}
∴ e^{-j2\pi kn/N} = e^{-j2\pi kn/N}e^{-j2\pi n} = e^{-j2\pi kn/N}
对称性：
e^{-j2\pi kn/N} = e^{-j2\pi (N-k)n/N}e^{-j2\pi n} = e^{-j2\pi (N-k)n/N}
缩放性：
e^{-j2\pi kn/N} = e^{-j2\pi (k/2)n/(N/2)}e^{-j2\pi kn/N}
FFT（快速傅里叶变换）
FFT算法，即快速傅里叶变换（Fast Fourier Transform）算法，是一种高效计算离散傅里叶变换（DFT）的方法。FFT算法的核心是将一个长度为N的离散傅里叶变换，通过利用旋转因子的周期性、对称性和缩放性实现蝶形迭代运算，从而分解为长度较短的离散傅里叶变换。
输入信号：x(n)，n=0,1,2,...,N-1
经过 X(k) = DFT{X(n)} = \sum_{n=0}^{N-1} X[n]e^{-j2\pi kn/N}，k=0,1,...,N-1，此公式可通过DFT计算得出频率分量。
奇偶部分关系式推导
我们重点理解基2FFT，我们将x(n),按下标分解为奇序列和偶序列然后分解计算，如下所示：
将x(n)分解奇偶序列后，DFT的公式亦可分解为奇偶序列两部分的求和运算，如下所示： 注： 此处
X(k) = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} + \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi (k+1/2)n/N}，k=0,1,...,N-1
其中：
- X[2n] 是偶序列
- X[2n+1] 是奇序列
则公式可变换为：
X(k) = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} + e^{-j2\pi k/N}\sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi kn/N}，k=0,1,...,N-1
根据旋转因子的缩放性，公式可进行如下变换：
X(k) = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/(N/2)} + e^{-j2\pi k/N}\sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi kn/(N/2)}，k=0,1,...,N-1
令：
- G(k) = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/(N/2)}
- H(k) = \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi kn/(N/2)}
则公式可简化为：X(k) = G(k) + e^{-j2\pi k/N}H(k)，k=0,1,2,...,N-1
验证共轭性
根据前面DFT的运算结果，我们可知X(k)前半段与后半段的结果存在共轭属性，因此我们进行公式的验证：
X(k)前 = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} + \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi (k+1/2)n/N}，k=0,1,2,...,N/2-1
X(k)后 = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi (k+N/2)n/N} + \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi (k+N/2+1/2)n/N}，k=0,1,2,...,N/2-1
根据旋转因子的周期性与对称性：
e^{-j2\pi (k+N/2)n/N} = e^{-j2\pi kn/N}e^{-j\pi n} = -e^{-j2\pi kn/N}
e^{-j2\pi (k+N/2+1/2)n/N} = e^{-j2\pi (k+1/2)n/N}e^{-j\pi n} = -e^{-j2\pi (k+1/2)n/N}
则：X(k)后 = -\sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} - \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi (k+1/2)n/N}，k=0,1,2,...,N/2-1
因此：
X(k) = \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} + \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi (k+1/2)n/N}，k=0,1,2,...,N/2-1
X(k+N/2) = -\sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} - \sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi (k+1/2)n/N}，k=0,1,2,...,N/2-1
FFT运算逻辑推导
已知根据DFT公式：X(k) = \sum_{n=0}^{N-1} X[n]e^{-j2\pi kn/N}，k=0,1,2,...,N-1
我们根据旋转因子特性进行不断基二分法，直到N=1.可做以下运算：
X(k) = \sum_{n=0}^{N-1} X[n]e^{-j2\pi kn/N}，k=0,1,2,...,N-1
第一次分解奇偶序列，带入FFT推导公式：
X(k) = X(k) + X(k+N/2)
= \sum_{n=0}^{N/2-1} X[2n]e^{-j2\pi kn/N} + e^{-j2\pi k/N}\sum_{n=0}^{N/2-1} X[2n+1]e^{-j2\pi kn/N}，k=0,1,2,...,N/ 2 - 1
第二次分解奇偶序列，带入FFT推导公式：
k = 0,1,2,..,N/4-1
= \sum_{n=0}^{N/4-1} X[4n]e^{-j2\pi kn/(N/2)} + e^{-j2\pi k/(N/2)}\sum_{n=0}^{N/4-1} X[4n+2]e^{-j2\pi kn/(N/2)} + e^{-j2\pi k/N}\sum_{n=0}^{N/4-1} X[4n+1]e^{-j2\pi kn/(N/2)} + e^{-j2\pi (k+1/2)/N}\sum_{n=0}^{N/4-1} X[4n+3]e^{-j2\pi kn/(N/2)}
第三次分解奇偶序列，带入FFT推导公式：
k = 0,1,2,..,N/8-1
= \sum_{n=0}^{N/8-1} X[8n]e^{-j2\pi kn/(N/4)} + e^{-j2\pi k/(N/4)}\sum_{n=0}^{N/8-1} X[8n+4]e^{-j2\pi kn/(N/4)} + e^{-j2\pi k/(N/2)}\sum_{n=0}^{N/8-1} X[8n+2]e^{-j2\pi kn/(N/4)} + e^{-j2\pi (k+1/2)/(N/2)}\sum_{n=0}^{N/8-1} X[8n+6]e^{-j2\pi kn/(N/4)} + e^{-j2\pi k/N}\sum_{n=0}^{N/8-1} X[8n+1]e^{-j2\pi kn/(N/4)} + e^{-j2\pi (k+1/2)/N}\sum_{n=0}^{N/8-1} X[8n+5]e^{-j2\pi kn/(N/4)} + e^{-j2\pi (k+1/4)/N}\sum_{n=0}^{N/8-1} X[8n+3]e^{-j2\pi kn/(N/4)} + e^{-j2\pi (k+3/4)/N}\sum_{n=0}^{N/8-1} X[8n+7]e^{-j2\pi kn/(N/4)}
后续以此规律进行基二分法迭代运算即可
FFT举例验证
Eg2: 假设输入信号为F(t)=cos(2Π3t)，以采样频率8Hz采集8个点，得到以下波形数据：
X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)
0.0 0.71 -1 0.71 0.0 -0.71 1 -0.71
基2蝶式FFT运算过程如下所示：
首先我们需要将输入序列不断进行奇偶提取形成奇偶队列，以此不断迭代直到最后每个奇偶序列只有一个数据，如下所示：
未分解 x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7)
第一次分解 x(0) x(2) x(4) x(6) x(1) x(3) x(5) x(7)
第二次分解 x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)
注: 红色--本轮的偶序列、紫色--为本轮的奇序列。
蝶形运算示意图如下所示：
由DFT运算过程我们可知：
我们以X(3)作为示例展示FFT运算过程：
由以上两计算式可以看出，根据旋转因子的周期性、缩放性以及对称性，我们亦可由DFT结果推导得出FFT计算后的X(3)的计算式。无论是DFT还是FFT计算得出的X(3)的结果是一致的。
注： 根据以上FFT运算蝶形图可知，完成傅里叶变换需进行NlogN次复数的乘法运算，在进行NlogN次复数的加法运算。综合来看，FFT的时间复杂度近似为O(NlogN)，回顾前面DFT的时间复杂度近似于O(N²)。假如N = 1024，FFT时间负责度为1024log2(1024) = 10240≈O(10^4)，DFT的时间复杂度为1024² = 1048576≈O(10^6)。相比之下降低了两个数量级。