2016-10-08 6 views
3

다차원 푸리에 변환에 fftw_plan_dft을 사용했습니다.fftw 전문가 인터페이스 사용 방법

fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in, 
         fftw_complex *out, int sign, unsigned flags); 

지금 내가 FFTW 전문가 인터페이스를 사용해야 할 것 같습니다, FFTW 64 비트 정수를 전달하려는.

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

하지만 howmany_rankhowmany_dims의 의미가 무엇인지 이해가 안 돼요. fftw_plan_guru_dft의 매뉴얼 말한다

이 두 기능

복잡한 데이터, 인터리브 분할 형식 다차원 DFT 각각 계획. 변환 차원은 차원 (howmany_rank, howmany_dims)의 다차원 벡터 (반복)에 대해 (순위, 희미하게) 표시됩니다. dims 및 howmany_dims는 각각 길이 순위와 howmany_rank의 fftw_iodim 배열을 가리켜 야합니다.

"다차원 벡터 (반복) 치수 (howmany_rank, howmany_dims)"는 무엇을 의미하는지 알고 있습니다. 예제를 주거나이 전문가 인터페이스를 사용하는 방법을 설명해 주시겠습니까?

+0

이에 pyfftw 태그 때문에, https://github.com/pyFFTW/pyFFTW/blob/ ([다음은 예입니다] 마스터/pyfftw/pyfftw.pyx # L1104). 이 줄은 전문가 인터페이스 함수를 가리키는 함수 포인터를 호출합니다. 위의 python을 읽고 매개 변수를 채우는 방법을 알아보십시오. –

+0

감사합니다. @Henry Gomersall. 이 라인이 전문가 인터페이스를 사용하고있는 것을 볼 수 있습니다. 이 python 코드를 읽기 위해 열심히 노력하고 있습니다. 그러나, 그것은 조금 오랫동안 점등되었습니다, 당신이'howmany_rank'와'howmany_dims'의 의미를 나에게 설명해도 괜찮습니까? –

+0

문서를 읽고 그것에 대해 열심히 생각해 보시기 바랍니다. 나는 그 짧은 지름길을 볼 수 없다. 당신은 임의의 변환을 지정하기 위해 다양한 정보 비트가 필요하다는 것을 알게 될 것이며 그 시점에서 분명해질 것입니다. 또는 링크 된 코드를 이해하는 데 약간의 노력을 기울여 모든 것을 명확하게 할 수 있습니다 (모양, 보폭 및 축을 FFTW 매개 변수로 변환하는 일반적인 매핑을 설명하기 때문에 좋은 예입니다). –

답변

1

다차원 배열의 크기와 스트라이드가 2^32보다 큰 경우 64 bit guru interface이 유용합니다.

복잡한 DTFs 복잡한를 만드는 함수의 원형은 다음

  • rank 치수의 FFTW 수행되는 변형, 즉 수의 등급이다

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

    .

  • dims은 크기가 rank 인 배열입니다. 각 차원 i의 경우 dims[i].n은 행의 크기이고 dims[i].is은 입력 배열의 행 사이의 보폭이고 dims[i].os은 출력 배열의 행 사이의 보폭입니다. 예를 들어 배열이 연속적으로 메모리에 있으면 the documentation of the guru interface은 반복을 사용할 것을 제안합니다 dims[i].is = n[i+1] * dims[i+1].is. 수행 할 변환 수와 시작점 사이의 간격은 howmany_rankhowmany_dims입니다.
  • howmany_rank은 특정 오프셋을 특징으로하는 변환 수를 지정합니다.
  • howmany_dims은 크기가 howmany_rank 인 배열입니다. i의 각 변환에 대해 howmany_dims[i].n은 계산할 변환 수이며, 각 입력은 howmany_dims[i].is 사이의 오프셋과 출력 사이의 오프셋은 howmany_dims[i].os입니다. 이 fftw_plan_dft_3d()과 같은 일을 수행하도록

다음 코드는 fftw_plan_guru64_dft()를 호출합니다. gcc main.c -o main -lfftw3 -lm -Wall에 의해 컴파일 될 수있다 : 예를 들어

#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 M = 12; 
    unsigned long int P = 14; 
    fftw_complex *in=fftw_malloc(N*M*P*sizeof(fftw_complex)); 
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    fftw_complex *out=fftw_malloc(N*M*P*sizeof(fftw_complex)); 
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    unsigned int i,j,k; 

    int rank=3; 
    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=P*M; 
    dims[0].os=P*M; 
    dims[1].n=M; 
    dims[1].is=P; 
    dims[1].os=P; 
    dims[2].n=P; 
    dims[2].is=1; 
    dims[2].os=1; 

    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=1; 
    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++){ 
     for(j=0;j<M;j++){ 
      for(k=0;k<P;k++){ 
       //printf("ijk\n"); 
       in[(i*M+j)*P+k]=30.+12.*sin(2*3.1415926535*i/((double)N))*sin(2*3.1415926535*j/((double)M))*sin(2*3.1415926535*k/((double)P))*I; 
      } 
     } 
    } 

    fftw_execute(p); 

    for (i = 0; i < N; i++){ 
     for (j = 0; j < M; j++){ 
      for (k = 0; k < P; k++){ 
       printf("result: %d %d %d %g %gI\n", i,j,k, creal(out[(i*M+j)*P+k]), cimag(out[(i*M+j)*P+k])); 
      } 
     } 
    } 


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

    free(dims); 
    free(howmany_dims); 

    return(0); 
} 

, 도사 인터페이스는 복잡한 3D 전계의 DFT를 계산하기 위해 사용될 수있다. 그리드의 각 지점에서 전기장은 크기가 3 인 벡터입니다. 따라서 전계를 4D 배열로 저장할 수 있습니다. 마지막 차원은 벡터의 구성 요소를 지정합니다.마지막으로, 전문가의 인터페이스는 한 번에 세 가지 차원 DFT들을 수행하는 데 사용할 수 있습니다 :

#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 M = 12; 
    unsigned long int P = 14; 
    unsigned long int DOF = 3; 
    fftw_complex *in=fftw_malloc(N*M*P*DOF*sizeof(fftw_complex)); 
    if(in==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    fftw_complex *out=fftw_malloc(N*M*P*DOF*sizeof(fftw_complex)); 
    if(out==NULL){fprintf(stderr,"malloc failed\n");exit(1);} 
    unsigned int i,j,k; 

    int rank=3; 
    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=P*M*DOF; 
    dims[0].os=P*M*DOF; 
    dims[1].n=M; 
    dims[1].is=P*DOF; 
    dims[1].os=P*DOF; 
    dims[2].n=P; 
    dims[2].is=DOF; 
    dims[2].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=3; 
    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++){ 
     for(j=0;j<M;j++){ 
      for(k=0;k<P;k++){ 
       //printf("ijk\n"); 
       in[((i*M+j)*P+k)*3]=30.+12.*sin(2*3.1415926535*i/((double)N))*sin(2*3.1415926535*j/((double)M))*sin(2*3.1415926535*k/((double)P))*I; 
       in[((i*M+j)*P+k)*3+1]=42.0; 
       in[((i*M+j)*P+k)*3+2]=1.0; 
      } 
     } 
    } 

    fftw_execute(p); 

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


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

    free(dims); 
    free(howmany_dims); 

    return(0); 
}