c++ - 关于 FFTW3 大师界面 : 3 simultaneous complex FFTs 的困惑

Question

我正在尝试使用 FFTW3 guru interface .但是，文档中的描述使我更加困惑而不是清晰。我想为 3 个标量复数 FFT 创建一个计划，每个长度 n，它们应该在内存中一个接一个地对齐(交错数组)。我不明白这是否意味着 rank == 3 或 howmany_rank == 3。我也不知道，数组 fftw_iodim *dims 和 const fftw_iodim *howmany_dims 需要具有哪个值。

同时 this stackexchange post似乎有一个很好的答案，它使用与文档相同的术语并且不会“在我的脑海中创建一个图像”关于 FFT 如何在内存中针对不同的 rank 和不同的 howmany_ranks。也就是说，答案对我要查找的内容没有帮助。

供引用，这是计划分配函数的定义:

fftw_plan fftw_plan_guru_dft(
      int rank, const fftw_iodim *dims,
      int howmany_rank, const fftw_iodim *howmany_dims,
      fftw_complex *in, fftw_complex *out,
      int sign, unsigned flags);

请注意，使用“标准界面”是不可能的，因为稍后我必须切换到 guru64(64 位)界面。

Answer 1

fftw_plan fftw_plan_guru_dft() 允许定义 1D-2D-3D-4D... DFT(rank=1,2,3,4...)应用多次使用位于 1D、2D、3D、4D... 网格上的起点(howmany_rank=1,2,3,4....

参数 rank 指定应用 DFT 的维数。这些维度的扩展和布局由参数 *dims 指定。执行DFT变换后，这些维度的索引对应于频率。这些参数定义了要应用的多维 DFT。

参数 howmany_rank 和 howmany_dims 指定多次应用上面定义的多维变换的起点。起点可以位于由 extends 和 strides 描述的多维数组上。参数 howmany_rank 描述起点数组的维数。

让我们考虑一个具有 3 个标量分量 u、v、w 的场，在横坐标 x_i=iL/N 的规则间隔点的直线上采样。该字段作为维度 (N,3) 的交错连续二维数组存储在内存中:

u(x_0) v(x_0) w(x_0) u(x_1) v(x_1) w(x_1) u(x_2) v(x_2) w(x_2)... u(x_N-1) v(x_N-1) w(x_N-1) 

将在空间坐标 x 上执行一维 DFT。因此，rank 等于 1。每个 1D DFT 的长度为 N，并且两个连续项(u_(x_0) 和 u( x_1)) 为 3。因此，dims[0].n=N、dims[0].is=3 和 dims[0 ].os=3 。 is 用于输入步幅，os 用于输出步幅。

一维 DFT 将执行多次，每个组件执行一次。由于这些 DFT 的起点 u(x_0)、v(x_0) 和 w(x_0) 是规则间隔的，因此这些位置起点定义一个维度为 1 的数组。因此 howmany_rank=1。此外，由于有 3 个连续的起点，起点数组的布局定义为 howmany_dims[0].n=3, howmany_dims[0].is=1 和 howmany_dims[0].os=1。

我对 How to use fftw Guru interface 的回答中提供的示例代码(很抱歉，它没有帮助你!)可以很容易地改编为:

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

int main(void){

    fftw_plan p;
    unsigned long int N = 10;
    unsigned long int DOF = 3;
    fftw_complex *in=fftw_malloc(N*DOF*sizeof(fftw_complex));
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    fftw_complex *out=fftw_malloc(N*DOF*sizeof(fftw_complex));
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    unsigned int i;

    int rank=1;
    fftw_iodim64 *dims=malloc(rank*sizeof(fftw_iodim64));
    if(dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    dims[0].n=N;
    dims[0].is=DOF;
    dims[0].os=DOF;

    int howmany_rank=1;
    fftw_iodim64 *howmany_dims=malloc(howmany_rank*sizeof(fftw_iodim64));
    if(howmany_dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    howmany_dims[0].n=DOF;
    howmany_dims[0].is=1;
    howmany_dims[0].os=1;

    printf("sizeof fftw complex %ld\n",sizeof(fftw_complex));
    printf("sizeof fftw_iodim64 %ld\n",sizeof(fftw_iodim64));
    printf("creating the plan\n");
    p=fftw_plan_guru64_dft(rank, dims,howmany_rank, howmany_dims,in, out,FFTW_FORWARD, FFTW_ESTIMATE);
    if (p==NULL){fprintf(stderr,"plan creation failed\n");exit(1);}
    printf("created the plan\n");

    for(i=0;i<N;i++){
                in[i*DOF]=30.+12.*sin(2*3.1415926535*i/((double)N));
                in[i*DOF+1]=42.0;
                in[i*DOF+2]=1.0;
    }

    fftw_execute(p);

    for (i = 0; i < N; i++){
        printf("result: %d || %g %gI | %g %gI | %g %gI\n", i, creal(out[i*DOF]), cimag(out[i*DOF]),creal(out[i*DOF+1]), cimag(out[i*DOF+1]),creal(out[i*DOF+2]), cimag(out[i*DOF+2]));
    }


    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

    free(dims);
    free(howmany_dims);

    return(0);
}

并由gcc main.c -o main -lfftw3 -lm -Wall编译。