c++ - FFTW3、cuFFT 和就地变换

Question

我正在尝试使用 FFTW3 进行实数到复数的 FFT 变换。到目前为止，我已经设法使用异地转换来完成它，但是我在实现它的就地版本时遇到了麻烦。我的印象是，就地转换你唯一需要改变的是:1)确保你的数据数组有足够的空间来容纳操作的复杂部分，2)当你创建计划时使用相同的输入和输出数据的地址，3)执行计划时，对输入和输出数据使用相同的地址。我已经做了所有这些事情，但我总是得到错误的结果。我正在对值为 [[1,1],[1,1]] 的 2x2 数组执行二维 FFT。预期结果(根据 Matlab)是一个 2x2 数组，其值为 [[4+0i, 0+0i], [0+0i, 0+0i]]。当我进行异地变换时，我得到了这个结果。但是当我进行就地转换时，我得到以下 [[2+0i, 0+0i],[2+0i,0+0i]]。我选择了 2x2 大小的 2D FFT，因为输入和输出数据的长度相同并且有助于调试。这是我的代码:

bool inplace = true; // true for in-place, false for out-of-place
int dim_size[] = {2,2};
int N[] = {2,2};
int data_length     = N[0]*(N[1]);      //  2 * (2)     = 4
int data_fft_length = N[0]*(N[1]/2+1);  //  2 * (2/2+1) = 4
float* h_data_r = nullptr;              //  fftw data array
fftwf_complex* h_data_c = nullptr;      //  fftw data array (only used in out-of-place tranforms)

//  allocate fftw memory
if(inplace) {
    h_data_r = (float*)fftwf_malloc(data_fft_length*sizeof(fftwf_complex));
    h_data_c = (fftwf_complex*)h_data_r;
} else {
    h_data_r = (float*)fftwf_malloc(data_length*sizeof(float));
    h_data_c = (fftwf_complex*)fftwf_malloc(data_fft_length*sizeof(fftwf_complex));
}

//  create plane
unsigned int flags = FFTW_MEASURE;
fftwf_plan m_plan = fftwf_plan_dft_r2c_2d(N[0],N[1],h_data_r,h_data_c,flags);

//  initialize data array
h_data_r[0] = 1;
h_data_r[1] = 1;
h_data_r[2] = 1;
h_data_r[3] = 1;

//  execute fft plan
fftwf_execute(m_plan);

std::cout << "result:" << std::endl;
for(int i = 0; i < data_fft_length; ++i)
    std::cout << "[" << i << "]: " << h_data_c[i][0] << " " << h_data_c[i][1]  << std::endl;

变量“就地”决定 FFT 变换是否就地。有人能告诉我我的代码有什么问题吗？代码非常简单。我没有做任何特别的事情。我只想要一个就地 FFT 实数到复数变换。如果您不愿意检查我的代码，但您有一个非常简单的 fftw3 就地 fft 转换代码，请随时复制粘贴它。

谢谢。

编辑 1:我进一步简化了代码，现在我使用 fftwf_plan_dft_r2c_2d() 来创建计划，使用 fftwf_execute() 来执行计划。问题依旧。

编辑 2:我将代码翻译成 cufft，它应该与 fftw3 具有几乎相同的语法。我对袖口有同样的问题。但默认情况下，cuFFT 启用了 FFTW 兼容模式 (CUFFT_COMPATIBILITY_FFTW_PADDING)。如果我使用标志 CUFFT_COMPATIBILITY_NATIVE 禁用 FFTW 兼容模式，则就地转换与 cuFFT 一起工作得很好。奇怪的是，根据 cuFFT 文档，CUFFT_COMPATIBILITY_FFTW_PADDING 应该在您进行批量转换时有所作为。就我而言，我不进行任何批量转换。我现在更糊涂了。

Answer 1

几个月前我设法解决了这个问题。结果我不得不使用高级计划创建界面并手动设置 inembed 和 onembed 指针的值。