3

我曾经用于fftw_plan_dft多维傅里叶变换。

fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in,
                        fftw_complex *out, int sign, unsigned flags);

现在我想将 64 位整数传递给 fftw，看来我需要使用 fftw guru 接口。

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

但我不明白是什么howmany_rank意思howmany_dims。的手册fftw_plan_guru_dft说：

这两个函数分别为交错和拆分格式规划了一个复杂数据、多维 DFT。变换维度由（rank, dims）在维度（howmany_rank, howmany_dims）的多维向量（循环）上给出。dims 和 howmany_dims 应该分别指向长度为 rank 和 howmany_rank 的 fftw_iodim 数组。

我知道什么是“维度的多维向量（循环）（howmany_rank，howmany_dims）”是什么意思。你能给我一个例子或解释一下如何使用这个大师界面吗？

4

1 回答 1

2

如果多维数组的大小和步幅大于 2^32，则64 位 guru 接口将变得有用。

创建复杂到复杂 DTF 的函数原型是：

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

在哪里：

  • rank是要执行的 FFTW 变换的秩，即维数。
  • dims是一个大小数组rank。对于每个维度i，dims[i].n是线的大小，dims[i].is是输入数组dims[i].os的线之间的步幅， 是输出数组的线之间的步幅。例如，如果数组在内存中是连续的，那么guru 接口的文档建议使用递归dims[i].is = n[i+1] * dims[i+1].is。要执行的变换次数和起点之间的偏移量由howmany_rank和给出howmany_dims。
  • howmany_rank指定具有特定偏移的变换的数量。
  • howmany_dims是一个大小数组howmany_rank。对于每个变换i，howmany_dims[i].n是要计算的变换的数量，每个变换都具有输入 之间的howmany_dims[i].is偏移量和输出之间的偏移量howmany_dims[i].os。

有关这些论点的更多详细信息，请参见关于 FFTW3 大师界面的混淆：3 个同时复杂的 FFT

以下代码调用fftw_plan_guru64_dft()，以便它执行与fftw_plan_dft_3d(). 它可以通过以下方式编译gcc main.c -o main -lfftw3 -lm -Wall：

#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 M = 12;
    unsigned long int P = 14;
    fftw_complex *in=fftw_malloc(N*M*P*sizeof(fftw_complex));
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    fftw_complex *out=fftw_malloc(N*M*P*sizeof(fftw_complex));
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    unsigned int i,j,k;

    int rank=3;
    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=P*M;
    dims[0].os=P*M;
    dims[1].n=M;
    dims[1].is=P;
    dims[1].os=P;
    dims[2].n=P;
    dims[2].is=1;
    dims[2].os=1;

    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=1;
    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++){
        for(j=0;j<M;j++){
            for(k=0;k<P;k++){
                //printf("ijk\n");
                in[(i*M+j)*P+k]=30.+12.*sin(2*3.1415926535*i/((double)N))*sin(2*3.1415926535*j/((double)M))*sin(2*3.1415926535*k/((double)P))*I;
            }
        }
    }

    fftw_execute(p);

    for (i = 0; i < N; i++){
        for (j = 0; j < M; j++){
            for (k = 0; k < P; k++){
                printf("result: %d %d %d %g %gI\n", i,j,k, creal(out[(i*M+j)*P+k]), cimag(out[(i*M+j)*P+k]));
            }
        }
    }


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

    free(dims);
    free(howmany_dims);

    return(0);
}

例如，guru 界面可用于计算复杂 3D 电场的 DFT。在网格的每个点，电场都是大小为 3 的向量。因此，我可以将电场存储为 4D 数组，最后一个维度指定向量的分量。最后，guru 界面可用于一次执行三个 3D DFT：

#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 M = 12;
    unsigned long int P = 14;
    unsigned long int DOF = 3;
    fftw_complex *in=fftw_malloc(N*M*P*DOF*sizeof(fftw_complex));
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    fftw_complex *out=fftw_malloc(N*M*P*DOF*sizeof(fftw_complex));
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    unsigned int i,j,k;

    int rank=3;
    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=P*M*DOF;
    dims[0].os=P*M*DOF;
    dims[1].n=M;
    dims[1].is=P*DOF;
    dims[1].os=P*DOF;
    dims[2].n=P;
    dims[2].is=DOF;
    dims[2].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++){
        for(j=0;j<M;j++){
            for(k=0;k<P;k++){
                //printf("ijk\n");
                in[((i*M+j)*P+k)*DOF]=30.+12.*sin(2*3.1415926535*i/((double)N))*sin(2*3.1415926535*j/((double)M))*sin(2*3.1415926535*k/((double)P))*I;
                in[((i*M+j)*P+k)*DOF+1]=42.0;
                in[((i*M+j)*P+k)*DOF+2]=1.0;
            }
        }
    }

    fftw_execute(p);

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


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

    free(dims);
    free(howmany_dims);

    return(0);
}
于 2016-10-14T17:16:22.070 回答