复杂数组的指针数学

Pointer Math with Complex Array

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

#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 个元素,盒子有 320 个。

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

但是,fftwf_complex 被定义为 typedef float fftw_complex[2],因为您需要同时保存复数的实部和虚部。它也被从 fftwf_complex * 转换为浮点数 *,并且我不确定这是做什么的,考虑到 typedef。

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

不同

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

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转换回浮点数组是可行的,因为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]
    }
  }
}