Qfyilyi

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.

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

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).

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)).

Of course there are negative components of the real part to the FFT of a limited Gaussian spectrum. You are just using the real part of the transform. So your plot is absolutely correct.

You seem to be mistaking the real part with the magnitude, which of course would not be negative. For that you would need to fftw_plan_dft_r2c_1d and then calculate the absolute values of the complex coefficients. Or you might be mistaking the Fourier transform with a limited DFT.

You might want to check here to convince yourself of the correctness of you calculation above:

http://docs.mantidproject.org/nightly/algorithms/FFT-v1.html

Please do keep in mind that the plots on the above page are shifted, so that the 0 frequency is in the middle of the spectrum.

Citing yourself, the nummeric integration of [r*exp(-(r^2))*sin(kr)]dr would have negative components for all k>1 if normalised to 0 for highest frequency.

TLDR: Your plot is absolute state of the art and inline with discrete and limited functional analysis.