FFTW3大师界面的困惑：3个同时复杂的FFT

Question

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

虽然似乎有一个很好的答案，但它使用了与文档相同的术语，并且没有 "create an image in my head" 不同 rank 的 FFT 如何在内存中对齐s 和不同的 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);

请注意，无法使用 "standard interface"，因为稍后我必须切换到 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 的场，在横坐标 spaced 点 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)

将在 space 坐标 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)有规律地spaced，所以这些起点的位置定义了一个维度为1的数组。因此howmany_rank=1.此外，由于有3个连续的起点，起点数组的布局定义为howmany_dims[0].n=3、howmany_dims[0].is=1和howmany_dims[0].os=1。

我对的回答中提供的示例代码（很抱歉，它没有帮助到您！）可以很容易地改编为：

#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编译。

FFTW3大师界面的困惑：3个同时复杂的FFT

Confusion about FFTW3 guru interface: 3 simultaneous complex FFTs

c++

fftw