11

列優先 (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、配列内の 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パラメータとdistパラメータを使用してインデックスを設定できるはずです。ドキュメントから理解できることから、変換される配列内のエントリはin + j*istride + k*idistwherej=0..n-1およびとしてインデックス付けされますk=0..howmany-1。（私の配列は1Dで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;
}
4

1 に答える 1

9

関数は文書化されているように機能します。ストライドの長さでエラーが発生しました。これは実際にhowmanyはこの場合のはずです。これを反映するように質問を更新しました。

FFTWのドキュメントは、例がないと理解するのがやや難しいと思います（私は読み書きができない可能性があります...）。そこで、の通常の使用法を比較して、より詳細な例を以下に投稿fftw_plan_dft_1dしfftw_plan_many_dftます。要約すると、として参照される連続したメモリブロックに格納されているhowmany長さの配列の場合、各変換の配列要素は次のようになります。Ninjk

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

次の2つのコードは同等です。最初の例では、いくつかの2D配列からの変換が明示的に行われ、2番目の例では、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);
于 2011-05-17T17:27:01.027 に答える