5

我目前正在开发一个必须实现 2D-FFT 的程序（用于互相关）。我用 CUDA 做了一个 1D FFT，它给了我正确的结果，我现在正在尝试实现一个 2D 版本。由于网上的示例和文档很少，我发现很难找出错误所在。

到目前为止，我一直只使用 cuFFT 手册。

无论如何，我创建了两个 5x5 数组并用 1 填充它们。我已将它们复制到 GPU 内存并进行了前向 FFT，将它们相乘，然后对结果进行了 ifft。这给了我一个值为 650 的 5x5 阵列。我希望仅在 5x5 阵列的一个插槽中获得值为 25 的 DC 信号。相反，我在整个数组中得到 650。

此外，在将信号复制到 GPU 内存后，我不允许打印信号的值。写作

cout << d_signal[1].x << endl;

给我一个访问违规。我在其他 cuda 程序中做过同样的事情，这不是问题。它与复杂变量的工作方式有关，还是人为错误？

如果有人对出了什么问题有任何指示，我将不胜感激。这是代码

   #include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <helper_functions.h>
#include <helper_cuda.h>

#include <ctime>
#include <time.h>
#include <stdio.h>
#include <iostream>
#include <math.h>
#include <cufft.h>
#include <fstream>

using namespace std;
typedef float2 Complex;





__global__ void ComplexMUL(Complex *a, Complex *b)
{
    int i = threadIdx.x;
    a[i].x = a[i].x * b[i].x - a[i].y*b[i].y;
    a[i].y = a[i].x * b[i].y + a[i].y*b[i].x;
}


int main()
{


    int N = 5;
    int SIZE = N*N;


    Complex *fg = new Complex[SIZE];
    for (int i = 0; i < SIZE; i++){
        fg[i].x = 1; 
        fg[i].y = 0;
    }
    Complex *fig = new Complex[SIZE];
    for (int i = 0; i < SIZE; i++){
        fig[i].x = 1; // 
        fig[i].y = 0;
    }
    for (int i = 0; i < 24; i=i+5)
    {
        cout << fg[i].x << " " << fg[i + 1].x << " " << fg[i + 2].x << " " << fg[i + 3].x << " " << fg[i + 4].x << endl;
    }
    cout << "----------------" << endl;
    for (int i = 0; i < 24; i = i + 5)
    {
        cout << fig[i].x << " " << fig[i + 1].x << " " << fig[i + 2].x << " " << fig[i + 3].x << " " << fig[i + 4].x << endl;
    }
    cout << "----------------" << endl;

    int mem_size = sizeof(Complex)* SIZE;


    cufftComplex *d_signal;
    checkCudaErrors(cudaMalloc((void **) &d_signal, mem_size)); 
    checkCudaErrors(cudaMemcpy(d_signal, fg, mem_size, cudaMemcpyHostToDevice));

    cufftComplex *d_filter_kernel;
    checkCudaErrors(cudaMalloc((void **)&d_filter_kernel, mem_size));
    checkCudaErrors(cudaMemcpy(d_filter_kernel, fig, mem_size, cudaMemcpyHostToDevice));

    // cout << d_signal[1].x << endl;
    // CUFFT plan
    cufftHandle plan;
    cufftPlan2d(&plan, N, N, CUFFT_C2C);

    // Transform signal and filter
    printf("Transforming signal cufftExecR2C\n");
    cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD);
    cufftExecC2C(plan, (cufftComplex *)d_filter_kernel, (cufftComplex *)d_filter_kernel, CUFFT_FORWARD);

    printf("Launching Complex multiplication<<< >>>\n");
    ComplexMUL <<< 32, 256 >> >(d_signal, d_filter_kernel);

    // Transform signal back
    printf("Transforming signal back cufftExecC2C\n");
    cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE);

    Complex *result = new Complex[SIZE];
    cudaMemcpy(result, d_signal, sizeof(Complex)*SIZE, cudaMemcpyDeviceToHost);

    for (int i = 0; i < SIZE; i=i+5)
    {
        cout << result[i].x << " " << result[i + 1].x << " " << result[i + 2].x << " " << result[i + 3].x << " " << result[i + 4].x << endl;
    }

    delete result, fg, fig;
    cufftDestroy(plan);
    //cufftDestroy(plan2);
    cudaFree(d_signal);
    cudaFree(d_filter_kernel);

}

上面的代码给出了以下终端输出：

1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
----------------
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
----------------
Transforming signal cufftExecR2C
Launching Complex multiplication<<< >>>
Transforming signal back cufftExecC2C

625 625 625 625 625
625 625 625 625 625
625 625 625 625 625
625 625 625 625 625
625 625 625 625 625
4

3 回答 3

2

这给了我一个值为 650 的 5x5 数组：它读取 625，即 5 5 5 5。您使用的卷积算法需要除以 N N。实际上，在cufft中，正向变换中没有归一化系数。因此，您的卷积不能是频域中两个字段的简单相乘。（有些人会称之为数学家 DFT 而不是物理学家 DFT）。

此外，在将信号复制到 GPU 内存后，我不允许打印信号的值：这是标准的 CUDA 行为。在设备上分配内存时，数据存在于设备内存地址空间中，不经额外努力就无法被 CPU 访问。搜索托管内存或零拷贝以从 PCI Express 的两侧访问数据（这在许多其他帖子中都有讨论）。

于 2016-04-27T13:16:00.000 回答
2

这里有几个问题：

  1. 对于乘法内核中输入数组的大小，您启动了太多线程，因此应该会因内存越界错误而失败。我很惊讶您没有收到任何类型的运行时错误。
  2. 我相信，您对 fft/fft - 点积 - ifft 序列的预期解决方案是不正确的。正确的解决方案是一个 5x5 矩阵，每个条目有 25 个。
  3. 正如 cuFFT 文档中明确描述的那样，该库执行非规范化FFT：

    cuFFT 执行未归一化的 FFT；也就是说，对输入数据集执行正向 FFT，然后对结果集执行反向 FFT，得到的数据等于输入，按元素数量缩放。通过数据集大小的倒数缩放任一变换，留给用户执行。

因此，根据我的估计，您的代码的正确输出解决方案应该是一个 5x5 矩阵，每个条目中有 625 个，这将被归一化为每个条目中有 25 个的 5x5 矩阵，即。预期的结果。我不明白 (1) 处的问题如何没有产生不同的结果，因为乘法内核应该失败。

TLDR；这里没什么好看的，离开...

于 2016-04-27T13:16:58.560 回答
0

就像已经提到的其他事情一样：我认为您的复数乘法内核没有做正确的事情。您正在覆盖a[i].x第一行，然后使用a[i].x第二行中的新值来计算a[i].y. a[i].x我认为您需要在覆盖之前先生成一个备份，例如：

float aReal_bk = a[i].x;
a[i].x = a[i].x * b[i].x - a[i].y * b[i].y;
a[i].y = aReal_bk * b[i].y + a[i].y * b[i].x;
于 2021-09-26T06:06:03.720 回答