一、FFT介绍
傅里叶变换是数字信号处理领域一个很重要的数学变换，它用来实现将信号从时域到频域的变换，在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域有广泛的应用。离散傅里叶变换(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）
快速傅里叶变换（英语：Fast Fourier Transform, FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法[1]。傅里叶分析将信号从原始域（通常是时间或空间）转换到頻域的表示或者逆过来转换。
FFT会通过把DFT矩阵分解为稀疏（大多为零）因子之积来快速计算此类变换。因此，它能够将计算DFT的复杂度从只用DFT定义计算需要的 O（n2），降低到O（nlog n），其中n为数据大小。
快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及，但早在1805年就已推导出来。1994年美國數學家吉尔伯特·斯特朗把FFT描述为“我们一生中最重要的数值算法”，它还被IEEE科学与工程计算期刊列入20世纪十大算法。
三、FFT的CPU实现
一维FFT基2算法的实现
我们使用按频率抽取的方法实现了一维FFT基2算法。算法的关键代码如下：
(1) 声明双精度复数的结构体：
struct Complex
{
double re; //复数的实部
double im; //复数的虚部
};
(2) 通过幂数获得快速傅里叶变换的长度，并初始化：
count = 1 << power; //power为幂数，count为快速傅里叶变换的长度
a = new Complex[count]; //a为初始化的数组，用来存放时域复数的值
b = new Complex[count]; //b为变换后存放结果的数组
memcpy(a, t, sizeof(Complex)*count); //初始化，将时域数据存放在a中，t为时域数据
(3) 计算旋转因子：
w = new Complex[count / 2];
for (i = 0; i<count / 2; i++)
{
angle = -i*pi * 2 / count;
w[i].re = cos(angle);
w[i].im = sin(angle);
}
(4) 采用频率分解法进行蝶形运算：
for (k = 0; k<power; k++)
{
for (j = 0; j<1 << k; j++)
{
bfsize = 1 << (power - k);
for (i = 0; i<bfsize / 2; i++)
{
p = j*bfsize;
b[i + p] = Add(a[i + p], a[i + p + bfsize / 2]); //Add()函数实现复数的加法
b[i + p + bfsize / 2] = Mul(Sub(a[i + p], a[i + p + bfsize / 2]), w[i*(1 << k)]);//Mul()函数实现复数的乘法，Sub()函数实现复数的减法
}
}
x = a;
a = b;
b = x;
}
(5) 蝶形运算全部完成后，对结果进行整序，从而得到正确的输出顺序：
for (j = 0; j<count; j++)
{
p = 0;
for (i = 0; i<power; i++)
{
if (j&(1 << i))
p += 1 << (power - i - 1);
}
f[j] = a[p];
}
二维FFT基2算法的实现
通过两次调用一维快速傅里叶变换即可实现二维快速傅里叶变换。关键代码如下：
(1) 计算进行傅里叶变换的宽度、高度，以及垂直方向上、水平方向上的迭代次数：
while (w * 2 <= width) //计算进行傅立叶变换的宽度（2的整数次方）
{
w *= 2; //w为变换的宽度
wp++; //wp为垂直方向上的迭代次数
}
while (h * 2 <= height)//计算进行傅立叶变换的高度（2的整数次方）
{
h *= 2; //h为变换的高度
hp++; //hp为水平方向上的迭代次数
}
(2) 两次调用一维快速傅里叶变换：
for (j = 0; j<h; j++) //在垂直方向上进行快速傅立叶变换
{
FFT1D(&t[w*j], &f[w*j], wp);
}
for (j = 0; j<h; j++) //转换变换结果
{
for (i = 0; i<w; i++)
{
t[j + h*i] = f[i + w*j];
}
}
for (j = 0; j<w; j++) //水平方向进行快速傅立叶变换
{
FFT1D(&t[j*h], &f[j*h], hp);
}
四、FFT的GPU实现
对一维或多维信号进行离散傅里叶变换的FFT变换自身具有可“分治”实现的特点，因此能高效地在GPU平台上实现。CUDA提供了封装好的CUFFT库，它提供了与CPU上的FFTW库相似的接口，能够让使用者轻易地挖掘GPU的强大浮点处理能力，又不用自己去实现专门的FFT内核函数。使用者通过调用CUFFT库的API函数，即可完成FFT变换。
常见的FFT库在功能上有很多不同。有些库采用了基2变换，只能处理长度为2的指数的FFT，而其他的一些可以处理任意长度的FFT。CUFFT4.0提供了以下功能：
(1) 对实数或复数进行一维、二维或三维离散傅里叶变换；
(2) 可以同时并行处理一批一维离散傅里叶变换；
(3) 对任意维的离散傅里叶变换，单精度最大长度可达到6400万，双精度最大长度可达到12800万（实际长度受限于GPU存储器的可用大小）；
(4) 对实数或复数进行的FFT，结果输出位置可以和输入位置相同（原地变换），也可以不同；
(5) 可在兼容硬件（GT200以及之后的GPU）上运行双精度的变换；
(6)支持流执行：数据传输过程中可以同时执行计算过程。
一维FFT算法的CUDA实现
关键代码如下：
#include <cufft.h> //CUFFT文件头
#define NX 1024
#define BATCH 1
cufftDoubleComplex *data; //显存数据指针
//在显存中分配空间
cudaMalloc((void**)&data, sizeof(cufftDoubleComplex)*NX*BATCH);
//创建CUFFT句柄
cufftHandle plan;
cufftPlan1d(&plan, NX, CUFFT_Z2Z, BATCH);
//执行CUFFT
cufftExecZ2Z(plan, data, data, CUFFT_FORWARD); //快速傅里叶正变换
//销毁句柄，并释放空间
cufftDestroy(plan);
cudaFree(data);
二维FFT算法的CUDA实现
二维FFT算法实现中，同一维FFT不同的是：
(1) 输入参数：没有BATCH，增加了NY。NX为行数，NY为列数；
(2) FFT的正变换的输入位置和输出位置不同。
关键代码如下：
#include <cufft.h> //CUFFT文件头
#define NX 1024
#define NY 1024
cufftDoubleComplex *idata, *odata; //显存数据指针
//在显存中分配空间
cudaMalloc((void**)&idata, sizeof(cufftDoubleComplex)*NX*NY);
cudaMalloc((void**)&odata, sizeof(cufftDoubleComplex)*NX*NY);
//创建CUFFT句柄
cufftHandle plan;
cufftPlan2d(&plan, NX, NY, CUFFT_Z2Z);
//执行CUFFT
cufftExecZ2Z(plan, idata, odata, CUFFT_FORWARD); //快速傅里叶正变换
//销毁句柄，并释放空间
cufftDestroy(plan);
cudaFree(idata);
cudaFree(odata);
三维FFT算法的CUDA实现
#define NX 64 #define NY 64 #define NZ 128 cufftHandle plan; cufftComplex *data1, *data2; cudaMalloc((void**)&data1, sizeof(cufftComplex)*NX*NY*NZ); cudaMalloc((void**)&data2, sizeof(cufftComplex)*NX*NY*NZ); /* Create a 3D FFT plan. */ cufftPlan3d(&plan, NX, NY, NZ, CUFFT_C2C); /* Transform the first signal in place. */ cufftExecC2C(plan, data1, data1, CUFFT_FORWARD); /* Transform the second signal using the same plan. */ cufftExecC2C(plan, data2, data2, CUFFT_FORWARD); /* Destroy the cuFFT plan. */ cufftDestroy(plan); cudaFree(data1); cudaFree(data2);
五、实例
#include <assert.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <math.h> // Include CUDA runtime and CUFFT #include <cuda_runtime.h> #include <cufft.h> // Helper functions for CUDA #include "device_launch_parameters.h" #define pi 3.1415926535 #define LENGTH 100 //signal sampling points int main() { // data gen float Data[LENGTH] = { 1,2,3,4 }; float fs = 1000000.000;//sampling frequency float f0 = 200000.00;// signal frequency for (int i = 0; i < LENGTH; i++) { Data[i] = 1.35*cos(2 * pi*f0*i / fs);//signal gen, } cufftComplex *CompData = (cufftComplex*)malloc(LENGTH * sizeof(cufftComplex));//allocate memory for the data in host int i; for (i = 0; i < LENGTH; i++) { CompData[i].x = Data[i]; CompData[i].y = 0; } cufftComplex *d_fftData; cudaMalloc((void**)&d_fftData, LENGTH * sizeof(cufftComplex));// allocate memory for the data in device cudaMemcpy(d_fftData, CompData, LENGTH * sizeof(cufftComplex), cudaMemcpyHostToDevice);// copy data from host to device cufftHandle plan;// cuda library function handle cufftPlan1d(&plan, LENGTH, CUFFT_C2C, 1);//declaration cufftExecC2C(plan, (cufftComplex*)d_fftData, (cufftComplex*)d_fftData, CUFFT_FORWARD);//execute cudaDeviceSynchronize();//wait to be done cudaMemcpy(CompData, d_fftData, LENGTH * sizeof(cufftComplex), cudaMemcpyDeviceToHost);// copy the result from device to host for (i = 0; i < LENGTH / 2; i++) { printf("i=%d f= %6.1fHz RealAmp=%3.1f ", i, fs*i / LENGTH, CompData[i].x*2.0 / LENGTH); printf("ImagAmp=+%3.1fi", CompData[i].y*2.0 / LENGTH); printf(" "); } cufftDestroy(plan); free(CompData); cudaFree(d_fftData); }