NIVIDIA高性能计算CUDA笔记（三） cuFFT的简介及实现案例

NIVIDIA高性能计算CUDA笔记（三） cuFFT的简介及实现案例

1. cuFFT库的简介（Introduction of cuFFT libaray）

​ Fourier变换是数字信号处理领域一个很重要的数学变换，它用来实现将信号实现将信号从时域到频域的变换，在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域有广泛的应用。离散傅里叶变换(Discrete Fourier Transform，DFT)是连续傅里叶变换在离散系统中的表示形式，由于DFT的计算量很大，因此在很长一段时间内其应用受到了很大的限制。20世纪60年代（1965年）由Cooley和Tukey提出了快速傅里叶变换(Fast Fourier Transform，FFT)算法，它是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上进行傅里叶变换一般需要做以下几步工作：

1. 在主机端，准备输入数据；
2. 在GPU设备端上分配内存，并将数据从主机复制到设备；（cudaMalloc,cudaMemcpy的接口 ）
3. 创建一个\(plan\), 调用函数\(cufftPlane1D/cufftPlane2D/cufftPlan3D\) 可以创建一个简单的Fourier变换。调用函数\(cufftPlanMany\) 则可以创建支持更多配置操作的变换计划。
   • \(cufftPlan1d()\): 针对单个1维信号
   • \(cufftPlan2d()\):针对单个2维信号
   • \(cufftPlan3d()\):针对单个3维信号
4. 执行\(plane\)。这一步可以使用\(cufftExecC2C()\)、\(cufftExecR2C()\)或\(cufftExecC2R()\)等函数完成上一步完成\(plane\)的计算任务。
5. 执行完成以下若不再需要该\(plan\)，则调用\(cufftDestroy()\)函数销毁该\(plan\) 及为其分配的计算资源。

3. cuFFT的傅里叶变换API接口类型（Fourier Transform Types）

​ \(cuFFT\) 库实现了三种不同类型的Fourier变换接口分为：\(C2C\)(复数变换到复数)，\(C2R\)(复数到实数), \(R2C\) (实数到复数)。本质上，这三种转换都可以被看做是复数域到复数域的变换，之所以这样划分，其最主要的考量是性能因素。例如，在一般的数字信号处理中，输入数据是一些离散的实数域上的采样点，这时候对它们做Fourier变换实际上就是\(R2C\)，根据埃尔米特对称性（Hermitian symmetry)，变换\(X_k=X_{N-k}^{*}\), \(*\) 代表共轭复数。\(cuFFT\) 的傅里叶变换则利用了这些冗余，将计算量降到最低。

​ 变换执行函数的单精度和双精度版本分别定义如下：

API	描述
\(cufftExecC2C()/cufftExecZ2Z()\)	单精度/双精度浮点数复数域到复数域的傅里叶变换
\(cufftExecR2C()/cufftExecD2Z()\)	单精度/双精度浮点数实数域到复数域的傅里叶变换（正变换）
\(cufftExecC2R()/cufftExecZ2D()\)	单精度/双精度浮点数复数域到实数域的傅里叶变换（逆变换）

4. 数据布局(Data Layout)

​ \(CUFFT\)库保含若干种数据类型，对于复数有\(cufftComplex/cufftDoubleComplex\) 两种数据类型，对于实数则分别有\(cufftReal/cufftDouble\) 两种数据类型 。

​ 根据转换结果的存储位置不同，\(FFT\)变换可分为就地变换(\(in-place\))和外部变换（\(out-place\))，前者直接在输入数据进行变换，而后者则会将变换后结果存入新的存储器地址。

​ 就地转换(\(in-place\)) 支持数据的两种布局：\(native\) 和 \(padded\)，前者用于获取最佳性能，而后者则用于与FFTW兼容。

​ 在\(padded\)布局中输出信号的开始地址与输入信号一样，换句话说，实数域到复数域变换的输入数据和复数域到实数域的输出数据必须被填充。在native布局中则没有填入要求。

​ 输入数据和输出数据的尺寸总结如下：

FFT type	input data size	output data size
\(C2C\)	\(X \space cufftComplex\)	\(X \space cufftComplex\)
\(C2R\)	\([\frac{X}{2}]+1 \space cufftComplex\)	\(X\) \(cufftReal\)
\(R2C^{*}\)	\(X\) \(cufftReal\)	\([\frac{X}{2}]+1 \space cufftComplex\)

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()来同时完成多个信号的fft；

cufftExecC2C():

• 第一个参数就是配置好的 cuFFT 句柄；
• 第二个参数为输入信号的首地址；
• 第三个参数为输出信号的首地址；
• 第四个参数为CUFFT_FORWARD表示执行的是\(fft\)正变换；CUFFT_INVERSE表示执行\(fft\)逆变换

需要注意的是，执行完\(fft\)之后，要对信号中的每个值乘以\(1/N\);

输出结果：