c++ - 使用 cuda cuFFT 从复数转换为实数时输出不正确

Question

我正在使用 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 时的示例输出

Answer 1

您缺少的信息是您不了解 C2C 转换与 C2R(或 R2C)预期的输入数据存在数据格式差异。

您应该先阅读 this section和 this section CUFFT 文档。

注意它说:

Each of those functions demands different input data layout

但是您将对 C2C 转换正确的输入数据直接传递给 C2R 转换。那是行不通的。

IMO 最直接的解决方案是将您的所有工作转换为 C2C 转换类型。 C2C 转换可以支持正向(例如“实到复杂”)和反向(例如“复杂到实”)。您使用的 C2R 转换类型也可以支持“从复杂到真实”，但是您将用于 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
$