FFT算法原理（按时间抽取的基2算法）

离散傅里叶变换DFT的定义为：其中，x[n]是N点的采样序列，X[k]为N点的DFT计算结果，，而称为旋转因子。假设x[n]的采样频率为fs，那么得出的X[k]的模值和幅角可以代表 f = k*fs/N处的信号强度（非0Hz处，为幅值的N/2倍，0Hz处为幅值的N倍）和相位信息。具体的频率分量和它的幅值Am之间的关系为：可以看到，离散傅里叶变换的分辨率为fs/N，采样获得的频谱范围为可以看到DFT

杰洛米

3498人浏览 · 2023-11-30 20:28:58

杰洛米 · 2023-11-30 20:28:58 发布

预备知识

离散傅里叶变换DFT的定义为：

$X\left[ k \right] = \sum\limits_{n = 0}^{N - 1} {x\left[ n \right]} \cdot {e^{ - j\frac{{2\pi }}{N}nk}} = \sum\limits_{n = 0}^{N - 1} {x\left[ n \right]} \cdot W_N^{nk},k = 0,1,...,N - 1$

其中，x[n]是N点的采样序列，X[k]为N点的DFT计算结果， ${W_N} = {e^{ - j\frac{{2\pi }}{N}}}$ ，而 $W_N^{nk} = {e^{ - j\frac{{2\pi }}{N}nk}}$ 称为旋转因子。

假设x[n]的采样频率为fs，那么得出的X[k]的模值和幅角可以代表 f = k*fs/N处的信号强度（非0Hz处，为幅值的N/2倍，0Hz处为幅值的N倍）和相位信息。具体的频率分量和它的幅值Am之间的关系为：

$X\left[ k \right] = \left\{ \begin{array}{l} N \cdot {A_m}\left( {f = 0} \right)\\ \frac{N}{2} \cdot {A_m}\left( {f = k \cdot \frac{{{f_{\rm{s}}}}}{N}} \right),k = 1,2,...,\frac{N}{2} - 1 \end{array} \right.$

可以看到，离散傅里叶变换的分辨率为fs/N，采样获得的频谱范围为

$f = 0,\frac{{{f_{\rm{s}}}}}{N},2 \cdot \frac{{{f_{\rm{s}}}}}{N},...,\left( {\frac{N}{2} - 1} \right) \cdot \frac{{{f_{\rm{s}}}}}{N}$

可以看到DFT无法获得高于一半采样频率以上的分量信息。如果你需要获取50Hz整数倍的频率分量，那就需要让50Hz为fs/N或者其整数倍。

例如基波频率为50Hz，N为128点，那么可以设置fs=6400Hz。另一方面，如果不是采样整数个周期的信号，则必然发生频谱泄露。

顺带一提，X[k]的幅角是采样信号在cos写法下的初相位，例如sin(2πf*t+30°)对应的分量的相角就是-60°，而不是30度。

总结：

X[k]的模值在k≠0时，等于f = k*fs/N对应信号幅值的N/2倍，在k=0时，代表直流分量幅值的N倍。

X[k]的幅角代表对应f = k*fs/N频率分量在cos写法下的初相角。

FFT算法原理（按时间抽取的基2算法）

在计算机中实现DFT需要非常大量的运算，所以我们一般不会直接使用DFT，而是使用FFT。FFT并不是近似计算，而是使用了一些数学技巧完成DFT，其计算结果和DFT完全一致，但是可以大大减少计算量。实际上，DFT可以分解为奇序列o[n]和偶序列e[n]的DFT变换：

$\begin{array}{l} X\left[ k \right] = \sum\limits_{n = 0}^{N - 1} {x\left[ n \right]} \cdot W_N^{nk}\\ = \sum\limits_{n = 0}^{N/2 - 1} {x\left[ {2n} \right]} \cdot W_N^{2nk} + \sum\limits_{n = 0}^{N/2 - 1} {x\left[ {2n + 1} \right]} \cdot W_N^{\left( {2n + 1} \right)k}\\ = \sum\limits_{n = 0}^{N/2 - 1} {e\left[ n \right]} \cdot W_N^{2nk} + W_N^k\sum\limits_{l = 0}^{N/2 - 1} {o\left[ n \right]} \cdot W_N^{2nk} \end{array}$

因为 $W_N^{2nk} = W_{N/2}^{nk}$ ，因此有

$\begin{array}{l} X\left[ k \right] = \sum\limits_{n = 0}^{N/2 - 1} {e\left[ n \right]} \cdot W_{N/2}^{nk} + W_N^k\sum\limits_{l = 0}^{N/2 - 1} {o\left[ n \right]} \cdot W_{N/2}^{nk}\\ {\rm{ }} = O\left[ k \right] + W_N^kE\left[ k \right],k = 0,...,\frac{N}{2} - 1 \end{array}$

也就是说，做N点的DFT可以通过将输入序列分为奇序列和偶序列，然后做2次N/2点的DFT即可。但是这样，只能得到X[k]的前半部分，后半部分需要使用到对称性 $W_N^{k + \frac{N}{2}} = - W_N^k$ ：

$\begin{array}{l} X\left[ {k + \frac{N}{2}} \right] = \sum\limits_{n = 0}^{N/2 - 1} {e\left[ n \right]} \cdot W_{N/2}^{n\left( {k + \frac{N}{2}} \right)} + W_N^{k + \frac{N}{2}}\sum\limits_{l = 0}^{N/2 - 1} {o\left[ n \right]} \cdot W_{N/2}^{n\left( {k + \frac{N}{2}} \right)}\\ {\rm{ }} = O\left[ k \right] - W_N^kE\left[ k \right],k = 0,...,\frac{N}{2} - 1 \end{array}$

其中用到了DFT变换结果的周期性 $O\left[ {k + \frac{N}{2}} \right] = O\left[ k \right],E\left[ {k + \frac{N}{2}} \right] = E\left[ k \right]$ 。那么DFT的结果就可以总结为：

$\left\{ \begin{array}{l} X\left[ k \right] = O\left[ k \right] + W_N^kE\left[ k \right]\\ X\left[ {k + \frac{N}{2}} \right] = O\left[ k \right] - W_N^kE\left[ k \right] \end{array} \right.,k = 0,...,\frac{N}{2} - 1$

因此，一个N点的DFT可以分解为两个N/2点的DFT然后再对结果进行组合；同理，O[k]和E[k]也可以一直往前迭代，直到将N点DFT分解为N/2次的2点DFT为止，这就是FFT算法的思路。显然，要做FFT就必须保证N是2的整数幂。

实际上，由于FFT算法的原理一般是从输出到输入的表述，即使你清楚了这个过程对你的编程也没什么帮助（-_-||），因此对于编程实现最有帮助的实际上是观察蝶形图，总结规律。

位倒序

FFT算法的过程包括位倒序和蝶形运算。蝶形图的输入端并不是直接输入采样序列x[n]，而是需要对x[n]进行重排，也就是位倒序。

记N = {2^M}，则位倒序的规则是：将数组下标0,1,…,N-1转换为M位的2进制，再将二进制下标倒序，如8位的位倒序过程就是：

000, 001, 010, 011, 100, 101, 110, 111 → 000, 100, 010, 110, 001, 101, 011, 111

0, 1, 2, 3, 4, 5, 6, 7 → 0, 4, 2, 6, 1, 5, 3, 7

这就是位倒序的过程。因此原有的序列将被排序为：

x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]

为了节省RAM空间，可以在程序里直接对原有的x[n]进行覆盖；对于蝶形运算，也可以进行一样的操作。

更进一步来说，实际上可以直接写好固定点数对应的位倒序下标数组，在存储FFT的输入数据时，直接按照位倒序的顺序x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]存入，这样可以节省FFT的运行时间。

蝶形运算

我们定义蝶形运算为：

那么根据之前的DFT向前的二分解原则，可以画出8点FFT的蝶形图为：

观察蝶形图，可以总结到以下规律（记 $N = {2^M}$ ）：

蝶形运算共有M层。
蝶形运算的组数逐层递减，第L层包含 $GroupNum = {2^{M - L}}$ 组蝶形运算（不同组之间蝶形不交叉），第一层包含N/2组蝶形运算，之后每一层减半。
蝶形运算的两个输入数据间隔和旋转因子 $W_N^k$ 数量相同, $FactorNum = {2^{L - 1}}$ ，都是逐层递增，第一层为1，之后每一层翻倍，最后一层为N/2。
每一层的旋转因子Factor都是从1开始递增，每一次的增量为 $W_N^{GroupNum}$ ，即旋转因子需要乘以这个数，称之为旋转因子倍率。

编程思路应为三层循环嵌套（伪代码如下）：

Factor_Num = 1;
Group_Num = N/2;
for (L = 0, L < M, L++)
{
      Factor = 1；
      计算当前层的旋转因子倍率；
      for (W = 0, W < W_Num, W++)
      {
            for (G = 0, G < Group_Num, G++)
            {
                  计算蝶形运算的对应下标；
                  更新下标对应的对应的两个数组值；
            }
            Factor=Factor乘以对应倍率;
      }
      Factor_Num = Factor_Num*2;
      Group_Num = Group_Num/2;
}