【问题标题】:Pointer Math with Complex Array复杂数组的指针数学
【发布时间】:2017-10-25 17:33:32
【问题描述】:

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

#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] = 某个数字。现在,我明白了,对于普通类型的数组,比如浮点数,数组 + 整数将给出数组 [整数],假设整数在数组的范围内。

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

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

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

【问题讨论】:

  • Cast 的优先级高于 +,因此加法将乘以 sizeof float。
  • 如果 R_FFT_INDEX 索引应该是相对于指针 boxfftwf_complex 元素的索引,则首先通过索引取消引用该元素,然后从结果中提取所需的组件fftwf_complex 对象。这是最简单的方法,实际上是正确的。
  • R_FFT_INDEX 是一个相对于 box 的指针,我认为,但它的构建方式是为了访问转换后存储在 box 中的真实数据。我仍然对它是如何工作的感到困惑。 @stark我不确定我是否理解你所说的加法将乘以sizeof(float)是什么意思?这是否意味着它实际上是元素 4 * 45(假设浮点数为 4 个字节)?
  • 当您将 45 添加到指向浮点数的指针时,您会将 45 * sizeof(float) 添加到指针中的数字。指针算法始终基于指向的对象的类型,这就是为什么 ++ 适用于遍历数组的指针。
  • @stark 我明白了。因此,由于 fftwf_complex 被定义为 typedef float fftwf_complex[2],其中 [0] 是实部,[1] 是虚部,当它被转换为 float * 时,*(float *)box + 45 真的会距离第一个地址 45 * 4 个字节?有没有办法告诉盒子的哪个元素对应,或者在这种情况下这个问题甚至没有意义?

标签: c arrays fftw


【解决方案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 转换回浮点数组是很实用的,因为 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),其中位于 k>N-1 的位置 i、j、k 的项目是没有意义的。见FFTW's Real-data DFT Array Format:

对于就地转换,由于复杂数据略大于实际数据,因此会出现一些复杂情况。在这种情况下,必须用额外的值填充真实数据的最终维度以适应复杂数据的大小——如果最后一个维度是偶数,则增加两个,如果是奇数,则增加一个。也就是说,真实数据的最后一维物理上必须包含 2 * (nd-1/2+1) 个双精度值(正好足以容纳复杂数据)。然而,这个物理数组大小不会改变逻辑数组大小——只有 nd-1 个值实际存储在最后一个维度中,而 nd-1 是传递给规划器的最后一个维度。

这就是引入真正的数组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]
    }
  }
}

【讨论】:

    猜你喜欢
    • 2020-03-04
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2012-11-17
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多