fortran complex to real fftw issue

Question

Je travaille actuellement sur un projet où j'ai besoin de mettre en œuvre une transformée de Fourier et une transformation inverse. Je testais un programme que j'ai modifié à partir d'un exemple en ligne; l'impression ou d'écriture sont normalement à des fins de débogage:

program testit 
    INCLUDE 'fftw3.f' 
    double complex out!, in 
    real in 
    parameter (N=100) 
    dimension in(N), out(N) 
    integer*8 p,p2 
    integer i,j 
    real x 
    real fact 
    write(*,*)"stuff in data" 
    OPEN(UNIT=12, FILE="input.txt", ACTION="write", STATUS="replace") 
    OPEN(UNIT=20, FILE="dftoutput.txt", ACTION="write", STATUS="replace") 
    x=0 
    in = 0 
    do i=1,N/2 
    in(i)=1 
    enddo 
    do i=1,N 
    write(*,"(f10.2,1x,f10.2)")in(i) 
    WRITE(12,*)real(in(i)) 
    enddo 
    write(*,*)"create plans" 
    call dfftw_plan_dft_r2c_1d(p ,N,in,out,FFTW_ESTIMATE) 
    call dfftw_plan_dft_c2r_1d(p2,N,in,out,FFTW_ESTIMATE) 
    write(*,*)"do it" 
    call dfftw_execute_dft_r2c(p,in,out) 
    do i=1,N 
    write(*,"(f12.4,1x,f12.4)")out(i) 
    WRITE(20,*)abs(out(i)) 
    enddo 
    write(*,*)"undo it" 
    call dfftw_execute_dft_c2r(p2,in,out) 
    fact=1.0/N 
    do i=1,N 
    write(*,)in(i) 
    write(*,)out(i) 
    enddo 
    write(*,*)"clean up" 
    call dfftw_destroy_plan(p,in,out) 
    call dfftw_destroy_plan(p2,in,out) 
    end program 

Le vrai à la transformation complexe fonctionne très bien. La transformation inverse donne de fausses valeurs et modifie les variables d'entrée et de sortie. Je ne sais pas quel est le problème et je n'ai pas pu trouver de réponses en ligne. L'aide est appréciée.

Merci d'avance!

Tchad W. Freer

Edit: Je me demandais aussi s'il y a une fonction similaire à fftshift() et ifftshift() de Matlab dans le paquet de FFTW.

Answer 1

Il y a un problème de précision: sur la plupart des ordinateurs, real fait référence à des nombres à virgule flottante de simple précision de 32 bits. real*8 peut être utilisé pour spécifier des nombres à virgule flottante en double précision, pour être cohérent avec double complex et dfftw_....

Pour une transformation FFTW réelle à complexe, si la taille de l'entrée est N real*8, le size of the output est N/2+1 double complex. Puisque l'entrée est réelle, coefficients of negative frequencies (higher than N/2+1) are conjugate of positive frequencies, et leur calcul est évité.

Selon la documentation de FFTW concernant planner flags,

FFTW_PRESERVE_INPUT spécifie qu'une transformation de hors place ne doit pas modifier son tableau d'entrée. Ceci est habituellement la valeur par défaut, sauf pour C2R et hc2r (c.-à-complexe à vrai) pour lequel FFTW_DESTROY_INPUT transforme la valeur par défaut est ...

Si vous souhaitez enregistrer les coefficients dans l'espace de Fourier out, soit copie le tableau, ou essayez d'utiliser FFTW_PRESERVE_INPUT. La première solution semble la meilleure solution si vous avez assez de mémoire.

L'ordre des arguments pour les fonctions FFTW est toujours origin,destination. Comme l'arrière de transformation est réalisée out-in:

call dfftw_plan_dft_c2r_1d(p2,N,out,in,FFTW_ESTIMATE) 
... 
call dfftw_execute_dft_c2r(p2,out,in) 

Voici un code basé sur le vôtre, qui effectue en avant et en arrière fft. Il est compilé par gfortran main.f90 -o main -lfftw3 -lm -Wall:

program testit 
INCLUDE 'fftw3.f' 
double complex out!, in 
real*8 in 
parameter (N=100) 
dimension in(N), out((N/2+1)) 
integer*8 p,p2 
integer i 
real x 
real fact 
write(*,*)"stuff in data" 
OPEN(UNIT=12, FILE="input.txt", ACTION="write", STATUS="replace") 
OPEN(UNIT=20, FILE="dftoutput.txt", ACTION="write", STATUS="replace") 
x=0 
in = 0 
do i=1,N/2 
    in(i)=1 
enddo 
do i=1,N 
    write(*,"(f10.2,1x/)")in(i) 
    WRITE(12,*)real(in(i)) 
enddo 
write(*,*)"create plans" 
call dfftw_plan_dft_r2c_1d(p ,N,in,out,FFTW_ESTIMATE) 
call dfftw_plan_dft_c2r_1d(p2,N,out,in,FFTW_ESTIMATE) 
write(*,*)"do it" 
call dfftw_execute_dft_r2c(p,in,out) 
do i=1,N/2+1 
    write(*,"(f12.4,1x,f12.4/)")out(i) 
    WRITE(20,*)abs(out(i)) 
enddo 
write(*,*)"undo it" 
call dfftw_execute_dft_c2r(p2,out,in) 
fact=1.0/N 
write(*,*)"input back, except for scaling" 
do i=1,N 
    write(*,"(f10.2,1x)")in(i) 
enddo 
write(*,*)"output may be modified by c2r...except if FFTW_PRESERVE_INPUT is set" 
do i=1,N/2+1 
    write(*,"(f10.2,1x,f10.2/)")out(i) 
enddo 
write(*,*)"clean up" 
call dfftw_destroy_plan(p,in,out) 
call dfftw_destroy_plan(p2,out,in) 
end program 

Après la fft avant in ->out et la fft arrière out ->in, in est identique à sa valeur initiale, à l'exception d'un facteur d'échelle N.