FFTW3学习笔记3：FFTW 和 CUFFT 的使用对比

一、流程

1.使用cufftHandle创建句柄

2.使用cufftPlan1d(),cufftPlan3d(),cufftPlan3d(),cufftPlanMany()对句柄进行配置，主要是配置句柄对应的信号长度，信号类型，在内存中的存储形式等信息。 

cufftPlan1d()：针对单个 1 维信号

cufftPlan2d()：针对单个 2 维信号

cufftPlan3d()：针对单个 3 维信号

cufftPlanMany()：针对多个信号同时进行 fft

3.使用cufftExec()函数执行 fft

4.使用cufftDestroy()函数释放 GPU 资源

 

二、单个 1 维信号的 fft

假设要执行 fft 的信号data_dev的长度为N，并且已经传输到 GPU 显存中，data_dev数据的类型为cufftComplex，可以用一下方式产生主机段的data_dev。

cufftComplex *data_Host = (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;
}

然后用cudaMemcpy()将主机端的data_host拷贝到设备端的data_dev，即可用下述方法执行 fft ：

cufftHandle plan; // 创建cuFFT句柄
cufftPlan1d(&plan, N, CUFFT_C2C, BATCH);
cufftExecC2C(plan, data_dev, data_dev, CUFFT_FORWARD); // 执行 cuFFT，正变换

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

三、代码实现

#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;
}

四、用fftw和cufft实现傅里叶变换

1.创建C++的文件命名为fftw.cpp，配置fftw环境（环境配置移步：这里），复制以下代码

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fftw3.h"
#include <windows.h>
#include <Eigen/Dense>
#include <iostream>  
#include <opencv2/core/eigen.hpp>
#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;
using namespace Eigen;

#define COLS 3
#define ROWS 3

#pragma comment(lib, "libfftw3-3.lib") // double版本
//#pragma comment(lib, "libfftw3f-3.lib")// float版本
// #pragma comment(lib, "libfftw3l-3.lib")// long double版本

extern "C"    void iteration_mat1();

/**********************************主函数****************************************/
int main()
{

    fftw_complex*result_temp_din, *result_temp_out;
    fftw_plan p;

    result_temp_din = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*COLS*ROWS);
    result_temp_out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*COLS*ROWS);
    cout << "fftw" << endl;

    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {

            result_temp_din[i + j*COLS][0] = (i+1)*(j+1);
            cout << result_temp_din[i + j*COLS][0] << " ";
            result_temp_din[i + j*COLS][1] = 0;
        }
    }

    //forward fft  
    p = fftw_plan_dft_2d(ROWS, COLS, result_temp_din, result_temp_out, FFTW_FORWARD, FFTW_ESTIMATE);
    
    fftw_execute(p);
    cout << endl;
    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {
            cout << result_temp_out[i + j*COLS][0] << " ";//实部  
            cout << result_temp_out[i + j*COLS][1] << endl;//虚部  
        }
    }

    cout << "cuda" << endl;
    iteration_mat1();
    system("pause");
    return 0;
}

2.创建cuda文件命名为cufft.cu，配置环境（环境配置移步：这里），复制以下代码

注： cufftPlan2d(&p, ROWS, COLS, CUFFT_C2C); 看清楚rows和cols，千万别出错！

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <cufft.h>
#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

#define COLS 3
#define ROWS 3

extern "C"    void iteration_mat1()
{
    cufftComplex *result_temp_din = (cufftComplex*)malloc(COLS*ROWS * sizeof(cufftComplex));
    cufftHandle p;
    //输入赋值数据
    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {
            result_temp_din[i + j*COLS].x = (i + 1)*(j + 1);
            cout << result_temp_din[i + j*COLS].x << " ";
            result_temp_din[i + j*COLS].y = 0;
        }
    }
    cout << endl;

    size_t pitch;

    cufftComplex *t_result_temp_din;
    cudaMallocPitch((void**)&t_result_temp_din, &pitch, COLS * sizeof(cufftComplex), ROWS);
     
    cufftComplex *t_result_temp_out;
    cudaMallocPitch((void**)&t_result_temp_out, &pitch, COLS * sizeof(cufftComplex), ROWS);

    //将值辅到Device
    //cudaMemcpy2D(t_result_temp_din, pitch, result_temp_din, COLS * sizeof(cufftComplex), COLS * sizeof(cufftComplex), ROWS, cudaMemcpyHostToDevice);
    cudaMemcpy(t_result_temp_din,result_temp_din,  ROWS * sizeof(cufftComplex)* COLS, cudaMemcpyHostToDevice);

    //forward fft  制定变换规则
    cufftPlan2d(&p, ROWS, COLS, CUFFT_C2C);

    //执行变换
    cufftExecC2C(p, (cufftComplex*)t_result_temp_din, (cufftComplex*)t_result_temp_out, CUFFT_FORWARD);

    //将值辅到host
    cudaMemcpy(result_temp_din,  t_result_temp_out, ROWS * sizeof(cufftComplex)* COLS, cudaMemcpyDeviceToHost);
    //cudaMemcpy2D(result_temp_din, pitch, t_result_temp_out, COLS * sizeof(cufftComplex), sizeof(cufftComplex)* ROWS, COLS, cudaMemcpyDeviceToHost);


    //提取实部和虚部
    for (size_t j = 0; j < ROWS; j++)
    {
        for (size_t i = 0; i < COLS; i++)
        {
            cout << result_temp_din[i + j*COLS].x << " ";//实部  
            cout << result_temp_din[i + j*COLS].y << endl;//虚部  
        }
    }

}

3.执行结果：