如何在转置数据数组上使用 fftw_plan_many_dft？

Question

我有一个以列主（Fortran 样式）格式存储的二维数据数组，我想对每行进行 FFT。我想避免转置数组（它不是方形的）。例如，我的数组

fftw_complex* data = new fftw_complex[21*256];

包含条目 [r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255]。

我可以使用 fftw_plan_many_dft 制定计划，以解决 data 数组中 21 个 FFT 中的每一个（如果它是行-major），例如 [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 和 dist 参数来设置索引。据我从文档中可以理解，要转换的数组中的条目索引为 in + j*istride + k*idist 其中 j=0..n-1代码>和k=0..howmany-1。 （我的数组是一维的，有多少）。但是，以下代码会产生一个 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;
}

Answer 1

该功能按文档工作。我在步幅长度方面犯了一个错误，在本例中实际上应该是 howmany。我已经更新了问题以反映这一点。

我发现如果没有示例，FFTW 的文档有点难以理解（我可能只是文盲......），因此我在下面发布了一个更详细的示例，将 fftw_plan_dft_1d 的常用用法与 <代码>fftw_plan_many_dft。回顾一下，如果有 howmany 个长度为 N 的数组存储在引用为 in 的连续内存块中，则数组元素 < code>j 对于每个变换 k 都是

*(in + j*istride + k*idist)

以下两段代码是等效的。在第一个中，从某些二维数组的转换是显式完成的，在第二个中，使用 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);