如何在转置的数据数组上使用fftw_plan_many_dft？
我有一个以列主要（Fortran风格）格式存储的2D数据数组，我想对每一行进行FFT。 我想避免转置数组（它不是正方形）。 例如，我的数组
fftw_complex* data = new fftw_complex[21*256]; 包含条目[r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255] 。
我可以使用fftw_plan_many_dft制定计划，以便在dataarrays中就地解决每个21个FFT，如果它是行 –主要的，例如[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节 ），该function的签名是
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和dist参数来设置索引。 根据我从文档中可以理解的，要转换的数组中的条目被索引为in + j*istride + k*idist ，其中j=0..n-1且k=0..howmany-1 。 （我的数组是1D，并且它们有howmany多）。 但是，以下代码会产生一个seg。 故障（ 编辑：步幅长度错误 ，请参阅下面的更新）：
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; }
该function按文档记录。 我用步幅长度做了一个错误，在这种情况下实际应该是howmany 。 我已更新问题以反映这一点。
我发现FFTW的文档在没有示例的情况下有点难以理解（我可能只是文盲……），所以我在下面发布一个更详细的例子，比较fftw_plan_dft_1d和fftw_plan_many_dft的常用用法。 回顾一下，在长度为N数组存储在连续的存储器块中的情况下，每个变换k的数组元素j是
*(in + j*istride + k*idist)
以下两段代码是等效的。 在第一个中，从一些2Darrays的转换是明确完成的，而在第二个中， fftw_plan_many_dft调用用于完成就地的所有操作。
明确复制
int N, howmany; // ... fftw_complex* data = (fftw_complex*) fftw_malloc(N*howmany*sizeof(fftw_complex)); // ... load data with howmany arrays of length N int istride, idist; // ... if data is column-major, set istride=howmany, idist=1 // if data is row-major, set istride=1, idist=N fftw_complex* in = (fftw_complex*) fftw_malloc(N*sizeof(fftw_complex)); fftw_complex* out = (fftw_complex*) fftw_malloc(N*sizeof(fftw_complex)); fftw_plan p = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_MEASURE); for (int k = 0; k < howmany; ++k) { for (int j = 0; j < N; ++j) { in[j] = data[j*istride + k*idist]; } fftw_execute(p); // do something with out }
计划很多
int N, howmany; // ... fftw_complex* data = (fftw_complex*) fftw_malloc(N*howmany*sizeof(fftw_complex)); // ... load data with howmany arrays of length N int istride, idist; // ... if data is column-major, set istride=howmany, idist=1 // if data is row-major, set istride=1, idist=N fftw_plan p = fftw_plan_many_dft(1,&N,howmany,data,NULL,howmany,1,data,NULL,howmany,1,FFTW_FORWARD,FFTW_MEASURE); fftw_execute(p);