我有一个以列为主(Fortran 风格)格式存储的二维数据数组,我想对每一行进行 FFT 。我想避免转置数组(它不是方形的)。例如,我的数组
fftw_complex* data = new fftw_complex[21*256];
包含条目[r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255]
。
如果它是行主要的,我可以fftw_plan_many_dft
用来制定计划来解决阵列中的 21 个 FFT 中的每一个,例如:data
[r0_val0, r0_val1,..., r0_val255, r1_val0,...,r20_val255]
int main() {
int N = 256;
int howmany = 21;
fftw_complex* data = new fftw_complex[N*howmany];
fftw_plan p;
// this plan is OK
p = fftw_plan_many_dft(1,&N,howmany,data,NULL,1,N,data,NULL,1,N,FFTW_FORWARD,FFTW_MEASURE);
// do stuff...
return 0;
}
根据文档(FFTW 手册的第 4.4.1 节),该函数的签名是
fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
fftw_complex *in, const int *inembed,
int istride, int idist,
fftw_complex *out, const int *onembed,
int ostride, int odist,
int sign, unsigned flags);
我应该能够使用stride
anddist
参数来设置索引。从我可以从文档中了解到,要转换的数组中的条目被索引为in + j*istride + k*idist
wherej=0..n-1
和k=0..howmany-1
. (我的数组是一维的,并且有howmany
)。但是,以下代码会导致段。错误(编辑:步幅错误,请参阅下面的更新):
int main() {
int N = 256;
int howmany = 21;
fftw_complex* data = new fftw_complex[N*howmany];
fftw_plan p;
// this call results in a seg. fault.
p = fftw_plan_many_dft(1,&N,howmany,data,NULL,N,1,data,NULL,N,1,FFTW_FORWARD,FFTW_MEASURE);
return 0;
}
更新:
我在选择步幅时出错了。正确的调用是(正确的步幅是howmany
,不是N
):
int main() {
int N = 256;
int howmany = 21;
fftw_complex* data = new fftw_complex[N*howmany];
fftw_plan p;
// OK
p = fftw_plan_many_dft(1,&N,howmany,data,NULL,howmany,1,data,NULL,howmany,1,FFTW_FORWARD,FFTW_MEASURE);
// do stuff
return 0;
}