cuFFT的简介及实现案例

1. cuFFT库的简介Introduction of cuFFT libaray​ Fourier变换是数字信号处理领域一个很重要的数学变换它用来实现将信号实现将信号从时域到频域的变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域有广泛的应用。离散傅里叶变换(Discrete Fourier TransformDFT)是连续傅里叶变换在离散系统中的表示形式由于DFT的计算量很大因此在很长一段时间内其应用受到了很大的限制。20世纪60年代1965年由Cooley和Tukey提出了快速傅里叶变换(Fast Fourier TransformFFT)算法它是DFT的快速算法使得离散傅里叶变换和卷积这类难度很大的计算工作的复杂度从N2量级降到了Nlog2N量级大大提高了DFT的运算速度从而使DFT在实际应用中得到了广泛的应用。​ cuFFT是NVIDIA提供的GPU加速的Fourier变换FFT库能极大提升涉及FFT计算的科学计算、信号处理和深度学习等任务的速度。下表概括了器主要特征和应用场景cuFFT的特征具体描述基本功能提供GPU加速的1D、2D、3D复数/实数FFT计算核心优势相比CPU实现利用GPU并行性可获得显著加速编程接口提供类似的FFTW的API便于熟悉CPU FFT的用户迁移高级功能支持批量执行、流异步、半/单/双精度、多GPU计算主要应用领域深度学习、计算物理学、医学成像、信号处理等2. 基于cuFFT库的Fourier变换步骤workflow of Fourier Transform based cuFFT在CUDA上进行傅里叶变换一般需要做以下几步工作在主机端准备输入数据在GPU设备端上分配内存并将数据从主机复制到设备cudaMalloc,cudaMemcpy的接口 创建一个plan, 调用函数cufftPlane1D/cufftPlane2D/cufftPlan3D1/2/3 可以创建一个简单的Fourier变换。调用函数cufftPlanMany 则可以创建支持更多配置操作的变换计划。cufftPlan1d()1(): 针对单个1维信号cufftPlan2d()2():针对单个2维信号cufftPlan3d()3():针对单个3维信号执行plane。这一步可以使用cufftExecC2C()2()、cufftExecR2C()2()或cufftExecC2R()2()等函数完成上一步完成plane的计算任务。执行完成以下若不再需要该plan则调用cufftDestroy()()函数销毁该plan 及为其分配的计算资源。3. cuFFT的傅里叶变换API接口类型Fourier Transform Types​ cuFFT 库实现了三种不同类型的Fourier变换接口分为C2C2(复数变换到复数)C2R2(复数到实数), R2C2 (实数到复数)。本质上这三种转换都可以被看做是复数域到复数域的变换之所以这样划分其最主要的考量是性能因素。例如在一般的数字信号处理中输入数据是一些离散的实数域上的采样点这时候对它们做Fourier变换实际上就是R2C2根据埃尔米特对称性Hermitian symmetry)变换XkX∗N−k−∗, ∗∗ 代表共轭复数。cuFFT 的傅里叶变换则利用了这些冗余将计算量降到最低。​ 变换执行函数的单精度和双精度版本分别定义如下API描述cufftExecC2C()/cufftExecZ2Z()2()/2()单精度/双精度浮点数复数域到复数域的傅里叶变换cufftExecR2C()/cufftExecD2Z()2()/2()单精度/双精度浮点数实数域到复数域的傅里叶变换正变换cufftExecC2R()/cufftExecZ2D()2()/2()单精度/双精度浮点数复数域到实数域的傅里叶变换逆变换4. 数据布局(Data Layout)​ CUFFT库保含若干种数据类型对于复数有cufftComplex/cufftDoubleComplex/ 两种数据类型对于实数则分别有cufftReal/cufftDouble/ 两种数据类型 。​ 根据转换结果的存储位置不同FFT变换可分为就地变换(in−place−)和外部变换out−place−)前者直接在输入数据进行变换而后者则会将变换后结果存入新的存储器地址。​ 就地转换(in−place−) 支持数据的两种布局native 和 padded前者用于获取最佳性能而后者则用于与FFTW兼容。​ 在padded布局中输出信号的开始地址与输入信号一样换句话说实数域到复数域变换的输入数据和复数域到实数域的输出数据必须被填充。在native布局中则没有填入要求。​ 输入数据和输出数据的尺寸总结如下FFT typeinput data sizeoutput data sizeC2C2X cufftComplex X cufftComplex C2R2[X2]1 cufftComplex[2]1 X cufftRealR2C∗2∗X cufftReal[X2]1 cufftComplex[2]1 5. 单个一维信号的FFT变换代码实现One Dimension SIgnal FFT Transfrom在本次测试代码中首先生成一维的随机信号利用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; }其中cufftPlan1d():第一个参数就是要配置的cuFFT句柄第二个参数就是要进行fft的信号的长度第三个CUFFT_C2C为要执行fft 的信号输入类型及输出类型复数CUFFT_C2R表示输入复数输出实数CUFFT_R2C表示输入实数输出复数CUFFT_R2R表示输入实数输出实数第四个参数BATCH表示要执行fft的信号的个数新版的已经使用cufftPlanMany()来同时完成多个信号的fftcufftExecC2C():第一个参数就是配置好的 cuFFT 句柄第二个参数为输入信号的首地址第三个参数为输出信号的首地址第四个参数为CUFFT_FORWARD表示执行的是fft正变换CUFFT_INVERSE表示执行fft逆变换需要注意的是执行完fft之后要对信号中的每个值乘以1/N1/;输出结果