GPU-FFT总结.doc

引言 傅里叶变换是数字信号处理领域一个很重要的数学变换，它用来实现将信号从时域到频域的变换，在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域有广泛的应用。离散傅里叶变换(Discrete Fourier Transform，DFT)是连续傅里叶变换在离散系统中的表示形式，由于DFT的计算量很大，因此在很长一段时间内其应用受到了很大的限制。20世纪60年代（1965年）由Cooley和Tukey提出了快速傅里叶变换(Fast Fourier Transform，FFT)算法，它是DFT的快速算法，使得离散傅里叶变换和卷积这类难度很大的计算工作的复杂度从N2量级降到了Nlog2N量级，大大提高了DFT的运算速度，从而使DFT在实际应用中得到了广泛的应用。 传统上，GPU只负责图形渲染，而大部分的处理都交给了CPU。自二十世纪九十年代开始，GPU的发展迅速。由于GPU具有强大的并行计算能力，加之其可编程能力的不断提高，GPU也用于通用计算，为科学计算的应用提供了新的选择。 2007年6月，NVIDIA公司推出了CUDA (Compute Unified Device Architecture)，CUDA 不需要借助图形学API，而是采用了类C语言进行开发。同时，CUDA采用了统一处理架构，降低了编程的难度，同时，NVIDIA GPU引入了片内共享存储器，提高了效率。这两项改进使CUDA架构更加适合进行GPU通用计算。 快速傅里叶变换（FFT） 设x(n) 为长度为N的有限长序列，则其离散傅里叶正变换(DFT)为： 其中，，称为旋转因子或蝶形因子。 从DFT的定义可知，在为复数序列的情况下，对某个k值，直接按上式计算需要N次复数乘法和(N-1)次复数加法。因此，对所有N个 k值，共需要N2次复数乘法和N(N-1)次复数加法。对于N值来说，直接计算它的DFT所需的计算量是很大的。不难发现旋转因子具有对称性和周期性。而FFT的基本思想是利用的周期性和对称性，使长序列的DFT分解为更小长度的DFT，利用这些小的DFT的计算来代替大的DFT计算，从而达到提高效率的目的。 一般而言，FFT的分解法可归为两类：按时间抽取的FFT(DIT-FFT)和按频率抽取的FFT(DIF-FFT)。DIT-FFT算法是时域内将每一级输入序列依次按奇/偶分成两个短序列进行计算，而DIF-FFT算法是频域内将每一级输入序列按奇/偶分成两个短序列进行计算。两者的区别是旋转因子出现的位置不同，但算法是一样的。在DIT-FFT算法中，旋转因子出现在输入端，而DIF-FFT算法中它出现在输出端。 从FFT的发展看，主要有两个方向：一个是针对N等于2的基数次幂的算法，如基2算法、基4算法、实因子算法和分裂基算法等；另一个是N不等于2 的整数次幂的算法，它是以Winograd为代表的一类算法（素因子算法，Winograd算法）。 了解了一维快速傅里叶变换，那么就很容易实现二维快速傅里叶变换。二维快速傅里叶变换的输入相当于一个N行N列的复数矩阵，对其每行进行一维快速傅里叶变换，再对其每列进行一维快速傅里叶变换，这样就达到了变换的目的。二维快速傅里叶变换主要应用于对图像等的快速傅里叶变换。 FFT的CPU实现 一维FFT基2算法的实现 我们使用按频率抽取的方法实现了一维FFT基2算法。算法的关键代码如下： 声明双精度复数的结构体： struct Complex { double re; //复数的实部 double im; //复数的虚部 }; 通过幂数获得快速傅里叶变换的长度，并初始化： count = 1power; //power为幂数，count为快速傅里叶变换的长度 a = new Complex [count]; //a为初始化的数组，用来存放时域复数的值 b = new Complex [count]; //b为变换后存放结果的数组 memcpy(a,t,sizeof(Complex)*count); //初始化，将时域数据存放在a中，t为时域数据 计算旋转因子： w = new Complex [count/2]; for( i=0; icount/2; i++ ) { angle = -i*pi*2/count; w[i].re = cos(angle); w[i].im = sin(angle); } 采用频率分解法进行蝶形运算： for( k=0; kpower; k++ ) { for( j=0; j1k; j++ ) { bfsize = 1(power-k); for( i=0; ibfsize/2; i++ ) { p = j*b