非方阵的错误 cuFFT 2D 卷积结果答案

【问题标题】：Wrong cuFFT 2D Convolution results with non square matrix非方阵的错误 cuFFT 2D 卷积结果
【发布时间】：2014-05-15 10:44:42
【问题描述】：

我正在尝试使用“FFT + point_wise_product + iFFT”方法执行 2D 卷积。使用 NxN 矩阵该方法运行良好，但是使用 非方阵 结果不正确。我已经阅读了整个 cuFFT 文档，寻找关于这种矩阵的行为的任何注释，在原地和原地 FFT 测试，但我忘记了一些东西。

我在 MATLAB 中用相同的矩阵测试了相同的算法，一切都是正确的。我向您展示了一个非常简化的代码，其中包含一个非常基本的过滤器，以确保清晰、它的输出和预期的输出。我究竟做错了什么？我还阅读了其他相关问题/答案，但没有一个能解决问题。非常感谢您。

const int W = 5;
const int H = 4;

//Signal, with just one value, for simplicity.
float A[H][W] = {{0,0,0,0,0},
        {0,1,0,0,0},
        {0,0,0,0,0},
        {0,0,0,0,0}};

//Central element of the kernel in the (0,0) position of the array.
float B[H][W] = {{0.5, 0.1,  0,    0,  0.2},
        {0  ,  0  ,  0,    0,  0},
        {0  ,  0  ,  0,    0,  0},
        {0  ,  0  ,  0,    0,  0}};


cufftReal* d_inA, *d_inB;
cufftComplex* d_outA, *d_outB;

size_t real_size = W * H * sizeof(cufftReal);
size_t complex_size = W * (H/2+1) * sizeof(cufftComplex);

cudaMalloc((void**)&d_inA, real_size);
cudaMalloc((void**)&d_inB, real_size);

cudaMalloc((void**)&d_outA, complex_size);
cudaMalloc((void**)&d_outB, complex_size);

cudaMemset(d_inA,0,real_size);
cudaMemset(d_inB,0,real_size);

cudaMemcpy(d_inA, A, real_size, cudaMemcpyHostToDevice);
cudaMemcpy(d_inB, B, real_size, cudaMemcpyHostToDevice);


cufftHandle fwplanA, fwplanB, bwplan;
cufftPlan2d(&fwplanA, W, H, CUFFT_R2C);
cufftPlan2d(&fwplanB, W, H, CUFFT_R2C);
cufftPlan2d(&bwplan, W, H, CUFFT_C2R);

cufftSetCompatibilityMode(fwplanA,CUFFT_COMPATIBILITY_NATIVE);
cufftSetCompatibilityMode(fwplanB,CUFFT_COMPATIBILITY_NATIVE);
cufftSetCompatibilityMode(bwplan,CUFFT_COMPATIBILITY_NATIVE);   

cufftExecR2C(fwplanA, d_inA, d_outA);
cufftExecR2C(fwplanB, d_inB, d_outB);

int blocksx = ceil((W*(H/2+1 )) / 256.0f);
dim3 threads(256);
dim3 grid(blocksx);
// One complex product for each thread, scaled by the inverse of the
// number of elements involved in the FFT
pointwise_product<<<grid, threads>>>(d_outA, d_outB, (W*(H/2+1)), 1.0f/(W*H));

cufftExecC2R(bwplan, d_outA, d_inA);


cufftReal* result = new cufftReal[W*2*(H/2+1)];
cudaMemcpy(result, d_inA, real_size,cudaMemcpyDeviceToHost);

// Print result...
// Free memory...

输出。注意位移值

-0.0  0.0   -0.0   -0.0   0.0   
0.0   0.5    0.1    0.0   -0.0  
0.2   0.0    0.0   -0.0   -0.0  
-0.0  0.0   -0.0   -0.0   -0.0

预期输出（MATLAB）

0         0         0         0         0
0.2000    0.5000    0.1000    0.0000    0.0000
0         0         0         0         0
0         0         0         0         0

【问题讨论】：

您是否在袖带计划中与列交换行？原型是cufftPlan2d(cufftHandle *plan, int nx, int ny, cufftType type)，其中nx是行数，ny是列数，所以应该是cufftPlan2d(&fwplanA, H, W, CUFFT_R2C);而不是cufftPlan2d(&fwplanA, W, H, CUFFT_R2C);。
那是错误......我不敢相信我以前没有看到它。我查看了 cufft 文档并说 nx=rows, ny=columns... 非常感谢！
我已为您的问题添加了一个答案，以便将该帖子从未回答列表中删除。

标签： c++ image-processing cuda fft convolution

【解决方案1】：

您正在与袖带计划中的列交换行。