关于FFTW3专家界面的困惑:3个同时进行的复杂FFT

赫尔普·德灵顿

我正在尝试使用FFTW3大师界面但是,文档中的描述使我更加困惑。我想为3个标量复数FFT创建一个计划,每个长度为n,应在内存中一个接一个地对齐(交织数组)。我不明白这是暗示rank == 3还是暗示howmany_rank == 3我也不知道,哪些值的阵列fftw_iodim *dims,并const fftw_iodim *howmany_dims需要有。

尽管这篇stackexchange帖子似乎有一个很好的答案,但是它使用了文档中的相同术语,并且没有“在我的脑海中创建图像”来说明FFT如何在内存中针对不同ranks和不同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位)接口。

弗朗西斯

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的尺寸数。这些尺寸的扩展和布局由arguments指定*dims在执行DFT变换后,这些维度的索引对应于频率。这些参数定义了要应用的多维DFT。

自变量howmany_rankhowmany_dims指定起点,以多次应用上述定义的多维变换。起点可以位于由延伸和步幅描述的多维数组上。参数howmany_rank描述起点数组的维数。

让我们考虑一个具有3个标量分量u,v,w的字段,该标量是在直线上以x_i = iL / N的规则间隔点采样的。该字段以大小为(N,3)的交错连续2D数组存储在内存中:

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=Ndims[0].is=3dims[0].os=3is用于输入步幅,os用于输出步幅。

一维DFT将执行多次,每个组件一次。由于起始这些的DFT的点u(x_0)v(x_0)并且w(x_0)被规则地间隔开,这些起始点的位置限定尺寸1。因此阵列howmany_rank=1此外,由于有3个连续的出发点,出发点的阵列的布局被定义为howmany_dims[0].n=3howmany_dims[0].is=1howmany_dims[0].os=1

我对如何使用fftw Guru界面的回答中提供的示例代码(很抱歉,它对您没有帮助!)可以很容易地调整为:

#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

本文收集自互联网,转载请注明来源。

如有侵权,请联系 [email protected] 删除。

编辑于
0

我来说两句

0 条评论
登录 后参与评论

相关文章