c - FFTW:3D 数据的 2D 傅里叶变换输出数组的大小

Question

我有一个真正的 Nx*Ny*Nz 维度的 3D 数组，想对每个 z 值进行二维傅里叶变换 using FFTW .这里的 z 索引是内存中变化最快的。目前以下代码按预期工作:

int Nx = 16; int Ny = 8; int Nz = 3;

// allocate memory
const int dims = Nx * Ny * Nz;

// input data (pre Fourier transform)
double *input = fftw_alloc_real(dims);              

// why is this the required output size?
const int outdims = Nx * (Ny/2 + 1) * Nz;           
// we want to perform the transform out of place
// (so seperate array for output)
fftw_complex *output = fftw_alloc_complex(outdims);    


// setup "plans" for forward and backward transforms
const int rank = 2; const int howmany = Nz;
const int istride = Nz; const int ostride = Nz;
const int idist = 1; const int odist = 1;
int n[] = {Nx, Ny};
int *inembed = NULL, *onembed = NULL;


fftw_plan fp = fftw_plan_many_dft_r2c(rank, n, howmany,
        input, inembed, istride, idist,
        output, onembed, ostride, odist,
        FFTW_PATIENT);

fftw_plan bp = fftw_plan_many_dft_c2r(rank, n, howmany,
        output, onembed, ostride, odist,
        input, inembed, istride, idist,
        FFTW_PATIENT);

据我了解，转换长度为 N 的一维序列需要 (N/2 + 1) 个复数值，所以如果我设置 outdims = (Nx/2 + 1)*( Ny/2 + 1)*Nz 正如人们对二维变换所期望的那样？

其次，我认为我可以使用以下方法访问从 qx = 0 到 Nx/2(含)的模式的实部和虚部:

#define outputRe(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * (qx))][0] )
#define outputIm(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * (qx))][1] )

编辑:Full code和 Makefile对于那些想玩的人。假设安装了 fftw 和 gsl。

EDIT2:如果我理解正确，索引(允许正频率和负频率)应该是(对于宏来说可能太乱了!):

#define outputRe(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * ( ((qx) >= 0) ? (qx) : (Nx + (qx)) ) )  ][0] )
#define outputIm(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * ( ((qx) >= 0) ? (qx) : (Nx + (qx)) ) )  ][1] )

for (int qx = -Nx/2; qx < Nx/2; ++qx)
    for (int qy = 0; qy <= Ny/2; ++qy)
        outputRe(qx, qy, d) = ...

其中 outputRe(-Nx/2, qy, d) 指向与 outputRe(Nx/2, qy, d) 相同的数据。在实践中，我可能只是循环遍历第一个索引并转换为频率，而不是反过来!

Answer 1

帮助澄清(关注 2D，因为它很容易扩展到 3D 数据的 2D 转换):

存储

Nx * Ny 数组在傅里叶变换后需要 Nx * (Ny/2 + 1) 个复杂元素。

首先，在 y 方向上，负频率可以从复共轭对称性(来自转换实数序列)中重建。 y 模式 ky 然后从 0 到 Ny/2 运行。 因此对于 y，我们需要 Ny/2 + 1 复数值。

接下来我们在 x 方向上进行变换，因为我们对复值 y 值进行操作，因此我们不能使用相同的对称假设。因此我们必须包括正频率和负频率，因此 x 模式 kx 从 -Nx/2 到 Nx/2 运行。但是 kx = -Nx/2 和 kx = Nx/2 是等价的，所以只存储一个(参见 here )。 因此对于 x 我们需要 Nx 复数值。

访问

正如 tir38 指出的 x 索引后转换从 0 运行到 Nx-1，但这并不意味着模式 kx 从 0 运行到 Nx-1。 FFTW 在数组的前半部分打包正频率，然后在后半部分打包负频率(以相反的顺序)，例如:

kx = 0, 1, 2, ..., Nx/2, -Nx/2 + 1, ..., -2, -1

我们可以考虑两种方式来访问这些元素。首先，正如 tir38 建议的那样，我们可以按顺序循环并从索引中计算出模式 kx:

for (int i = 0; i < Nx; i++)
{
    // produces the list of kxs above
    int kx = (i <= Nx/2) ? i : i - Nx;

    // here we index with i, but with the knowledge that the mode is kx
    outputRe(i, ...) = some function of kx 
}

或者我们可以遍历模式 kx 并转换为索引:

for (int kx = -Nx/2; kx < Nx/2; kx++)
{
    // work out index from mode kx
    int i = (kx >= 0) ? i : Nx + i;

    // here we index with i, but with the knowledge that the mode is kx
    outputRe(i, ...) = some function of kx 
}

两种类型的索引以及其余代码 can found here .