详细的频域滤波学习笔记（4）–基2时间抽取FFT算法August的博客-

  对于N个数，直接进行DFT计算时需要N²次复数乘法，以及N（N-1）次复数加法，这使得计算效率非常低下，这里介绍一种提高计算效率的算法，即基于时间抽取的快速傅里叶变换。
1.算法原理
  设输入序列长度为N=2^M（M为整数），将该序列按照时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也成为Coolkey-Tukey算法。其中即2表示：N=2^M,M为整数。若不满足这个条件，可以认为地加上若干零值（加零补长）使其达到N=2^m。

2.算法步骤
①分组，变量置换

  先将x(n)按n地奇偶分为两组，变量置换：
     当n=偶数时，令n=2r;
     当n=奇数时，令n=2r+1;
得到

②将上式代入DFT中，有以下结果：

由于旋转因子地性质

可以对X（k）进行变形，可以得到

X₁(k),X₂(k)只有N/2个点，以N/2为周期；而X(k)却有N个点，以N为周期。要用X₁(k),X₂(k)表达全部的X(k)的值，还必须要利用旋转因子的周期性：


又根据旋转因子的对称性：

因此可以得出下式：

该式子称之为蝶形运算，蝶形运算流图符号如下所示：

对蝶形运算进行说明：（1）左边两路为输出；（2）右边两路为输出；（3）中间以一个小圆表示加、减运算（右上为相加输出，右下为相减输出）。一次蝶形运算需要1次复乘，2次复加。
采用蝶形运算可以大大减少运算量，有下表做出比较：

例子：求N=2³=8点FFT变换
  按照基2时间抽取FFT算法，先将N=8点的DFT分解成2个4点DFT，再将4点的DFT分解成2个2点的DFT，再继续分解。
  首先分解为两个4点的DFT，有如下分解图：

再分解为2点的DFT，有如下示意图：

最后分解为1点的DFT，有如下示意图：

以上就是整个的过程的推导，在实际代码实现中需要注意第一级蝶形运算时,x(k)打乱了顺序，可以注意到x(1)与x(4)，x(3)与x(6)进行了替换，不难发现1的二进制与4的二进制数刚好倒序，因此在进行DFT前需要对原先的数组进行位逆序置换，这里贴上整个基2时间抽取的代码：

//位逆序置换 complex<double>*Rev(complex<double>*Data, int length) { int a = 1, num = 0; while (a <= length) {    a <<= 1;    num++; }   vector<int>temp(length); for (int i = 0; i < length; i++) {    temp[i] = 0; } for (int i = 0; i<length; i++) { int m = i; for (int j = 0; j < num - 1; j++) {          temp[i] <<= 1;     temp[i] |= (m & 1);     m >>= 1; } if (i < temp[i]) {     complex<double>  b = Data[i];     Data[i] = Data[temp[i]];     Data[temp[i]] = b; } } return Data; } //基2时间抽取FFT，sign=-1时为正变换，sign=1时为逆变换 void dit2(complex<double>*Data, int Log2N, int sign) { int i, j, k, step, length;   complex<double> wn, temp, deltawn;   length = 1 << Log2N;   complex<double>*Data = Rev(Data, length); for (i = 0; i<length; i += 2) {    temp = Data[i];    Data[i] = Data[i] + Data[i + 1];    Data[i + 1] = temp - Data[i + 1]; } for (i = 2; i <= Log2N; i++) {    wn = 1; step = 1 << i; deltawn = complex<double>(cos(2.0*PI / step), sin(sign*2.0*PI / step)); for (j = 0; j<step / 2; j++) { for (i = 0; i<length / step; i++) {      temp = Data[i*step + step / 2 + j] * wn;      Data[i*step + step / 2 + j] = Data[i*step + j] - temp;      Data[i*step + j] = Data[i*step + j] + temp; }     wn = wn*deltawn; } } if (sign == 1) for (i = 0; i<length; i++)     Data[i] /= length; } 

参考资料：https://www.bilibili.com/video/BV1W7411c7Kc?p=2（讲的非常仔细地一个视频）