5

我有一个维度为 Nx*Ny*Nz 的真实 3D 数组,并希望使用 FFTW对每个 z 值进行 2D 傅里叶变换。这里 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);

outdims = (Nx/2 + 1)*(Ny/2 + 1)*Nz据我了解,转换长度为 N 的 1D 序列需要 (N/2 + 1) 个复数值,那么如果我设置为 2D 转换所期望的那样,为什么上面的代码会崩溃?

其次,我认为我可以qx = 0 to 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] )

编辑:为那些想要玩的人提供完整的代码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)。在实践中,我可能只是循环第一个索引并转换为频率,而不是反过来!

4

2 回答 2

5

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

贮存

Nx * Ny数组需要Nx * (Ny / 2 + 1)经过傅里叶变换后的复数元素。

首先,在 y 方向上,负频率可以从复共轭对称性(来自变换实数序列)中重建。然后 y 模式ky0 to Ny/2包容开始运行。所以对于 y 我们需要Ny/2 + 1复杂的值。

接下来,我们在 x 方向进行变换,我们不能使用相同的对称假设,因为我们正在对复值 y 值进行操作。因此,我们必须包括正频率和负频率,因此 x 模式kx-Nx/2 to Nx/2包容开始。但是kx = -Nx/2kx = Nx/2是等价的,所以只存储一个(见这里)。所以对于 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 
}

可以在此处找到这两种类型的索引以及其余代码。

于 2013-07-03T08:48:31.630 回答
1

第一个问题:

通过对每个 z 值进行 2D FFT,您基本上是在进行 Nz 数的 2D FFT。对称性仅在一维上。因此,对于单个 [Nx x Ny] 2D FFT,您将获得 Nx * (Ny/2 +1) 的输出。因此,由于您正在执行这些 FFT 的 Nz,您将拥有 [Nx x (Ny/2 +1+) x Nz] 3D 输出。

第二个问题

是的,这应该是输出的存储方式。我知道不是每个人都可以使用 Matlab,但是当我开始使用 FFTW 时,我总是将小矩阵与 C++(或你的情况下的 C)和 Matlab 进行比较

更新(基于您的第二次编辑)

它们仍然从零开始索引,因此:

for (int qx = -Nx/2; qx < Nx/2; ++qx) => for (int gx = 0; gx < Nx; gx++)

对称性在 y 轴上,因此 x 数据不是多余的:

output(0, a, b) != output(Nx-1, a ,b)

其中 a, b 是一些 y 和 z 值;

于 2013-06-28T02:16:13.973 回答