详细的频域滤波学习笔记（4）--基2时间抽取FFT

  对于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（讲的非常仔细地一个视频）