c - 复杂数组的指针数学

Question

我有这段代码，其中包含一些我无法理解的指针数学:

#include <stdlib.h>
#include <complex.h>
#include <fftw3.h>

int main(void)
{
    int i, j, k;
    int N, N2;
    fftwf_complex *box;
    fftwf_plan plan;
    float *smoothed_box;

    // Allocate memory for arrays (Ns are set elsewhere and properly,
    // I've just left it out for clarity)
    box = (fftwf_complex *)fftwf_malloc(N * sizeof(fftwf_complex));
    smoothed_box = (float *)malloc(N2 * sizeof(float));

    // Create complex data and fill box with it. Do FFT. Box has the
    // Hermitian symmetry that complex data has when doing FFTs with 
    // real data
    plan = fftwf_plan_dft_c2r_3d(N,N,N,box,(float *)box,
         FFTW_ESTIMATE);
    ...

    // end fft

    // Now do the loop I don't understand
    for(i = 0; i < N2; i++)
    {
        for(j = 0; j < N2; j++)
        {
            for(k = 0; k < N2; k++)
            {
                smoothed_box[R_INDEX(i,j,k)] = *((float *)box + 
                    R_FFT_INDEX(i*f + 0.5, j*f + 0.5, k*f +0.5))/V;
            }
        }
    }

    // Do other stuff
    ...

    return 0;
}

其中 f 和 V 只是代码中其他地方设置的一些数字，与这个特定问题无关。此外，函数 R_FFT_INDEX 和 R_INDEX 也无关紧要。重要的是，对于第一次循环迭代，当 i=j=k=0 时，R_INDEX = 0 且 R_FFT_INDEX=45。 smoothed_box 有 8 个元素，box 有 320 个。

因此，在 gdb 中，当我在循环后打印 smoothed_box[0] 时，我得到 smoothed_box[0] = 一些数字。现在，我明白了，对于普通类型的数组，比如 float ，数组 + 整数将给出数组[整数]，假设整数在数组的范围内。

但是，fftwf_complex 被定义为 typedef float fftw_complex[2]，因为您需要同时保存复数的实部和虚部。它也被从 fftwf_complex * 转换为 float *，鉴于 typedef，我不确定这是做什么的。

我所知道的是，当我在 gdb 中打印 box[45] 时，我得到 box[45] = 一些不是 smoothed_box[0] * V 的复数。即使我打印 *((float *)box + 45)/V，我得到的数字与 smoothed_box[0] 不同。

所以，我只是想知道是否有人可以向我解释在上述循环中完成的指针数学运算？谢谢，感谢您的宝贵时间!

Answer 1

box 被分配为 N fftwf_complex 的数组。然后在 box 上执行使用 N,N,N 的后向 3D c2r fftw 变换，需要 N*N*(N/2+1) fftwf_complex。参见 http://www.fftw.org/fftw3_doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format因此，此代码可能会在到达指针算法之前触发未定义的行为，例如段错误...

将 box 转换回 float 数组是可行的，因为 DFT 是就地执行的。事实上，box 在创建 fftwf_plan 时被使用了两次。 box既是复数的输入数组，也是实数的输出数组:

plan = fftwf_plan_dft_c2r_3d(N,N,N,box,(float *)box,
     FFTW_ESTIMATE);

一旦 fftwf_execute(plan); 被调用，box 最好被视为一个实数数组。尽管如此，此数组的大小为 N*N*2*(N/2+1)，其中位于位置 i、j、k 且 k>N-1 的项是无意义的。参见 FFTW's Real-data DFT Array Format :

For an in-place transform, some complications arise since the complex data is slightly larger than the real data. In this case, the final dimension of the real data must be padded with extra values to accommodate the size of the complex data—two extra if the last dimension is even and one if it is odd. That is, the last dimension of the real data must physically contain 2 * (nd-1/2+1) double values (exactly enough to hold the complex data). This physical array size does not, however, change the logical array size—only nd-1 values are actually stored in the last dimension, and nd-1 is the last dimension passed to the planner.

这就是为什么要引入真正的数组 smoothed_box 的原因，尽管应该是 N*N*N 数组。如果 smoothed_box 是一个大小为 N*N*N 的数组，则可以执行以下转换:

for(i=0;i<N;i++){
  for(j=0;j<N;j++){
    for(k=0;k<N;k++){
     smoothed_box[(i*N+j)*N+k]=((float *)box)[(i*N+j)*(2*(N/2+1))+k]
    }
  }
}