Evaluating the fast Fourier transform of Gaussian function in FORTRAN using FFTW3 library

  • A+
Category:Languages

I am trying to write a FORTRAN code to evaluate the fast Fourier transform of the Gaussian function f(r)=exp(-(r^2)) using FFTW3 library. As everyone knows, the Fourier transform of the Gaussian function is another Gaussian function.

I consider evaluating the Fourier-transform integral of the Gaussian function in the spherical coordinate.

Hence the resulting integral can be simplified to be integral of [r*exp(-(r^2))*sin(kr)]dr.

I wrote the following FORTRAN code to evaluate the discrete SINE transform DST which is the discrete Fourier transform DFT using a PURELY real input array. DST is performed by C_FFTW_RODFT00 existing in FFTW3, taking into account that the discrete values in position space are r=i*delta (i=1,2,...,1024), and the input array for DST is the function r*exp(-(r^2)) NOT the Gaussian. The sine function in the integral of [r*exp(-(r^2))*sin(kr)]dr resulting from the INTEGRATION over the SPHERICAL coordinates, and it is NOT the imaginary part of exp(ik.r) that appears when taking the analytic Fourier transform in general.

However, the result is not a Gaussian function in the momentum space.

Module FFTW3  use, intrinsic :: iso_c_binding include 'fftw3.f03' end module    program sine_FFT_transform use FFTW3 implicit none integer, parameter :: dp=selected_real_kind(8)  real(kind=dp), parameter :: pi=acos(-1.0_dp) integer, parameter :: n=1024  real(kind=dp) :: delta, k real(kind=dp) :: numerical_F_transform integer :: i type(C_PTR) ::  my_plan real(C_DOUBLE), dimension(1024) :: y real(C_DOUBLE), dimension(1024) :: yy, yk integer(C_FFTW_R2R_KIND) :: C_FFTW_RODFT00  my_plan= fftw_plan_r2r_1d(1024,y,yy,FFTW_FORWARD, FFTW_ESTIMATE)  delta=0.0125_dp do i=1, n        !inserting the input one-dimension position function y(i)= 2*(delta)*(i-1)*exp(-((i-1)*delta)**2)  ! I multiplied by 2 due to the definition of C_FFTW_RODFT00 in FFTW3 end do  call fftw_execute_r2r(my_plan, y,yy)    do i=2, n k = (i-1)*pi/n/delta  yk(i) = 4*pi*delta*yy(i)/2  !I divide by 2 due to the definition of                              !C_FFTW_RODFT00 numerical_F_transform=yk(i)/k write(11,*) i,k,numerical_F_transform end do call fftw_destroy_plan(my_plan)  end program  

Executing the previous code gives the following plot which is not for Gaussian function. Evaluating the fast Fourier transform of Gaussian function in FORTRAN using FFTW3 library Can anyone help me understand what the problem is? I guess the problem is mainly due to FFTW3. Maybe I did not use it properly especially concerning the boundary conditions.

 


Looking at the related pages in the FFTW site (Real-to-Real Transforms, transform kinds, Real-odd DFT (DST)) and the header file for Fortran, it seems that FFTW expects FFTW_RODFT00 etc rather than FFTW_FORWARD for specifying the kind of real-to-real transform. For example,

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE ) my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE ) 

performs the "type-I" discrete sine transform (DST-I) shown in the above page. This modification seems to fix the problem (i.e., makes the Fourier transform a Gaussian with positive values).


The following is a slightly modified version of OP's code to experiment the above modification:

! ... only the modified part is shown... real(dp) :: delta, k, r, fftw, num, ana integer :: i, j, n type(C_PTR) ::  my_plan real(C_DOUBLE), allocatable :: y(:), yy(:)  delta = 0.0125_dp ; n = 1024   ! rmax = 12.8 ! delta = 0.1_dp    ; n = 128    ! rmax = 12.8 ! delta = 0.2_dp    ; n = 64    ! rmax = 12.8 ! delta = 0.4_dp    ; n = 32    ! rmax = 12.8  allocate( y( n ), yy( n ) )  ! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE ) my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )  ! Loop over r-grid do i = 1, n     r = i * delta              ! (2-a)     y( i )= r * exp( -r**2 ) end do  call fftw_execute_r2r( my_plan, y, yy )  ! Loop over k-grid do i = 1, n      ! Result of FFTW     k = i * pi / ((n + 1) * delta)    ! (2-b)     fftw = 4 * pi * delta * yy( i ) / k / 2   ! the last 2 due to RODFT00      ! Numerical result via quadrature     num = 0     do j = 1, n         r = j * delta         num = num + r * exp( -r**2 ) * sin( k * r )     enddo     num = num * 4 * pi * delta / k      ! Analytical result     ana = sqrt( pi )**3 * exp( -k**2 / 4 )      ! Output     write(10,*) k, fftw     write(20,*) k, num     write(30,*) k, ana end do 

Compile (with gfortran-8.2 + FFTW3.3.8 + OSX10.11):

$ gfortran -fcheck=all -Wall sine.f90 -I/usr/local/Cellar/fftw/3.3.8/include -L/usr/local/Cellar/fftw/3.3.8/lib -lfftw3 

If we use FFTW_FORWARD as in the original code, we get

Evaluating the fast Fourier transform of Gaussian function in FORTRAN using FFTW3 library

which has a negative lobe (where fort.10, fort.20, and fort.30 correspond to FFTW, quadrature, and analytical results). Modifying the code to use FFTW_RODFT00 changes the result as below, so the modification seems to be working (but please see below for the grid definition).

Evaluating the fast Fourier transform of Gaussian function in FORTRAN using FFTW3 library


Additional notes

  • I have slightly modified the grid definition for r and k in my code (Lines (2-a) and (2-b)), which is found to improve the accuracy. But I'm still not sure whether the above definition matches the definition used by FFTW, so please read the manual for details...
  • The fftw3.f03 header file gives the interface for fftw_plan_r2r_1d

    type(C_PTR) function fftw_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftw_plan_r2r_1d')   import   integer(C_INT), value :: n   real(C_DOUBLE), dimension(*), intent(out) :: in   real(C_DOUBLE), dimension(*), intent(out) :: out   integer(C_FFTW_R2R_KIND), value :: kind   integer(C_INT), value :: flags end function fftw_plan_r2r_1d 
  • (Because of no Tex support, this part is very ugly...) The integral of 4 pi r^2 * exp(-r^2) * sin(kr)/(kr) for r = 0 -> infinite is pi^(3/2) * exp(-k^2 / 4) (obtained from Wolfram Alpha or by noting that this is actually a 3-D Fourier transform of exp(-(x^2 + y^2 + z^2)) by exp(-i*(k1 x + k2 y + k3 z)) with k =(k1,k2,k3)). So, although a bit counter-intuitive, the result becomes a positive Gaussian.

  • I guess the r-grid can be chosen much coarser (e.g. delta up to 0.4), which gives almost the same accuracy as long as it covers the frequency domain of the transformed function (here exp(-r^2)).

Comment

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: