使用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;? //創(chuàng)建句柄

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為要執(zhí)行 fft 的信號輸入類型及輸出類型都為復數；

//第四個參數BATCH表示要執(zhí)行 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表示執(zhí)行的是 fft 正變換；CUFFT_INVERSE表示執(zhí)行 fft 逆變換。

//!!!需要注意的是，執(zhí)行完逆 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