1

cuda バージョン 7.5cufftを使用して、FFT と逆 FFT を実行しています。関数を使用して逆 FFT を実行すると問題が発生しcufftExecC2R(.,.)ます。

実際、 で a を使用するbatch_size = 1と、cufftPlan1d(,)正しい結果が得られます。ただし、バッチ サイズを大きくすると、結果が正しくありません。

これを説明するために、サンプルの最小限のコードを貼り付けています。ざっと作成しただけなので、汚いコードは無視してください。

  #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

画像1

使用時の出力例BATCH_SIZE = 2 画像2

4

1 に答える 1

5

あなたが見逃している情報は、C2C 変換と C2R (または R2C) で期待される入力データのデータ形式の違いがあることを理解していないことです。

このセクションとCUFFT ドキュメントのこのセクションを読むことから始めてください。

それが言うことに注意してください:

これらの関数はそれぞれ、異なる入力データ レイアウトを必要とします

ただし、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
$
于 2016-10-22T23:24:36.513 に答える