Calcul de la FFT 3D et de la FFT inverse en C

Question

Je souhaite calculer les transformations FFT et arrière pour vérifier si elles fonctionnent de la même manière. J'ai une application de grande matrice 3D dans mon code j'ai essayé de le tester avec 4 * 4 * 4 matrice et voici mon code

`

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


int main() 
{ 
int N = 4; //Dimension of matrix 
unsigned int seed = 1; 
double *in = (double*)malloc(sizeof(double)*N*N*N); 
fftw_complex *out = fftw_malloc(sizeof(fftw_complex)*N*N*N); 
double *out1 = (double*)malloc(sizeof(double)*N*N*N); 

fftw_plan plan_backward; 
fftw_plan plan_forward; 

srand (seed); 

for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      in[i*(N*N) + j*N + k] = rand (); 
     } 
    } 
} 

printf(" Given matrix in\n"); 


for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      printf("%f\n", in[i*(N*N) + j*N + k]); 
     } 
    } 
} 


printf("\n"); 

plan_backward = fftw_plan_dft_r2c_3d (N, N, N, in, out, FFTW_ESTIMATE); 

fftw_execute (plan_backward); 

fftw_destroy_plan (plan_backward); 

printf("out matrix\n"); 

for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      printf("%f\t%f\n", creal(out[i*(N*N) + j*N + k]), cimag(out[i*(N*N) + j*N + k])); 
     } 
    } 
} 

printf("\n"); 

plan_forward = fftw_plan_dft_c2r_3d (N, N, N, out, out1, FFTW_ESTIMATE); 

fftw_execute (plan_forward); 

fftw_destroy_plan (plan_forward); 

printf("out1 matrix\n"); 


for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      printf("%f\n", out1[i*(N*N) + j*N + k]); 
     } 
    } 
} 

fftw_free(in); 
free(out); 
fftw_free(out1); 

return 0; 

}` 

Apparemment, mes résultats transformés ne sont pas les mêmes. Je ne comprends pas ce qui ne va pas?

Answer 1

Votre FFT n'est pas normalisée. Il y a un facteur constant entre votre entrée et votre sortie.

Jetez un oeil here

Ces transformations sont non normalisée, donc un R2C suivi d'un C2R Transform (ou vice versa) se traduira par les données d'origine mis à l'échelle par le nombre d'éléments qui sont des données réelles, le produit des dimensions (logiques) des données réelles.

Le facteur doit donc être N * N * N. Juste divisez vos données par ce facteur et vous devriez récupérer les mêmes données que votre entrée.

Answer 2

Comme décrit dans le manuel de FFTW3:

Ces transformées sont non normalisés, de sorte qu'un R2C suivie d'une C2R transformée (ou vice versa) se traduira par les données d'origine mises à l'échelle par le nombre d'éléments de données réelles C'est-à-dire le produit des dimensions (logiques) des données réelles.

Le nombre de dimensions réelles est 4^3 dans votre cas et en divisant le premier résultat 115474520512 par ce numéro rend la première entrée 115474520512/(4^3) = 1804289383.