¿Cómo uso fftw_plan_many_dft en una matriz de datos transpuesta?

Question

Tengo una matriz 2D de datos almacenada en formato de columna principal (estilo Fortran) y me gustaría tomar la FFT de cada una. fila.Me gustaría evitar la transposición de la matriz (no es cuadrada).Por ejemplo, mi matriz

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

contiene entradas [r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255].

Puedo usar fftw_plan_many_dft hacer un plan para resolver cada una de las 21 FFT existentes en el data matriz si es fila-mayor, p.e. [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;
}

Según la documentación (sección 4.4.1 del manual FFTW), la firma de la función es

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

y debería poder usar el stride y dist parámetros para configurar la indexación.Por lo que puedo entender de la documentación, las entradas en la matriz a transformar están indexadas como in + j*istride + k*idist dónde j=0..n-1 y k=0..howmany-1.(Mis matrices son 1D y hay howmany de ellos).Sin embargo, el siguiente código da como resultado un seg.falla (editar: la longitud de la zancada es equivocado, ver actualización a continuación):

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

Actualizar:

Cometí un error al elegir la longitud de la zancada.La llamada correcta es (la longitud de zancada correcta es howmany, no 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

La función funciona según lo documentado.Cometí un error con la longitud de la zancada, que en realidad debería ser howmany en este caso.He actualizado la pregunta para reflejar esto.

Creo que la documentación de FFTW es algo difícil de comprender sin ejemplos (podría ser analfabeto...), así que publicaré un ejemplo más detallado a continuación, comparando el uso habitual de fftw_plan_dft_1d con fftw_plan_many_dft.En resumen, en el caso de howmany matrices con longitud N que se almacenan en un bloque contiguo de memoria al que se hace referencia como in, los elementos de la matriz j para cada transformación k son

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

Los siguientes dos fragmentos de código son equivalentes.En el primero, la conversión de alguna matriz 2D se realiza explícitamente, y en el segundo la fftw_plan_many_dft La llamada se utiliza para hacer todo en el lugar.

Copia explícita

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
}

Planifica muchos

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