Путаница в интерфейсе гуру FFTW3: 3 одновременных сложных БПФ - PullRequest
Новая
       15

Путаница в интерфейсе гуру FFTW3: 3 одновременных сложных БПФ

1 голос
/ 23 октября 2019

Я пытаюсь использовать интерфейс гуру FFTW3 . Тем не менее, описание в документации вызывает у меня больше путаницы, чем ясность. Я хочу создать план для 3 скалярных комплексных БПФ, каждое из которых имеет длину n, которые должны быть выровнены в памяти один за другим (чередующиеся массивы). Я не понимаю, подразумевает ли это rank == 3 или howmany_rank == 3. Я также не знаю, какие значения должны иметь массивы fftw_iodim *dims и const fftw_iodim *howmany_dims.

Хотя этот стековый обменный пост , кажется, имеет хороший ответ, он использует ту же терминологиюдокументации и не «создает образ в моей голове» того, как FFT будут выровнены в памяти для разных rank с и разных howmany_rank с. То есть ответ не помогает с тем, что я ищу.

Для справки, это определение функции распределения плана: 

fftw_plan fftw_plan_guru_dft(
      int rank, const fftw_iodim *dims,
      int howmany_rank, const fftw_iodim *howmany_dims,
      fftw_complex *in, fftw_complex *out,
      int sign, unsigned flags);

Обратите внимание, что используется «стандартный интерфейс». "невозможно, так как позже я должен переключиться на интерфейс guru64 (64 бит).

1 Ответ

1 голос
/ 25 октября 2019

fftw_plan fftw_plan_guru_dft() позволяет определить 1D-2D-3D-4D ... DFT (rank = 1,2,3,4 ...), примененные несколько раз с использованием начальных точек, расположенных на 1D, 2D, 3D, 4D ... grid (howmany_rank = 1,2,3,4 ....

Аргумент rank указывает число измерений, к которым применяется DFT. Расширения и компоновкаэти размеры задаются аргументом *dims. После выполнения преобразования DFT индексы этих измерений соответствуют частотам. Эти аргументы определяют многомерный DFT, который будет применяться.

Аргумент howmany_rank и howmany_dims определяет начальные точки для многократного применения многомерного преобразования, определенного выше. Начальные точки могут быть расположены на многомерном массиве, описываемом extends и step. Аргумент howmany_rank описатьчисло измерений массива начальных точек.

Рассмотрим поле, содержащее 3 скалярных компонента u, v, w, отобранные на прямой в равномерно расположенной точке абсциссы x_i = iL / N. Полехранится в памяти как чередующийся непрерывный двумерный массив измерения (N, 3): 

u(x_0) v(x_0) w(x_0) u(x_1) v(x_1) w(x_1) u(x_2) v(x_2) w(x_2)... u(x_N-1) v(x_N-1) w(x_N-1)

1D DFT должно выполняться по пространственной координате x. Следовательно, rank равно 1. Каждый 1D DFT имеет длину N, а расстояние (или шаг) между двумя последовательными элементами (u_(x_0) и u(x_1)) равно 3. Следовательно, dims[0].n=N, dims[0].is=3 и dims[0].os=3. is для входного шага и os для выходного шага.

1D DFT должен выполняться несколько раз, один раз для каждого компонента. Поскольку начальные точки этих ДПФ u(x_0), v(x_0) и w(x_0) расположены на регулярном расстоянии, позиции этих начальных точек определяют массив измерения 1. Следовательно, howmany_rank=1. Более того, поскольку есть 3 смежные начальные точки, расположение массива начальных точек определяется как howmany_dims[0].n=3, howmany_dims[0].is=1 и howmany_dims[0].os=1.

Пример кода, приведенный в моем ответе на Как использовать интерфейс fftw Guru (извините, что он вам не помог!) Может быть легко адаптирован как: 

#include<stdlib.h>
#include<complex.h>
#include<math.h>
#include<fftw3.h>

int main(void){

    fftw_plan p;
    unsigned long int N = 10;
    unsigned long int DOF = 3;
    fftw_complex *in=fftw_malloc(N*DOF*sizeof(fftw_complex));
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    fftw_complex *out=fftw_malloc(N*DOF*sizeof(fftw_complex));
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    unsigned int i;

    int rank=1;
    fftw_iodim64 *dims=malloc(rank*sizeof(fftw_iodim64));
    if(dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    dims[0].n=N;
    dims[0].is=DOF;
    dims[0].os=DOF;

    int howmany_rank=1;
    fftw_iodim64 *howmany_dims=malloc(howmany_rank*sizeof(fftw_iodim64));
    if(howmany_dims==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    howmany_dims[0].n=DOF;
    howmany_dims[0].is=1;
    howmany_dims[0].os=1;

    printf("sizeof fftw complex %ld\n",sizeof(fftw_complex));
    printf("sizeof fftw_iodim64 %ld\n",sizeof(fftw_iodim64));
    printf("creating the plan\n");
    p=fftw_plan_guru64_dft(rank, dims,howmany_rank, howmany_dims,in, out,FFTW_FORWARD, FFTW_ESTIMATE);
    if (p==NULL){fprintf(stderr,"plan creation failed\n");exit(1);}
    printf("created the plan\n");

    for(i=0;i<N;i++){
                in[i*DOF]=30.+12.*sin(2*3.1415926535*i/((double)N));
                in[i*DOF+1]=42.0;
                in[i*DOF+2]=1.0;
    }

    fftw_execute(p);

    for (i = 0; i < N; i++){
        printf("result: %d || %g %gI | %g %gI | %g %gI\n", i, creal(out[i*DOF]), cimag(out[i*DOF]),creal(out[i*DOF+1]), cimag(out[i*DOF+1]),creal(out[i*DOF+2]), cimag(out[i*DOF+2]));
    }


    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

    free(dims);
    free(howmany_dims);

    return(0);
}

И скомпилировано gcc main.c -o main -lfftw3 -lm -Wall.

...