使用 cuda cuFFT 从复数转换为实数时输出不正确
Incorrect output when transforming from complex to real number using cuda cuFFT
我正在使用 cuda 版本 7.5 cufft 来执行一些 FFT 和逆 FFT。
使用 cufftExecC2R(.,.) 函数执行逆 FFT 时遇到问题。
实际上，当我在 cufftPlan1d(,) 中使用 batch_size = 1 时，我得到了正确的结果。但是，当我增加批量大小时，结果不正确。
我正在粘贴一个示例最小代码来说明这一点。请忽略我刚刚快速创建的代码的脏代码。
#include <cufft.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctime>
#include <iostream>
typedef float2 Complex;
void iTest(int argc, char** argv);
#define SIGNAL_SIZE 9
#define BATCH_SIZE 2
int main(int argc, char** argv) {
iTest(argc, argv);
return 0;
}
void iProcess(Complex *x, double *y, size_t n) {
cufftComplex *deviceData;
cudaMalloc(reinterpret_cast<void**>(&deviceData),
SIGNAL_SIZE * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(deviceData, x, SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyHostToDevice);
cufftResult cufftStatus;
cufftHandle handle;
cufftStatus = cufftPlan1d(&handle, SIGNAL_SIZE, CUFFT_C2C, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
cufftComplex *d_complex;
cudaMalloc(reinterpret_cast<void**>(&d_complex),
sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2C(handle, deviceData, d_complex, CUFFT_FORWARD);
if (cufftStatus != cudaSuccess) {
printf("cufftExecR2C failed!");
}
cufftComplex *hostOutputData = (cufftComplex*)malloc(
(SIGNAL_SIZE) * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(hostOutputData, d_complex,
SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting COMPLEX" << "\n";
for (int j = 0; j < (SIGNAL_SIZE) * BATCH_SIZE; j++)
printf("%i \t %f \t %f\n", j, hostOutputData[j].x, hostOutputData[j].y);
//! convert complex to real
cufftHandle c2r_handle;
cufftStatus = cufftPlan1d(&c2r_handle, SIGNAL_SIZE, CUFFT_C2R, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
cufftReal *d_odata;
cudaMalloc(reinterpret_cast<void**>(&d_odata),
sizeof(cufftReal) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2R(c2r_handle, d_complex, d_odata);
cufftReal odata[SIGNAL_SIZE * BATCH_SIZE];
cudaMemcpy(odata, d_odata, sizeof(cufftReal) * SIGNAL_SIZE * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting REAL" << "\n";
for (int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; i++) {
std::cout << i << " \t" << odata[i]/(SIGNAL_SIZE) << "\n";
}
cufftDestroy(handle);
cudaFree(deviceData);
}
void iTest(int argc, char** argv) {
Complex* h_signal = reinterpret_cast<Complex*>(
malloc(sizeof(Complex) * SIGNAL_SIZE * BATCH_SIZE));
std::cout << "\nPrinting INPUT" << "\n";
for (unsigned int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; ++i) {
h_signal[i].x = rand() / static_cast<float>(RAND_MAX);
h_signal[i].y = 0;
std::cout << i << "\t" << h_signal[i].x << "\n";
}
std::cout << "\n";
double y[SIGNAL_SIZE * BATCH_SIZE];
iProcess(h_signal, y, 1);
}
我无法找出代码中的错误在哪里以及缺少哪些信息。
使用 BATCH_SIZE = 1
时的示例输出
使用 BATCH_SIZE = 2 时的示例输出
您缺少的信息是您不了解 C2C 转换与 C2R（或 R2C）预期的输入数据存在数据格式差异。
您应该先阅读 this section and this section CUFFT 文档。
注意上面写着：
Each of those functions demands different input data layout
但是您将对 C2C 转换正确的输入数据直接传递给 C2R 转换。那不行。
IMO 最直接的解决方案是将您的所有工作转换为 C2C 转换类型。 C2C 变换可以支持正向（例如 "real-to-complex"）和反向（例如 "complex-to-real"）。您使用的 C2R 转换类型也可以支持 "complex-to-real"，但是您将用于 C2R 的数据排列与您将用于具有指定反向路径的 C2C 的数据排列不同 ，否则 相同的变换 。你没有考虑到这一点。
这里是一个工作示例，显示了代码的修改版本，它对正向和反向路径使用 C2C，并正确地再现了批量大小为 2 的输入：
$ cat t19.cu
#include <cufft.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctime>
#include <iostream>
typedef float2 Complex;
void iTest(int argc, char** argv);
#define SIGNAL_SIZE 9
#define BATCH_SIZE 2
int main(int argc, char** argv) {
iTest(argc, argv);
return 0;
}
void iProcess(Complex *x, double *y, size_t n) {
cufftComplex *deviceData;
cudaMalloc(reinterpret_cast<void**>(&deviceData),
SIGNAL_SIZE * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(deviceData, x, SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyHostToDevice);
cufftResult cufftStatus;
cufftHandle handle;
cufftStatus = cufftPlan1d(&handle, SIGNAL_SIZE, CUFFT_C2C, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
cufftComplex *d_complex;
cudaMalloc(reinterpret_cast<void**>(&d_complex),
sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2C(handle, deviceData, d_complex, CUFFT_FORWARD);
if (cufftStatus != cudaSuccess) {
printf("cufftExecR2C failed!");
}
cufftComplex *hostOutputData = (cufftComplex*)malloc(
(SIGNAL_SIZE) * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(hostOutputData, d_complex,
SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting COMPLEX" << "\n";
for (int j = 0; j < (SIGNAL_SIZE) * BATCH_SIZE; j++)
printf("%i \t %f \t %f\n", j, hostOutputData[j].x, hostOutputData[j].y);
//! convert complex to real
/* cufftHandle c2r_handle;
cufftStatus = cufftPlan1d(&c2r_handle, SIGNAL_SIZE, CUFFT_C2R, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
*/
cufftComplex *d_odata;
cudaMalloc(reinterpret_cast<void**>(&d_odata),
sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2C(handle, d_complex, d_odata, CUFFT_INVERSE);
cufftComplex odata[SIGNAL_SIZE * BATCH_SIZE];
cudaMemcpy(odata, d_odata, sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting REAL" << "\n";
for (int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; i++) {
std::cout << i << " \t" << odata[i].x/(SIGNAL_SIZE) << "\n";
}
cufftDestroy(handle);
cudaFree(deviceData);
}
void iTest(int argc, char** argv) {
Complex* h_signal = reinterpret_cast<Complex*>(
malloc(sizeof(Complex) * SIGNAL_SIZE * BATCH_SIZE));
std::cout << "\nPrinting INPUT" << "\n";
for (unsigned int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; ++i) {
h_signal[i].x = rand() / static_cast<float>(RAND_MAX);
h_signal[i].y = 0;
std::cout << i << "\t" << h_signal[i].x << "\n";
}
std::cout << "\n";
double y[SIGNAL_SIZE * BATCH_SIZE];
iProcess(h_signal, y, 1);
}
$ nvcc -arch=sm_61 -o t19 t19.cu -lcufft
t19.cu: In function ‘void iProcess(Complex*, double*, size_t)’:
t19.cu:34:32: warning: comparison between ‘cufftResult {aka enum cufftResult_t}’ and ‘enum cudaError’ [-Wenum-compare]
if (cufftStatus != cudaSuccess) {
^
t19.cu:43:32: warning: comparison between ‘cufftResult {aka enum cufftResult_t}’ and ‘enum cudaError’ [-Wenum-compare]
if (cufftStatus != cudaSuccess) {
^
$ cuda-memcheck ./t19
========= CUDA-MEMCHECK
Printing INPUT
0 0.840188
1 0.394383
2 0.783099
3 0.79844
4 0.911647
5 0.197551
6 0.335223
7 0.76823
8 0.277775
9 0.55397
10 0.477397
11 0.628871
12 0.364784
13 0.513401
14 0.95223
15 0.916195
16 0.635712
17 0.717297
Printing COMPLEX
0 5.306536 0.000000
1 0.015338 -0.734991
2 -0.218001 0.740248
3 0.307508 -0.706533
4 1.022732 0.271765
5 1.022732 -0.271765
6 0.307508 0.706533
7 -0.218001 -0.740248
8 0.015338 0.734991
9 5.759857 0.000000
10 -0.328981 0.788566
11 0.055356 -0.521014
12 -0.127504 0.581872
13 0.014066 0.123027
14 0.014066 -0.123027
15 -0.127504 -0.581872
16 0.055356 0.521014
17 -0.328981 -0.788566
Printing REAL
0 0.840188
1 0.394383
2 0.783099
3 0.79844
4 0.911647
5 0.197551
6 0.335223
7 0.76823
8 0.277775
9 0.55397
10 0.477397
11 0.628871
12 0.364784
13 0.513401
14 0.95223
15 0.916195
16 0.635712
17 0.717297
========= ERROR SUMMARY: 0 errors
$
我正在使用 cuda 版本 7.5 cufft 来执行一些 FFT 和逆 FFT。
使用 cufftExecC2R(.,.) 函数执行逆 FFT 时遇到问题。
实际上，当我在 cufftPlan1d(,) 中使用 batch_size = 1 时，我得到了正确的结果。但是，当我增加批量大小时，结果不正确。
我正在粘贴一个示例最小代码来说明这一点。请忽略我刚刚快速创建的代码的脏代码。
#include <cufft.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctime>
#include <iostream>
typedef float2 Complex;
void iTest(int argc, char** argv);
#define SIGNAL_SIZE 9
#define BATCH_SIZE 2
int main(int argc, char** argv) {
iTest(argc, argv);
return 0;
}
void iProcess(Complex *x, double *y, size_t n) {
cufftComplex *deviceData;
cudaMalloc(reinterpret_cast<void**>(&deviceData),
SIGNAL_SIZE * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(deviceData, x, SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyHostToDevice);
cufftResult cufftStatus;
cufftHandle handle;
cufftStatus = cufftPlan1d(&handle, SIGNAL_SIZE, CUFFT_C2C, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
cufftComplex *d_complex;
cudaMalloc(reinterpret_cast<void**>(&d_complex),
sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2C(handle, deviceData, d_complex, CUFFT_FORWARD);
if (cufftStatus != cudaSuccess) {
printf("cufftExecR2C failed!");
}
cufftComplex *hostOutputData = (cufftComplex*)malloc(
(SIGNAL_SIZE) * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(hostOutputData, d_complex,
SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting COMPLEX" << "\n";
for (int j = 0; j < (SIGNAL_SIZE) * BATCH_SIZE; j++)
printf("%i \t %f \t %f\n", j, hostOutputData[j].x, hostOutputData[j].y);
//! convert complex to real
cufftHandle c2r_handle;
cufftStatus = cufftPlan1d(&c2r_handle, SIGNAL_SIZE, CUFFT_C2R, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
cufftReal *d_odata;
cudaMalloc(reinterpret_cast<void**>(&d_odata),
sizeof(cufftReal) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2R(c2r_handle, d_complex, d_odata);
cufftReal odata[SIGNAL_SIZE * BATCH_SIZE];
cudaMemcpy(odata, d_odata, sizeof(cufftReal) * SIGNAL_SIZE * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting REAL" << "\n";
for (int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; i++) {
std::cout << i << " \t" << odata[i]/(SIGNAL_SIZE) << "\n";
}
cufftDestroy(handle);
cudaFree(deviceData);
}
void iTest(int argc, char** argv) {
Complex* h_signal = reinterpret_cast<Complex*>(
malloc(sizeof(Complex) * SIGNAL_SIZE * BATCH_SIZE));
std::cout << "\nPrinting INPUT" << "\n";
for (unsigned int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; ++i) {
h_signal[i].x = rand() / static_cast<float>(RAND_MAX);
h_signal[i].y = 0;
std::cout << i << "\t" << h_signal[i].x << "\n";
}
std::cout << "\n";
double y[SIGNAL_SIZE * BATCH_SIZE];
iProcess(h_signal, y, 1);
}
我无法找出代码中的错误在哪里以及缺少哪些信息。
使用 BATCH_SIZE = 1
使用 BATCH_SIZE = 2 时的示例输出
您缺少的信息是您不了解 C2C 转换与 C2R（或 R2C）预期的输入数据存在数据格式差异。
您应该先阅读 this section and this section CUFFT 文档。
注意上面写着：
Each of those functions demands different input data layout
但是您将对 C2C 转换正确的输入数据直接传递给 C2R 转换。那不行。
IMO 最直接的解决方案是将您的所有工作转换为 C2C 转换类型。 C2C 变换可以支持正向（例如 "real-to-complex"）和反向（例如 "complex-to-real"）。您使用的 C2R 转换类型也可以支持 "complex-to-real"，但是您将用于 C2R 的数据排列与您将用于具有指定反向路径的 C2C 的数据排列不同 ，否则 相同的变换 。你没有考虑到这一点。
这里是一个工作示例，显示了代码的修改版本，它对正向和反向路径使用 C2C，并正确地再现了批量大小为 2 的输入：
$ cat t19.cu
#include <cufft.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctime>
#include <iostream>
typedef float2 Complex;
void iTest(int argc, char** argv);
#define SIGNAL_SIZE 9
#define BATCH_SIZE 2
int main(int argc, char** argv) {
iTest(argc, argv);
return 0;
}
void iProcess(Complex *x, double *y, size_t n) {
cufftComplex *deviceData;
cudaMalloc(reinterpret_cast<void**>(&deviceData),
SIGNAL_SIZE * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(deviceData, x, SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyHostToDevice);
cufftResult cufftStatus;
cufftHandle handle;
cufftStatus = cufftPlan1d(&handle, SIGNAL_SIZE, CUFFT_C2C, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
cufftComplex *d_complex;
cudaMalloc(reinterpret_cast<void**>(&d_complex),
sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2C(handle, deviceData, d_complex, CUFFT_FORWARD);
if (cufftStatus != cudaSuccess) {
printf("cufftExecR2C failed!");
}
cufftComplex *hostOutputData = (cufftComplex*)malloc(
(SIGNAL_SIZE) * BATCH_SIZE * sizeof(cufftComplex));
cudaMemcpy(hostOutputData, d_complex,
SIGNAL_SIZE * sizeof(cufftComplex) * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting COMPLEX" << "\n";
for (int j = 0; j < (SIGNAL_SIZE) * BATCH_SIZE; j++)
printf("%i \t %f \t %f\n", j, hostOutputData[j].x, hostOutputData[j].y);
//! convert complex to real
/* cufftHandle c2r_handle;
cufftStatus = cufftPlan1d(&c2r_handle, SIGNAL_SIZE, CUFFT_C2R, BATCH_SIZE);
if (cufftStatus != cudaSuccess) {
printf("cufftPlan1d failed!");
}
*/
cufftComplex *d_odata;
cudaMalloc(reinterpret_cast<void**>(&d_odata),
sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE);
cufftStatus = cufftExecC2C(handle, d_complex, d_odata, CUFFT_INVERSE);
cufftComplex odata[SIGNAL_SIZE * BATCH_SIZE];
cudaMemcpy(odata, d_odata, sizeof(cufftComplex) * SIGNAL_SIZE * BATCH_SIZE,
cudaMemcpyDeviceToHost);
std::cout << "\nPrinting REAL" << "\n";
for (int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; i++) {
std::cout << i << " \t" << odata[i].x/(SIGNAL_SIZE) << "\n";
}
cufftDestroy(handle);
cudaFree(deviceData);
}
void iTest(int argc, char** argv) {
Complex* h_signal = reinterpret_cast<Complex*>(
malloc(sizeof(Complex) * SIGNAL_SIZE * BATCH_SIZE));
std::cout << "\nPrinting INPUT" << "\n";
for (unsigned int i = 0; i < SIGNAL_SIZE * BATCH_SIZE; ++i) {
h_signal[i].x = rand() / static_cast<float>(RAND_MAX);
h_signal[i].y = 0;
std::cout << i << "\t" << h_signal[i].x << "\n";
}
std::cout << "\n";
double y[SIGNAL_SIZE * BATCH_SIZE];
iProcess(h_signal, y, 1);
}
$ nvcc -arch=sm_61 -o t19 t19.cu -lcufft
t19.cu: In function ‘void iProcess(Complex*, double*, size_t)’:
t19.cu:34:32: warning: comparison between ‘cufftResult {aka enum cufftResult_t}’ and ‘enum cudaError’ [-Wenum-compare]
if (cufftStatus != cudaSuccess) {
^
t19.cu:43:32: warning: comparison between ‘cufftResult {aka enum cufftResult_t}’ and ‘enum cudaError’ [-Wenum-compare]
if (cufftStatus != cudaSuccess) {
^
$ cuda-memcheck ./t19
========= CUDA-MEMCHECK
Printing INPUT
0 0.840188
1 0.394383
2 0.783099
3 0.79844
4 0.911647
5 0.197551
6 0.335223
7 0.76823
8 0.277775
9 0.55397
10 0.477397
11 0.628871
12 0.364784
13 0.513401
14 0.95223
15 0.916195
16 0.635712
17 0.717297
Printing COMPLEX
0 5.306536 0.000000
1 0.015338 -0.734991
2 -0.218001 0.740248
3 0.307508 -0.706533
4 1.022732 0.271765
5 1.022732 -0.271765
6 0.307508 0.706533
7 -0.218001 -0.740248
8 0.015338 0.734991
9 5.759857 0.000000
10 -0.328981 0.788566
11 0.055356 -0.521014
12 -0.127504 0.581872
13 0.014066 0.123027
14 0.014066 -0.123027
15 -0.127504 -0.581872
16 0.055356 0.521014
17 -0.328981 -0.788566
Printing REAL
0 0.840188
1 0.394383
2 0.783099
3 0.79844
4 0.911647
5 0.197551
6 0.335223
7 0.76823
8 0.277775
9 0.55397
10 0.477397
11 0.628871
12 0.364784
13 0.513401
14 0.95223
15 0.916195
16 0.635712
17 0.717297
========= ERROR SUMMARY: 0 errors
$