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); // 初始为0cudaMemcpy(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 3extern "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);//将值辅到hostcudaMemcpy(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.执行结果：

转载于:https://www.cnblogs.com/aiguona/p/9485606.html