I am trying to write a FORTRAN code to evaluate the fast Fourier transform of the Gaussian function
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
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).
- 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...
fftw3.f03header file gives the interface for
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.
deltaup to 0.4), which gives almost the same accuracy as long as it covers the frequency domain of the transformed function (here