# 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. 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). • 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)`).