我正在尝试使用FFTW3大师界面。但是,文档中的描述使我更加困惑。我想为3个标量复数FFT创建一个计划,每个长度为n
,应在内存中一个接一个地对齐(交织数组)。我不明白这是暗示rank == 3
还是暗示howmany_rank == 3
。我也不知道,哪些值的阵列fftw_iodim *dims
,并const fftw_iodim *howmany_dims
需要有。
尽管这篇stackexchange帖子似乎有一个很好的答案,但是它使用了文档中的相同术语,并且没有“在我的脑海中创建图像”来说明FFT如何在内存中针对不同rank
s和不同howmany_rank
s对齐。也就是说,答案无济于事。
供参考,以下是计划分配功能的定义:
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_rank
并howmany_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=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
。
我对如何使用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] 删除。
我来说两句