Как мне использовать fftw_plan_many_dft для транспонированного массива данных? - PullRequest
Новая
       25

Как мне использовать fftw_plan_many_dft для транспонированного массива данных?

10 голосов
/ 16 мая 2011

У меня есть двумерный массив данных, сохраненный в формате основных столбцов (в стиле Фортрана), и я хотел бы взять БПФ каждой строки .Я хотел бы избежать транспонирования массива (он не квадратный).Например, мой массив 

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

содержит записи [r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255].

Я могу использовать fftw_plan_many_dft, чтобы составить план для решения каждого из 21 БПФ на месте в 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;
}

Согласно документации ( раздел 4.4.1 руководства FFTW ),подпись для функции 

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).Однако следующий код приводит к сегментированию.ошибка ( редактировать: длина шага неправильно , см. обновление ниже): 

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;
}

1 Ответ

8 голосов
/ 17 мая 2011

Функция работает как описано.Я сделал ошибку с длиной шага, которая в данном случае должна быть howmany.Я обновил вопрос, чтобы отразить это.

Я считаю, что документацию для FFTW довольно сложно понять без примеров (я мог бы быть просто неграмотным ...), поэтому я публикую более подробный пример ниже, сравнивая обычное использование fftw_plan_dft_1d сfftw_plan_many_dft.Напомним, что в случае howmany массивов длиной N, которые хранятся в непрерывном блоке памяти, указанном как in, элементами массива j для каждого преобразования k являются 

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

Следующие два фрагмента кода эквивалентны.В первом преобразование из некоторого 2D-массива выполняется явно, а во втором вызов 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
}

Plan Many 

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);
...