scipy.fftpack.fft2 和 cufft 之间的结果差异

Question

现在，我正在将我的 python 脚本移植到 CUDA 程序。 在我的 python 脚本中，使用了 scipy.fftpack.fft2。 为了验证 cufft 的结果，我使用 cufft 编写了示例程序。 但是，scipy.fftpack.fft2 和 cufft 之间似乎存在差异。

有什么建议吗？

python 脚本：

def test2():
   g = [18,19,19,23,24,24,23,24,24]
   g = numpy.array(g)
   g.shape = [3,3]
   G = fft2(g)

   print "---------------"
   print g
   print G
   return 

python脚本的结果：

   ---------------
    [[18 19 19]
     [23 24 24]
     [23 24 24]]
    [[ 198.+0.j   -3.+0.j   -3.+0.j]
     [ -15.+0.j    0.+0.j    0.+0.j]
     [ -15.+0.j    0.+0.j    0.+0.j]]

cuda 程序：

        cufftHandle plan;
        int nRows = 3;
        int nCols = 3;
        cufftPlan2d(&plan, nRows, nCols, CUFFT_R2C);
        float h_in[9] = {18,19,19,23,24,24,23,24,24};
        float* d_in;
        cudaMalloc(&d_in, sizeof(cufftComplex)*9); 
        cufftComplex* d_freq;
        cudaMalloc(&d_freq, sizeof(cufftComplex)*9); 
        cudaMemcpy(d_in,h_in,sizeof( cufftComplex)*9,cudaMemcpyHostToDevice);
        cufftExecR2C(inverse_plan, d_in, d_freq);
        cufftComplex* h_freq = (float2*)malloc(sizeof( cufftComplex)*9);    
        cudaMemcpy(h_freq,d_freq,sizeof( cufftComplex)*9,cudaMemcpyDeviceToHost);
        for(int i=0; i<9; i++) {
        printf("%i %f %f\n", i, h_freq[i].x, h_freq[i].y);
        }

cuda 程序的结果：

0 198.000000 -0.000001
1 -2.999996 -0.000001
2 -15.000000 0.000000
3 -0.000000 0.000000
4 -15.000000 0.000000
5 -0.000000 0.000000
6 497922732955248410000000000000.000000 8589934592.000000
7 572199135312371230000000000000.000000 8589934592.000000
8 -0.000000 0.000000

Answer 1

我不是袖带专家，但命名方式会泄露发生了什么：

• 在 numpy 中，您正在运行完整的 2D FFT。因为你的输入是真实的，所以输出是对称的，正如你所看到的：每一行（或每列）中的最后一项都等于前一项。

• 您可以利用这一点更快地运行 FFT，而在 numpy 中，这是通过 rfft2 函数实现的：

  >>> np.fft.rfft2(g)
  array([[ 198.+0.j,   -3.+0.j],
         [ -15.+0.j,    0.+0.j],
         [ -15.+0.j,    0.+0.j]])
  

• 我的猜测是您的CUFFT_R2C 计划名称中的R2C 的意思是“真实到复杂”，因此您要求的是np.rfft2 的等价物。如果您将数组中未使用的最后 3 项放在一边，结果几乎相同，除了舍入错误，以及您的 CUDA 实现使用 32 位浮点数，而不是 numpy 默认使用的 64 位浮点数.

• 快速的谷歌搜索显示CUFFT_C2C 和cufftExecR2C 是有效的袖带标识符。使用这些应该会产生您所追求的正确结果。为了更接近重现，重构您的代码并使用Z2Z 版本，它适用于double，而不是float。