使用CUDA计算傅里叶变换
使用cuda中自带的cufft库计算fft
再cuda设备中无法计算cucomplex与double、float的乘法，无法计算cufftcomplex与double、float的乘法，所以在实际计算中先试用double之间的乘法计算。然后再将结果使用cufftcomplex的构造函数来构造复数变量进行计算。
我目前没找到cuda中复数与double的直接乘法计算，如果有大佬看到可以教我一下
代码原文：
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "cufft.h"
#include <string>
#include <stdlib.h>
#include <ctime>
#include <Windows.h>
#include <mmSystem.h>
using namespace std;
#include<iostream>
#include <math.h>
//这是一个用于检查cuda函数是否正常运行的函数，直接复制粘贴就好
#define CHECK(call)\
{\
if ((call) != cudaSuccess)\
{\
printf("Error: %s:%d, ", __FILE__, __LINE__);\
printf("code:%d, reason: %s\n", (call), cudaGetErrorString(cudaGetLastError()));\
exit(1);\
}\
}
__global__ void com(cufftComplex* in, cufftComplex* out, int m) //cuda计算
{
int ii = threadIdx.x;
//int jj = threadIdx.y;
double a1 = 2.0;
double a2 = 3.0;
double a3 = a1 * a2;
cufftComplex aa = { a3, a3 };
cufftComplex bb = { 1, 2 };
out[ii] = cuCmulf(aa, bb);
}
const double PI = 3.141592653;
int main() {
const int NX = 16; //进行傅里叶变换的数据长度
const int fs = 4000000; //频率
double T = 2e-6; //周期
const int BATCH = 1;
int bei = fs / NX;
double timebei = T / NX;
double aaa = 1.0 / NX;
double Lk = 1 / T;
cufftHandle plan; //创建句柄
cufftComplex* data;
cufftComplex* data_cpu;
double t[NX]; //时间单位
data_cpu = (cufftComplex*)malloc(sizeof(cufftComplex) * NX);
if (data_cpu == NULL) return -1;
CHECK(cudaMalloc((void**)&data, sizeof(cufftComplex) * NX));
CHECK(cufftPlan1d(&plan, NX, CUFFT_C2C, 1)); //对句柄进行配置，主要是配置句柄对应的信号长度，信号类型，在内存中的存储形式等信息。
//第一个参数就是要配置的 cuFFT 句柄；
//第二个参数为要进行 fft 的信号的长度；
//第三个CUFFT_C2C为要执行 fft 的信号输入类型及输出类型都为复数；
//第四个参数BATCH表示要执行 fft 的信号的个数，新版的已经使用cufftPlanMany()来同时完成多个信号的 fft
//输入数据
//for (int i = 0; i < NX; ++i) //给一个三角函数时间信号，用于傅里叶变换
//{
// t[i] = i * timebei;
// //printf("%.15lf\n", t[i]);
// data_cpu[i].x = cos(PI * 8 * 0.25 * fs / (3 * T * T) * pow(t[i] - T / 2.0, 3) + 2 * PI * 0.25 * fs * t[i]);
// data_cpu[i].y = sin(PI * 8 * 0.25 * fs / (3 * T * T) * pow(t[i] - T / 2.0, 3) + 2 * PI * 0.25 * fs * t[i]);
//}
for (int i = 0; i < NX; ++i) //给一个好设定的时域信号，用于检验傅里叶变换
{
data_cpu[i].x =(float) i;
data_cpu[i].y = 0;
printf("(%.15lf,%.15lf)\n", data_cpu[i].x, data_cpu[i].y);
}
//for (int i = 0; i < NX; ++i) //输出时间信号的时间序列和模值
//{
// //printf("%f %f\n", data_cpu[i].x, data_cpu[i].y);
// double mo = sqrt(data_cpu[i].x * data_cpu[i].x + data_cpu[i].y * data_cpu[i].y);
// double time = i*timebei;
// printf("%.15f %.20f\n", time, mo);
//}
//数据传输cpu->gpu
CHECK(cudaMemcpy(data, data_cpu, sizeof(cufftComplex) * NX * BATCH, cudaMemcpyHostToDevice));
CHECK(cudaDeviceSynchronize());
clock_t ct1, ct2,ct3,ct4; //添加计时
ct1 = clock();
CHECK(cufftExecC2C(plan, data, data, CUFFT_FORWARD));
//CHECK(cufftExecC2C(plan, data, data, CUFFT_INVERSE) != CUFFT_SUCCESS);
//第一个参数就是配置好的 cuFFT 句柄；
//第二个参数为输入信号的首地址；
//第三个参数为输出信号的首地址；
//第四个参数CUFFT_FORWARD表示执行的是 fft 正变换；CUFFT_INVERSE表示执行 fft 逆变换。
//!!!需要注意的是，执行完逆 fft 之后，要对信号中的每个值乘以 1/N
CHECK(cudaDeviceSynchronize());
ct2 = clock();
CHECK(cudaMemcpy(data_cpu, data, sizeof(cufftComplex) * NX, cudaMemcpyDeviceToHost));
CHECK(cudaDeviceSynchronize());
//for (int i = 0; i < NX; ++i) //这个不是输出傅里叶变换，这个是输出频域信号的
//{
// //printf("%f %f\n", data_cpu[i].x, data_cpu[i].y);
// double mo = sqrt(data_cpu[i].x * data_cpu[i].x + data_cpu[i].y * data_cpu[i].y);
// //double Fre = i*F;
// printf("%d %.15f\n", i * bei, mo * aaa);
//}
for (int i = 0; i < NX; ++i) //这个输出的是傅里叶变换的实部和虚部
{
printf("(%.15lf,%.15lf)\n", data_cpu[i].x, data_cpu[i].y);
}
ct3 = clock();
com << <1, 16 >> > (data, data, NX); //测试计算的函数
CHECK(cudaDeviceSynchronize());
ct4 = clock();
//数据传输gpu->cpu
CHECK(cudaMemcpy(data_cpu, data, sizeof(cufftComplex) * NX, cudaMemcpyDeviceToHost));
CHECK(cudaDeviceSynchronize());
printf("\n");
for (int i = 0; i < NX; ++i)
{
printf("(%.15lf,%.15lf)\n", data_cpu[i].x,data_cpu[i].y);
}
double cputime = ((double)(ct2 - ct1)) / CLOCKS_PER_SEC;
double cputime2 = ((double)(ct4 - ct3)) / CLOCKS_PER_SEC;
cout << "ct1=" << ct1 << endl;
cout << "ct2=" << ct2 << endl;
cout << "ct3=" << ct3 << endl;
cout << "ct4=" << ct4 << endl;
cout << "傅里叶变换时间=" << cputime << endl;
cout << "计算操作的时间=" << cputime2 << endl;
//printf("变换用时: %dms\n", End - Start);
cufftDestroy(plan); //释放 GPU 资源
cudaFree(data);
//printf("CUFFT_FORWARD:\n");
system("pause");
return 0;
}
运行结果，16个一个单位，前面一部分忘记分开了
代码参考来源
https://blog.csdn.net/lqbird/article/details/127054520?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522167591186216800217024470%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=167591186216800217024470&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduend~default-1-127054520-null-null.142^v73^insert_down4,201^v4^add_ask,239^v1^insert_chatgpt&utm_term=CUFFT&spm=1018.2226.3001.4187
