快速傅氏变换之旅(一) 概念简介 3

原创 2012年03月31日 17:16:43

 

1)

蝶形变换普通的FFT算法称为基2的FFT算法,这种算法的核心是蝶形变换  长度为n=2^k1的变换共需要做   k1   *   n/2   次蝶形变换,(如上图所示)若需变换数据表示为一个复数数组c[],则每次蝶形变换有2个输入   c[i],c[i+s],两个输出:c[i],c[i+s],s成为翅间距。   每个变换的基本算法是:
   
    t=wr   *   c[i+s];  
    c[i+s]=c[i]-t;
    c[i]=c[i]+t;

      前面说过,长度为n=2^k1的变换共需要做   k1   *   n/2次变换,实际的程序是一个3层循环,共需要k1*k2*(k3/2)次变换(k2*k3/2=n/2)。前面的wr是w的整数次方,w=e^(-   2*PI/k3)   (k3=2,4,8,16...n,PI是圆周率),也成为旋转因子,例如n=32的变换需要log2(32)=5趟变换:
      第1趟变换需要16*1次变换,翅间距是1,     若w=e^(-2*PI/2),   则wr=w^1
      第2趟变换需要8*2次变换,   翅间距是2,     若w=e^(-2*PI/4),   则wr=w^1,w^2
      第3趟变换需要4*2次变换,   翅间距是4,     若w=e^(-2*PI/8),   则wr=w^1,w^2,w^3,w^4  
      第4趟变换需要2*8次变换,   翅间距是8,     若w=e^(-2*PI/16),则wr=w^1,w^2,w^3,w^4,w^5,w^6,w^7,w^8
      第5趟变换需要16*1次变换,翅间距是16,   若w=e^(-2*PI/32),则wr=w^1,w^2,w^3,w^4,w^5...w^15,w^16

void  fft(struct complex c_IO[],int m)
{
	int L, B, j, k,  p, q;
	for(k = 0; k < SAMPLECOUNT; k++)
	{
	  g_out[k] = c_IO[k];
	}
	for(L = 1.0; L <= m; L++)
	{
		B=(int)pow(2.0,	L-1);
		for(j=0;j<=B-1;j++)
		{
			p=j*(int)pow(2.0, m-L);
			q=(int)pow(2.0, L);
			for(k=j;k<=SAMPLECOUNT-1;k=k+q)
			{
				temp[k]=complex_add(g_out[k],complex_mult(W_n_k(SAMPLECOUNT,p) ,g_out[k+B] ) );
				temp[k+B]=complex_remove(g_out[k],complex_mult(W_n_k(SAMPLECOUNT,p) ,g_out[k+B] ) );
				g_out[k]=temp[k];
				g_out[k+B]=temp[k+B];
			}
		}
	}
}

 

//*************************************************************************
// *
// * 函数名称:
// *   FFT()
// *
// * 参数:
// *   complex<double> * TD- 指向时域数组的指针
// *   complex<double> * FD- 指向频域数组的指针
// *   r-2的幂数,即迭代次数
// *
// * 返回值:
// *   无。
// *
// * 说明:
// *   该函数用来实现快速付立叶变换。
// *
// ************************************************************************/
void FFT(	complex<double> *TD,  complex<double> *FD, 
			complex<double> *X1,  complex<double> *X2,
			complex<double> *Omega,	  int r)
{
	// 付立叶变换点数
	long count;

	// 循环变量
	int i,j,k;

	// 中间变量
	int bfsize,p;

	// 角度
	double angle;

	complex<double> *X;

	// 计算付立叶变换点数
	count = 1 << r;

	// 分配运算所需存储器
	//Omega  = new complex<double>[count / 2];
	//X1 = new complex<double>[count];
	//X2 = new complex<double>[count];

	// 计算加权系数
	for(i = 0; i < count / 2; i++)
	{
		angle = -i * 3.1415926 * 2 / count;
		Omega[i] = complex<double>(cos(angle), sin(angle));
	}

	// 将时域点写入X1
	memcpy(X1, TD, sizeof(complex<double>) * count);

	// 采用蝶形算法进行快速付立叶变换
	for(k = 0; k < r; k++)
	{
		for(j = 0; j < 1 << k; j++)
		{
			bfsize = 1 << (r-k);
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * Omega[i * (1<<k)];
			}
		}
		X  = X1;
		X1 = X2;
		X2 = X;
	}

	// 重新排序
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j&(1<<i))
			{
				p+=1<<(r-i-1);
			}
		}
		FD[j]=X1[p];
	}

	// 释放内存
	//delete Omega;
	//delete X1;
	//delete X2;
}



2)复数数组排序,在基2的蝶形变换中,复数数组需要重新排序,c[i]   要放置到数组c的第   reverse(c[i])   的位置,m=reverse(n)   函数的算法是这样的:

若   n的   k位2进制的为b[],   b[k-1],B[k-2],...b[2],b[1],b[0],(   b[i]   等于1或者0,b[0]为最低bit).   则m=reverse(n)的2进制的为   b[0],b[1],b[2],b[3],...b[k-1]   (b[k-1]为最低bit).

2.更复杂的变换算法:基2的蝶形变换算法不止一种,它可分为2类,一类为基2时分傅立叶变换,另一类为基2频分傅立叶变换上例的变为基2时分算法,在每一趟变换中,翅间距依次变大,第一趟为2,最后一趟为n/2,数组重排在变换之前进行,基2频分算法正好相反,翅间距依次缩小,起于n/2,止于2,数组重排在蝶形变换之后进行。   在<傅立叶变换>一书中,提到3种基2时分变换,3种基2频分变换。上述算法称为基2时分FFT第二种算法。我在看你的这个程序的同时,还看到朱志刚写的一个FFT程序,这个程序的算法是基2时分FFT第一种算法,它比经典的算法更复杂,需要对wr数组进行逆序排列。

3.更复杂的FFT算法,除了基2   的FFT算法外,还有更加复杂的基4算法,基8算法,甚至基3,基5算法,纯粹的基4算法只能计算长度为4^k的变换,但它比基2的算法速度更高。为了提高速度,很多FFT算法使用混合基算法。如我看到的2个效率很高程序均使用了混合基算法。第一个程序是来自:http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html,它使用了基2,基4,甚至基8混合算法,共提供了6同样功能的算法。可对长度为2^k的序列做FFT,这个程序的写的很复杂,我现在尚不能完全看懂。另一个程序来自:http://hjem.get2net.dk/jjn/fft.htm。相对于前者,这个程序相对简单一点。它使用了基2,基3,基4,基5,基8,基10   混合算法,几乎可以计算任意长度的FFT。具体的,当序列长度n为2,3,5,7,11,13,17,19,23,29,31,37等小素数时,或者n的最大素因数小于等于37时,可计算这个序列的FFT。

相关文章推荐

快速傅氏变换之旅(一) 概念简介 1

FFT(Fast Fourier Transformation),即为快速傅氏变换,是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。它对傅氏...

快速傅氏变换之旅(五) 嵌入式中的FFT(最好选择DSP或FPGA)

欢迎大家提意见讨论。 转载请标明是引用于 http://blog.csdn.net/chenyujing1234  例子代码:(编译工具:VS2005)     ==============...

快速傅氏变换之旅(四) 第一个别人的例子解读

转载请标明是引用于 http://blog.csdn.net/chenyujing1234  例子代码:(编译工具:VS2005)  http://www.rayfile.com/zh-cn/fi...

快速傅氏变换

  • 2012年09月15日 10:21
  • 2KB
  • 下载

【快速学习】学习cocos2dx之旅(3)

3.创建自己的hellowordScene 右键添加 选择新建项 自己创建一个头文件和一个cpp文件 头文件中先打出如上图代码,这时防止重复包含头文件的过程 紧接着包含cocos2d...

八位大牛带你见证H3CIE认证之旅

  • 2014年11月24日 16:02
  • 1.02MB
  • 下载

cocos2dx3.x游戏开发之旅1-17章源代码

  • 2017年08月04日 09:34
  • 33.27MB
  • 下载

JAVA进阶之旅(二)——认识Class类,反射的概念,Constructor,Field,Method,反射Main方法,数组的反射和实践

JAVA进阶之旅(二)——认识Class类,反射的概念,Constructor,Fiald,Method,反射Main方法,数组的反射和实践 我们继续聊JAVA,这次比较有意思,那就是反射了 一....

Cocos2d-x 3.x游戏开发之旅

  • 2017年06月01日 08:39
  • 84.78MB
  • 下载
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:快速傅氏变换之旅(一) 概念简介 3
举报原因:
原因补充:

(最多只允许输入30个字)