Erreurs avec des appels FFTW répétés

Question

J'ai un problème étrange que je ne peux pas résoudre. J'ai fait cela comme un exemple simple qui démontre le problème. J'ai une onde sinusoïdale définie entre [0, 2 * pi]. Je prends la transformée de Fourier en utilisant FFTW. Ensuite, j'ai une boucle for où je prends à plusieurs reprises la transformée de Fourier inverse. Dans chaque itération, je prends la moyenne de ma solution et j'imprime les résultats. Je m'attends à ce que la moyenne reste la même à chaque itération car il n'y a pas de changement de solution, y. Cependant, quand je prends N = 256 et d'autres valeurs paires de N, je note que la moyenne croît comme s'il y avait des erreurs numériques. Cependant, si je choisis, disons, N = 255 ou N = 257, ce n'est pas le cas et j'obtiens ce qui est attendu (avg = 0.0 pour chaque itération).

code:

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

int main(void) 
{ 
    int N = 256; 
    double dx = 2.0 * M_PI/(double)N, dt = 1.0e-3; 
    double *x, *y; 

    x = (double *) malloc (sizeof (double) * N); 
    y = (double *) malloc (sizeof (double) * N); 

    // initial conditions 
    for (int i = 0; i < N; i++) { 
    x[i] = (double)i * dx; 
    y[i] = sin(x[i]); 
    } 

    fftw_complex yhat[N/2 + 1]; 
    fftw_plan fftwplan, fftwplan2;  

    // forward plan 
    fftwplan = fftw_plan_dft_r2c_1d(N, y, yhat, FFTW_ESTIMATE); 
    fftw_execute(fftwplan); 

    // set N/2th mode to zero if N is even 
    if (N % 2 < 1.0e-13) { 
    yhat[N/2][0] = 0.0; 
    yhat[N/2][1] = 0.0; 
    } 

    // backward plan 
    fftwplan2 = fftw_plan_dft_c2r_1d(N, yhat, y, FFTW_ESTIMATE); 

    for (int i = 0; i < 50; i++) { 
    // yhat to y 
    fftw_execute(fftwplan2); 

    // rescale 
    for (int j = 0; j < N; j++) { 
     y[j] = y[j]/(double)N; 
    } 

    double avg = 0.0; 
    for (int j = 0; j < N; j++) { 
     avg += y[j]; 
    } 
    printf("%.15f\n", avg/N); 
    } 

    fftw_destroy_plan(fftwplan); 
    fftw_destroy_plan(fftwplan2); 
    void fftw_cleanup(void); 
    free(x); 
    free(y); 
    return 0; 
} 

sortie N = 256:

0.000000000000000 
0.000000000000000 
0.000000000000000 
-0.000000000000000 
0.000000000000000 
0.000000000000022 
-0.000000000000007 
-0.000000000000039 
0.000000000000161 
-0.000000000000314 
0.000000000000369 
0.000000000004775 
-0.000000000007390 
-0.000000000079126 
-0.000000000009457 
-0.000000000462023 
0.000000000900855 
-0.000000000196451 
0.000000000931323 
-0.000000009895302 
0.000000039348379 
0.000000133179128 
0.000000260770321 
-0.000003233551979 
0.000008285045624 
-0.000016331672668 
0.000067450106144 
-0.000166893005371 
0.001059055328369 
-0.002521514892578 
0.005493164062500 
-0.029907226562500 
0.093383789062500 
-0.339111328125000 
1.208251953125000 
-3.937500000000000 
13.654296875000000 
-43.812500000000000 
161.109375000000000 
-479.250000000000000 
1785.500000000000000 
-5369.000000000000000 
19376.000000000000000 
-66372.000000000000000 
221104.000000000000000 
-753792.000000000000000 
2387712.000000000000000 
-8603776.000000000000000 
29706240.000000000000000 
-96833536.000000000000000 

Toutes les idées?

Answer 1

libfftw a l'habitude odieuse de modifier ses entrées. Sauvegardez yhat si vous voulez faire des transformations inverses répétées.

OTOH, c'est pervers, mais pourquoi répétez-vous la même opération si vous ne vous attendez pas à des résultats différents? (Bien que ce soit le cas)

Answer 2

Comme indiqué dans les commentaires: «Si vous souhaitez conserver les données d'entrée inchangé, utilisez l'indicateur FFTW_PRESERVE_INPUT par http://www.fftw.org/doc/Planner-Flags.html »

Par exemple:

// backward plan 
fftwplan2 = fftw_plan_dft_c2r_1d(N, yhat, y, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);