如何使用FFTW3库执行原位FFT ?

Question

我想使用FFTW3库执行3D数组的r2c和c2r就地FFT。数组的大小为(Nx,Ny,Nz+2)，刚好可以存储NxNy(Nz/2+1)个复数。我希望得到的复数的实数(Re)和虚数(Im)部分连续存储在数组中，即Re, Im, Re, Im, Re, Im，…，等等。

遵循FFTW手册，我猜效率不高，我想出了以下代码:

#include <iostream>
#include <fftw3.h>
int main()
{
const int Nx {4};
const int Ny {4};
const int Nz {4};
const int n[] = {Nx, Ny, Nz};
const int inembed[] = {Nx, Ny, Nz+2};
const int onembed[] = {Nx, Ny, Nz/2+1};
double tmp[] {0, 1, 2, 3}; // Some garbage input
double* in = (double*) fftw_malloc( sizeof(double) * Nx*Ny*(Nz+2));
fftw_complex* out = reinterpret_cast<fftw_complex*>(in);
fftw_plan r2c_3D_a = fftw_plan_many_dft_r2c(3, n, 1, in, inembed, 1, 0, out, onembed, 1, 0, FFTW_MEASURE);
fftw_plan c2r_3D_a = fftw_plan_many_dft_c2r(3, n, 1, out, onembed, 1, 0, in, inembed, 1, 0, FFTW_MEASURE);
// Input array 
for (int i = 0; i < Nx; i++){
for (int j = 0; j < Ny; j++){
for (int k = 0; k < Nz+2; k++){
unsigned int offset;
offset = k + (Nz+2) * (j + i * Ny);
in[offset] = k < Nz ? tmp[k] : 0;
std::cout << in[offset] << " ";
}
std::cout << std::endl;
}
}
// R2C FFT
fftw_execute_dft_r2c(r2c_3D_a, in, out);
// Normalization
for (int i = 0; i < Nx*Ny*(Nz+2); i++){
in[i] = in[i] / (Nx*Ny*Nz);
}

// C2R FFT
//fftw_execute_dft_c2r(c2r_3D_a, out, in);

std::cout << std::endl;
for (int i = 0; i < Nx; i++){
for (int j = 0; j < Ny; j++){
for (int k = 0; k < Nz+2; k++){
unsigned int offset;
offset = k + (Nz+2) * (j + i * Ny);
std::cout << in[offset] << " ";
}
std::cout << std::endl;
}
}
fftw_destroy_plan(r2c_3D_a);
fftw_destroy_plan(c2r_3D_a);
free(in);
return 0;
}

给我期望的输入:

0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0 
0 1 2 3 0 0

在r2c FFT之后，每一行应该包含Nz/2+1 = 3个复数，存储为Re, Im, Re, Im, Re, Im。相反，我得到了这个意想不到的输出:

1.5 0 -0.5 0.5 -0.5 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0 
0 0 0 0 0 0

然而，FFT对，即r2c之后的c2r变换重建了初始数据:

0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0 
0 1 2 3 -0.5 0

注意，最后两列只是填充，在转换之后是不使用的。老实说，我卡住了，有人能帮我一下吗?

Answer 1

您的"意外输出"是预期的。你的输入数据在两个维度上是恒定的，所以输出频率在这些维度上都是零，除了直流分量。如果您更改为可变数据，您将得到非零结果。