cuFFT作为CUDA最基础的库之一是NVIDIA提供的GPU加速的Fourier变换FFT库能极大提升涉及FFT计算的科学计算、信号处理和深度学习等任务的速度。一、傅里叶变换和快速傅里叶变换​ Fourier变换是数字信号处理领域一个很重要的数学变换它用来实现将信号实现将信号从时域到频域的变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域有广泛的应用。离散傅里叶变换(Discrete Fourier TransformDFT)是连续傅里叶变换在离散系统中的表示形式由于DFT的计算量很大因此在很长一段时间内其应用受到了很大的限制。快速傅里叶变换(Fast Fourier TransformFFT)算法由Cooley和Tukey在20世纪60年代1965年提出它是DFT的快速算法使得离散傅里叶变换和卷积这类难度很大的计算工作的复杂度从N2量级降到了Nlog2N量级大大提高了DFT的运算速度从而使DFT在实际应用中得到了广泛的应用。FFT具体流程单个DFT分成了看起来相似两个更小的DFT。一个奇一个偶。从上面的对称性特点来看我们只需对每个子问题作一半的计算量只要我们的小傅里叶变换是偶倍数就可以再作分治直到分解出来的子问题小到无法通过分治提高效率接近极限时这个递归是 O(n log2n) 级的。这个递归算法能在python里快速实现当子问题被分解到合适大小时再用回原本那种“慢方法”。二、FFT在CUDA中的实现对一维或多维信号进行离散傅里叶变换的FFT变换自身具有可“分治”实现的特点因此能高效地在GPU平台上实现。CUDA提供了封装好的cuFFT库能够让使用者轻易地挖掘GPU的强大浮点处理能力又不用自己去实现专门的FFT内核函数。使用者通过调用CUFFT库的API函数即可完成FFT变换以下是一套基本调用流程在主机端准备输入数据在GPU设备端上分配内存按需求提取数据或将数据从主机复制到设备创建一个句柄plan, 调用函数cufftPlane1D/cufftPlane2D/cufftPlan3D 可以创建一个简单的Fourier变换。调用函数cufftPlanMany 则可以创建支持更多配置操作的变换计划。执行plan。这一步可以使用cufftExecC2C()、cufftExecR2C()或cufftExecC2R()等函数完成上一步完成plane的计算任务。执行完成以下若不再需要该plan则调用cufftDestroy()函数销毁该plan 及为其分配的计算资源。//基本流程正变换 #include cufft.h //CUFFT文件头 #define NX 1024 //要执行FFT的信号长度 #define BATCH 1 //要执行FFT的信号个数 cufftDoubleComplex *data; //显存数据指针 //在显存中分配空间 cudaMalloc ((void**) data, sizeof (cufftDoubleComplex) *NX*BATCH) ; //创建CUFFT句柄 //cufftPlan1d()针对单个 1 维信号 //cufftPlan2d()针对单个 2 维信号 //cufftPlan3d()针对单个 3 维信号 //cufftPlanMany()针对多个信号同时进行 fft cufftHandle plan; cufftPlanld (plan, NX, CUFFT_C2C, BATCH) ; //(句柄名信号长度输入类型输出类型信号个数) //执行CUFET cufftExecC2C (plan,data,data, CUFFT_FORWARD); //(句柄名输入信号首地址输出信号首地址正/逆变换) //销毁句柄并释放空间 cufftDestroy(plan); cudaFree (data); //进阶流程首先生成一维的随机信号利用cufft 先进行正变换然后逆变换并判定逆变换后结果与原输入信号判断是否相等。 #include iostream #include time.h #include cuda_runtime.h #include device_launch_parameters.h #include cufft.h #define NX 3335 // 有效数据个数 #define N 5335 // 补0之后的数据长度 #define BATCH 1 #define BLOCK_SIZE 1024 using std::cout; using std::endl; /** * 功能判断两个 cufftComplex 数组的是否相等 * 输入idataA 输入数组A的头指针 * 输入idataB 输出数组B的头指针 * 输入size 数组的元素个数 * 返回true | false */ bool IsEqual(cufftComplex *idataA, cufftComplex *idataB, const int size) { for (int i 0; i size; i) { if (abs(idataA[i].x - idataB[i].x) 0.000001 || abs(idataA[i].y - idataB[i].y) 0.000001) return false; } return true; } /** * 功能实现 cufftComplex 数组的尺度缩放也就是乘以一个数 * 输入idata 输入数组的头指针 * 输出odata 输出数组的头指针 * 输入size 数组的元素个数 * 输入scale 缩放尺度 */ static __global__ void cufftComplexScale(cufftComplex *idata, cufftComplex *odata, const int size, float scale) { const int threadID blockIdx.x * blockDim.x threadIdx.x; if (threadID size) { odata[threadID].x idata[threadID].x * scale; odata[threadID].y idata[threadID].y * scale; } } int main() { cufftComplex *data_dev; // 设备端数据头指针 cufftComplex *data_Host (cufftComplex*)malloc(NX*BATCH * sizeof(cufftComplex)); // 主机端数据头指针 cufftComplex *resultFFT (cufftComplex*)malloc(N*BATCH * sizeof(cufftComplex)); // 正变换的结果 cufftComplex *resultIFFT (cufftComplex*)malloc(NX*BATCH * sizeof(cufftComplex)); // 先正变换后逆变换的结果 // 初始数据 for (int i 0; i NX; i) { data_Host[i].x float((rand() * rand()) % NX) / NX; data_Host[i].y float((rand() * rand()) % NX) / NX; } dim3 dimBlock(BLOCK_SIZE); // 线程块 dim3 dimGrid((NX BLOCK_SIZE - 1) / dimBlock.x); // 线程格 cufftHandle plan; // 创建cuFFT句柄 cufftPlan1d(plan, N, CUFFT_C2C, BATCH); // 计时 clock_t start, stop; double duration; start clock(); cudaMalloc((void**)data_dev, sizeof(cufftComplex)*N*BATCH); // 开辟设备内存 cudaMemset(data_dev, 0, sizeof(cufftComplex)*N*BATCH); // 初始为0 cudaMemcpy(data_dev, data_Host, NX * sizeof(cufftComplex), cudaMemcpyHostToDevice); // 从主机内存拷贝到设备内存 cufftExecC2C(plan, data_dev, data_dev, CUFFT_FORWARD); // 执行 cuFFT正变换 cudaMemcpy(resultFFT, data_dev, N * sizeof(cufftComplex), cudaMemcpyDeviceToHost); // 从设备内存拷贝到主机内存 cufftExecC2C(plan, data_dev, data_dev, CUFFT_INVERSE); // 执行 cuFFT逆变换 cufftComplexScale dimGrid, dimBlock (data_dev, data_dev, N, 1.0f / N); // 乘以系数 cudaMemcpy(resultIFFT, data_dev, NX * sizeof(cufftComplex), cudaMemcpyDeviceToHost); // 从设备内存拷贝到主机内存 stop clock(); duration (double)(stop - start) * 1000 / CLOCKS_PER_SEC; cout 时间为 duration ms endl; cufftDestroy(plan); // 销毁句柄 cudaFree(data_dev); // 释放空间 cout IsEqual(data_Host, resultIFFT, NX) endl; return 0; }