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

Question

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.

Showery

Hi, the plot of your function is missing which makes it difficult to answer, as does being sat on a bus. But the 2 slightly strange things I can see are the one roygvib mentions, that you create the plan with yk but use yy in the transform, and that for the frequencies while your version is not wrong it is more usual to have the second half of the range as negative frequencies - that way if you have done the FFT right the plot will look lie a gaussian, as you have done it it will not (it will be 2 "1/2 Gaussians", one at either end of the range) – Anglim

Thank roygvib and Ian Bush. Your answers were helpful. there is no problem with the negative frequencies, they can be generated once the positive frequencies are correct. After editing the code, there is still negative transforms in the graph which is below the x axis, and this does not make the plot to be Gaussian. The problem still exists !!!! – Showery

Answer 1

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

Ticonderoga

This is the correct answer. I really did not have time to dig into this, but now seeing your change, it is indeed necessary. – Selfridge

Thanks a lot @roygvib. I am really grateful for you. Your answer is correct and perfect. – Showery

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

Roethke