fftw 다차원 분석

libttfw 에서는 2차원 3차원 DFT를 지원하는데, 내가 생각하는 그 2차원 3차원이 맞나 모르겠다.

센서가 3축이면 3차원은 맞긴한데.. 희소행렬 형태로 삽입이 되어야 하나?

2.2 Complex Multi-Dimensional DFTs Multi-dimensional transforms work much the same way as one-dimensional transforms: you allocate arrays of fftw_complex (preferably using fftw_malloc), create an fftw_plan, execute it as many times as you want with fftw_execute(plan), and clean up with fftw_destroy_plan(plan) (and fftw_free). FFTW provides two routines for creating plans for 2d and 3d transforms, and one routine for creating plans of arbitrary dimensionality. The 2d and 3d routines have the following signature: fftw_plan fftw_plan_dft_2d(int n0, int n1,                            fftw_complex *in, fftw_complex *out,                            int sign, unsigned flags); fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2,                            fftw_complex *in, fftw_complex *out,                            int sign, unsigned flags);

[링크 : https://www.fftw.org/fftw3_doc/Complex-Multi_002dDimensional-DFTs.html]

 

3.2.1 Row-major Format The multi-dimensional arrays passed to fftw_plan_dft etcetera are expected to be stored as a single contiguous block in row-major order (sometimes called “C order”). Basically, this means that as you step through adjacent memory locations, the first dimension’s index varies most slowly and the last dimension’s index varies most quickly. To be more explicit, let us consider an array of rank d whose dimensions are n0 × n1 × n2 × … × nd-1 . Now, we specify a location in the array by a sequence of d (zero-based) indices, one for each dimension: (i0, i1, i2,..., id-1). If the array is stored in row-major order, then this element is located at the position id-1 + nd-1 * (id-2 + nd-2 * (... + n1 * i0)). Note that, for the ordinary complex DFT, each element of the array must be of type fftw_complex; i.e. a (real, imaginary) pair of (double-precision) numbers. In the advanced FFTW interface, the physical dimensions n from which the indices are computed can be different from (larger than) the logical dimensions of the transform to be computed, in order to transform a subset of a larger array. Note also that, in the advanced interface, the expression above is multiplied by a stride to get the actual array index—this is useful in situations where each element of the multi-dimensional array is actually a data structure (or another array), and you just want to transform a single field. In the basic interface, however, the stride is 1.

[링크 : https://www.fftw.org/fftw3_doc/Row_002dmajor-Format.html]

 

A multi-dimensional array whose size is declared at compile time in C is already in row-major order. You don’t have to do anything special to transform it. For example: {      fftw_complex data[N0][N1][N2];      fftw_plan plan;      ...      plan = fftw_plan_dft_3d(N0, N1, N2, &data[0][0][0], &data[0][0][0],                              FFTW_FORWARD, FFTW_ESTIMATE);      ... } This will plan a 3d in-place transform of size N0 x N1 x N2. Notice how we took the address of the zero-th element to pass to the planner (we could also have used a typecast).

[링크 : https://www.fftw.org/fftw3_doc/Fixed_002dsize-Arrays-in-C.html]

[링크 : https://www.fftw.org/fftw3_doc/Multi_002ddimensional-Array-Format.html]

 

+

에라이 모르겠다 ㅠㅠ

switch (sz->rnk) {  case 1:       if (p->sign < 0) {    if (verbose > 2) printf("using plan_dft_r2c_1d\n");    return FFTW(plan_dft_r2c_1d)(sz->dims[0].n,  (bench_real *) p->in,  (bench_complex *) p->out, flags);       }       else {    if (verbose > 2) printf("using plan_dft_c2r_1d\n");    return FFTW(plan_dft_c2r_1d)(sz->dims[0].n,  (bench_complex *) p->in,  (bench_real *) p->out, flags);       }       break;  case 2:       if (p->sign < 0) {    if (verbose > 2) printf("using plan_dft_r2c_2d\n");    return FFTW(plan_dft_r2c_2d)(sz->dims[0].n, sz->dims[1].n, (bench_real *) p->in,  (bench_complex *) p->out, flags);       }       else {    if (verbose > 2) printf("using plan_dft_c2r_2d\n");    return FFTW(plan_dft_c2r_2d)(sz->dims[0].n, sz->dims[1].n, (bench_complex *) p->in,  (bench_real *) p->out, flags);       }       break;  case 3:       if (p->sign < 0) {    if (verbose > 2) printf("using plan_dft_r2c_3d\n");    return FFTW(plan_dft_r2c_3d)( sz->dims[0].n, sz->dims[1].n, sz->dims[2].n, (bench_real *) p->in, (bench_complex *) p->out, flags);       }       else {    if (verbose > 2) printf("using plan_dft_c2r_3d\n");    return FFTW(plan_dft_c2r_3d)( sz->dims[0].n, sz->dims[1].n, sz->dims[2].n, (bench_complex *) p->in, (bench_real *) p->out, flags);       }       break;  default: {       int *n = mkn(sz);       if (p->sign < 0) {    if (verbose > 2) printf("using plan_dft_r2c\n");    pln = FFTW(plan_dft_r2c)(sz->rnk, n,     (bench_real *) p->in,      (bench_complex *) p->out,     flags);       }       else {    if (verbose > 2) printf("using plan_dft_c2r\n");    pln = FFTW(plan_dft_c2r)(sz->rnk, n,     (bench_complex *) p->in,      (bench_real *) p->out,     flags);       }       bench_free(n);       return pln;  }      }

[링크 : https://github.com/FFTW/fftw3/blob/master/tests/bench.c]

 

+

mpi 예제도 있음.

그런데 3d로 한다고 하니 어마어마하게 데이터가 필요할 것 같은데

단순하게 3축 데이터라면 엄청 큰 공간에 딱 3줄 넣는 식이 될텐데 어떻게 해야 할까?

data = (FFTW_COMPLEX *) FFTW_MALLOC(sizeof(FFTW_COMPLEX) * NX*NY*NZ*sizeof(FFTW_COMPLEX));          /* create plan for forward DFT */     gettimeofday(&tv_start, NULL);     plan = FFTW_PLAN_DFT_3D(NX, NY, NZ, data, data,     FFTW_FORWARD, FFTW_ESTIMATE);     gettimeofday(&tv_stop, NULL);     deltaT =  tv_stop.tv_sec  - tv_start.tv_sec +       1e-6 * (tv_stop.tv_usec - tv_start.tv_usec);     printf("Plan init time (c2c)  : %f sec\n", deltaT);          /* initialize data */     for (i = 0; i < NX; ++i) {       for (j = 0; j < NY; ++j) {         for (k = 0; k < NZ; ++k) {           data[(i*NY + j)*NZ + k][0] = (i*NY + j)*NZ + k;           data[(i*NY + j)*NZ + k][1] = 0;         }       }     }

[링크 : https://github.com/pkestene/simpleFFTW/blob/master/c/fftw_test_3d_serial.c]

 

+

암만봐도 2d나 3d fft는 어떻게 써야 할지 감도 안와서 사용 시도 자체를 포기!