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

   日期:2020-05-30     浏览:217    评论:0    
核心提示:  对于N个数,直接进行DFT计算时需要N2次复数乘法,以及N(N-1)次复数加法,这使得计算效率非常低下,这里介绍一种提高计算效率的算法,即基于时间抽取的快速傅里叶变换。1.算法原理  设输入序列长度为N=2M(M为整数),将该序列按照时间顺序的奇偶分解为越来越短的子序列,称为基2按时间抽取的FFT算法。也成为Coolkey-Tukey算法。其中即2表示:N=2M,M为整数。若不满足这个条件,可以认为地加上若干零值(加零补长)使其达到N=2m。2.算法步骤①分组,变量置换  先将x(n)按n地

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

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

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

②将上式代入DFT中,有以下结果:

由于旋转因子地性质

可以对X(k)进行变形,可以得到

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


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

因此可以得出下式:

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

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

例子:求N=23=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(讲的非常仔细地一个视频)

 
打赏
 本文转载自:网络 
所有权利归属于原作者,如文章来源标示错误或侵犯了您的权利请联系微信13520258486
更多>最近资讯中心
更多>最新资讯中心
更多>相关资讯中心
0相关评论

推荐图文
推荐资讯中心
点击排行
最新信息
新手指南
采购商服务
供应商服务
交易安全
关注我们
手机网站:
新浪微博:
微信关注:

13520258486

周一至周五 9:00-18:00
(其他时间联系在线客服)

24小时在线客服