cuFFT неверные результаты только при запуске из комплекса

Question

Мне помогли раньше в this Ответьте, чтобы реализовать трансформацию на месте, и она работает хорошо, но ТОЛЬКО, если я начну с реальных данных. Если я начинаю с сложных данных, результаты после IFT + FFT ошибочны, и это происходит только в версии на месте, у меня есть отличные результаты с внеплотной версией этого преобразования. Это код:

#include <stdio.h> 
#include <stdlib.h> 
#include <cuda_runtime.h> 
#include <complex.h> 
#include <cuComplex.h> 
#include <cufft.h> 
#include <cufftXt.h> 

#define N 4 
#define N_PAD (2*(N/2+1)) 

void print_3D_Real(double *array){ 
    printf("\nPrinting 3D real matrix \n"); 
    unsigned long int idx; 
    for (int z = 0; z < N; z++){ 
     printf("---------------------------------------------------------------------------- plane %d below\n", z); 
     for (int x = 0; x < N; x++){ 
      for (int y = 0; y < N; y++){ 
       idx = z + N_PAD * (y + x * N); 
       printf("%.3f \t", array[idx]); 
      } 
      printf("\n"); 
     } 
    } 
} 

void print_3D_Comp(cuDoubleComplex *array){ 
    printf("\nPrinting 3D complex matrix \n"); 
    unsigned long int idx; 
    for (int z = 0; z < (N/2+1); z++){ 
     printf("---------------------------------------------------------------------------- plane %d below\n", z); 
     for (int x = 0; x < N; x++){ 
      for (int y = 0; y < N; y++){ 
       idx = z + (N/2+1) * (y + x * N); 
       printf("%+.3f%+.3fi \t", array[idx].x, array[idx].y); 
      } 
      printf("\n"); 
     } 
    } 
} 

// Main function 
int main(int argc, char **argv){ 
    CU_ERR_CHECK(cudaSetDevice(0)); 
    unsigned long int idx, in_mem_size, out_mem_size; 
    cuDoubleComplex *in = NULL, *d_in = NULL; 
    double *out = NULL, *d_out = NULL; 
    cufftHandle plan_r2c, plan_c2r; 
    in_mem_size = sizeof(cuDoubleComplex) * N*N*(N/2+1); 
    out_mem_size = in_mem_size; 

    in = (cuDoubleComplex *) malloc (in_mem_size); 
    out = (double *) in; 

    cudaMalloc((void **)&d_in, in_mem_size); 
    d_out = (double *) d_in; 

    cufftPlan3d(&plan_c2r, N, N, N, CUFFT_Z2D); 
    cufftPlan3d(&plan_r2c, N, N, N, CUFFT_D2Z); 

    memset(in, 0, in_mem_size); 
    memset(out, 0, out_mem_size); 

    // Initial complex data 
    for (int x = 0; x < N; x++){ 
     for (int y = 0; y < N; y++){ 
      for (int z = 0; z < (N/2+1); z++){ 
       idx = z + (N/2+1) * (y + x * N); 
       in[idx].x = idx; 
      } 
     } 
    } 
    print_3D_Comp(in); 
    cudaMemcpy(d_in, in, in_mem_size, cudaMemcpyHostToDevice); 

    cufftExecZ2D(plan_c2r, (cufftDoubleComplex *)d_in, (cufftDoubleReal *)d_out); 

    cudaMemcpy(out, d_out, out_mem_size, cudaMemcpyDeviceToHost); 
    // Normalisation 
    for (int i = 0; i < N*N*N_PAD; i++) 
     out[i] /= (N*N*N); 
    print_3D_Real(out); 
    cudaMemcpy(d_out, out, out_mem_size, cudaMemcpyHostToDevice); 

    cufftExecD2Z(plan_r2c, (cufftDoubleReal *)d_out, (cufftDoubleComplex *)d_in); 

    cudaMemcpy(in, d_in, in_mem_size, cudaMemcpyDeviceToHost)); 

    print_3D_Comp(in); 

    cudaDeviceReset(); 
    return 0; 
} 

Выход моей программы на этом pastebin.

Может ли кто-нибудь направить меня по правильному пути? Заранее большое спасибо.

Answer 1

Прежде всего, ваш код не компилируется.

В своем наиболее общем определении преобразование Фурье выполняет преобразование из одного сложного домена в другой комплексный домен, и эта операция должна быть reversible.

Однако C2R и R2C являются особыми случаями с предположением, что сигнал полностью представим в одном из двух доменов («область времени») как чисто вещественный сигнал (все мнимые компоненты равны нулю).

Однако должно быть очевидно, что будут некоторые сложные «частотные» представления домена, которые не могут быть представлены сигналом чисто реального времени. Если счетчик был прав (любой сигнал сложной частотной области может быть представлен как сигнал чисто реального времени), то БПФ не может быть обратимым для сигнала сложной временной области (так как все наборы данных в частотной области сопоставляются с данными в реальном масштабе времени в реальном времени устанавливает.)

Поэтому вы не можете выбирать произвольные данные в частотной области и ожидать, что она правильно отобразится в сигнал чисто реального времени. (*)

В качестве демонстрации, измените входные данные установить на следующее:

  in[idx].x = (idx)?0:1; 

, и я верю, что вы получите «проходной» тестовый случай.

Кроме того, ваше утверждение о том, что «у меня есть отличные результаты с неправильной версией этого преобразования», я считаю, что не поддерживается, если вы действительно используете этот конкретный набор данных, как указано в вашем вопросе. Если вы не согласны, напишите полный код, демонстрирующий ваш тестовый пример с преобразованием вне места, которое в остальном идентично вашему опубликованному коду.

Наконец, мы можем проверить это с помощью fftw. Преобразование вашей программы использовать FFTW вместо CUFFT производит точно такой же вывод:

$ cat t355.cpp 
#include <stdio.h> 
#include <stdlib.h> 
#include <fftw3.h> 
#include <string.h> 

#define N 4 
#define N_PAD (2*(N/2+1)) 

void print_3D_Real(double *array){ 
    printf("\nPrinting 3D real matrix \n"); 
    unsigned long int idx; 
    for (int z = 0; z < N; z++){ 
     printf("---------------------------------------------------------------------------- plane %d below\n", z); 
     for (int x = 0; x < N; x++){ 
      for (int y = 0; y < N; y++){ 
       idx = z + N_PAD * (y + x * N); 
       printf("%.3f \t", array[idx]); 
      } 
      printf("\n"); 
     } 
    } 
} 

void print_3D_Comp(fftw_complex *array){ 
    printf("\nPrinting 3D complex matrix \n"); 
    unsigned long int idx; 
    for (int z = 0; z < (N/2+1); z++){ 
     printf("---------------------------------------------------------------------------- plane %d below\n", z); 
     for (int x = 0; x < N; x++){ 
      for (int y = 0; y < N; y++){ 
       idx = z + (N/2+1) * (y + x * N); 
       printf("%+.3f%+.3fi \t", array[idx][0], array[idx][1]); 
      } 
      printf("\n"); 
     } 
    } 
} 

// Main function 
int main(int argc, char **argv){ 
    unsigned long int idx, in_mem_size, out_mem_size; 
    fftw_complex *in = NULL; 
    double *out = NULL; 
    in_mem_size = sizeof(fftw_complex) * N*N*(N/2+1); 
    out_mem_size = in_mem_size; 

    in = (fftw_complex *) malloc (in_mem_size); 
    out = (double *) in; 

    memset(in, 0, in_mem_size); 
    memset(out, 0, out_mem_size); 

    // Initial complex data 
    for (int x = 0; x < N; x++){ 
     for (int y = 0; y < N; y++){ 
      for (int z = 0; z < (N/2+1); z++){ 
       idx = z + (N/2+1) * (y + x * N); 
       in[idx][0] = idx; 
      } 
     } 
    } 
    print_3D_Comp(in); 

    fftw_plan plan_c2r = fftw_plan_dft_c2r_3d(N, N, N, in, out, FFTW_ESTIMATE); 
    fftw_plan plan_r2c = fftw_plan_dft_r2c_3d(N, N, N, out, in, FFTW_ESTIMATE); 

    fftw_execute(plan_c2r); 
    // Normalisation 
    for (int i = 0; i < N*N*N_PAD; i++) 
     out[i] /= (N*N*N); 
    print_3D_Real(out); 
    fftw_execute(plan_r2c); 



    print_3D_Comp(in); 

    return 0; 
} 
$ g++ t355.cpp -o t355 -lfftw3 
$ ./t355 

Printing 3D complex matrix 
---------------------------------------------------------------------------- plane 0 below 
+0.000+0.000i +3.000+0.000i +6.000+0.000i +9.000+0.000i 
+12.000+0.000i +15.000+0.000i +18.000+0.000i +21.000+0.000i 
+24.000+0.000i +27.000+0.000i +30.000+0.000i +33.000+0.000i 
+36.000+0.000i +39.000+0.000i +42.000+0.000i +45.000+0.000i 
---------------------------------------------------------------------------- plane 1 below 
+1.000+0.000i +4.000+0.000i +7.000+0.000i +10.000+0.000i 
+13.000+0.000i +16.000+0.000i +19.000+0.000i +22.000+0.000i 
+25.000+0.000i +28.000+0.000i +31.000+0.000i +34.000+0.000i 
+37.000+0.000i +40.000+0.000i +43.000+0.000i +46.000+0.000i 
---------------------------------------------------------------------------- plane 2 below 
+2.000+0.000i +5.000+0.000i +8.000+0.000i +11.000+0.000i 
+14.000+0.000i +17.000+0.000i +20.000+0.000i +23.000+0.000i 
+26.000+0.000i +29.000+0.000i +32.000+0.000i +35.000+0.000i 
+38.000+0.000i +41.000+0.000i +44.000+0.000i +47.000+0.000i 

Printing 3D real matrix 
---------------------------------------------------------------------------- plane 0 below 
23.500 -1.500 -1.500 -1.500 
-6.000 0.000 0.000 0.000 
-6.000 0.000 0.000 0.000 
-6.000 0.000 0.000 0.000 
---------------------------------------------------------------------------- plane 1 below 
-0.500 0.750 0.000 -0.750 
3.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 
-3.000 0.000 0.000 0.000 
---------------------------------------------------------------------------- plane 2 below 
0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 
---------------------------------------------------------------------------- plane 3 below 
-0.500 -0.750 0.000 0.750 
-3.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 
3.000 0.000 0.000 0.000 

Printing 3D complex matrix 
---------------------------------------------------------------------------- plane 0 below 
+0.000+0.000i +6.000+0.000i +6.000+0.000i +6.000+0.000i 
+24.000+0.000i +30.000+0.000i +30.000+0.000i +30.000+0.000i 
+24.000+0.000i +30.000+0.000i +30.000+0.000i +30.000+0.000i 
+24.000+0.000i +30.000+0.000i +30.000+0.000i +30.000+0.000i 
---------------------------------------------------------------------------- plane 1 below 
+1.000+0.000i +4.000+0.000i +7.000+0.000i +10.000+0.000i 
+13.000+0.000i +16.000+0.000i +19.000+0.000i +22.000+0.000i 
+25.000+0.000i +28.000+0.000i +31.000+0.000i +34.000+0.000i 
+37.000+0.000i +40.000+0.000i +43.000+0.000i +46.000+0.000i 
---------------------------------------------------------------------------- plane 2 below 
+2.000+0.000i +8.000+0.000i +8.000+0.000i +8.000+0.000i 
+26.000+0.000i +32.000+0.000i +32.000+0.000i +32.000+0.000i 
+26.000+0.000i +32.000+0.000i +32.000+0.000i +32.000+0.000i 
+26.000+0.000i +32.000+0.000i +32.000+0.000i +32.000+0.000i 
$ 

(*) Вы можете возразить, если вы хотите, то, что комплексно-сопряженная функция симметрия C2R и r2c преобразований должны учитывают правильное отображение всех возможных сложных «частотных» сигналов домена в уникальные, чисто реальные «временные» сигналы домена. Я не утверждаю, что это не так, с двумя точками данных:

1. Пример кода в этом вопросе.
2. Поскольку комплексное пространство в преобразовании C2R или R2C численно больше реального пространства (в (2*(N/2+1))/N), то разумно, что не может быть уникального отображения всех возможных комплексных сигналов 1: 1 в уникальный реальный сигналы. И единственное отображение 1: 1 необходимо для полной обратимости.

Для дополнительного фона на возможность отсутствия симметрии в случайных данных, обратите внимание на дискуссии вокруг CUFFT_COMPATIBILITY_FFTW_ASYMMETRIC в cufft documentation.