 A+
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 Fouriertransform 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 onedimension position function y(i)= 2*(delta)*(i1)*exp(((i1)*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 = (i1)*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 (RealtoReal Transforms, transform kinds, Realodd 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 realtoreal 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 "typeI" discrete sine transform (DSTI) 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 rgrid do i = 1, n r = i * delta ! (2a) y( i )= r * exp( r**2 ) end do call fftw_execute_r2r( my_plan, y, yy ) ! Loop over kgrid do i = 1, n ! Result of FFTW k = i * pi / ((n + 1) * delta) ! (2b) 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 gfortran8.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 (2a) and (2b)), 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 forfftw_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 ispi^(3/2) * exp(k^2 / 4)
(obtained from Wolfram Alpha or by noting that this is actually a 3D 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 counterintuitive, the result becomes a positive Gaussian.  I guess the rgrid 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 (hereexp(r^2)
).