0

I'm trying to get the 3D Fourier Transform of the gaussian function e^(-r^(2)/2) in python using the numpy.fft library. I've attempted using different ffts from the library with different inputs, shifting the results with np.fft.fftshift, trying to find a multiplicative factor and many other things, the last thing I tried was using the 1D fft function, and then cubing the result, here's the corresponding source code:

import numpy as np

R = float(10)
N = float(100)

y= np.dtype(np.float64)

dr = R/N

def F(x):
    return np.exp(-((x*dr)**2)/2)

Frange = np.arange(1,int(N)+1)
y = np.zeros((int(N)))

i = 0
while i<int(N):

    y[i] = F(Frange[i])
    i += 1

y = y/3

y_fft = np.fft.fftshift(np.abs(np.fft.fft(y)))**3

print (y_fft)

The first values I get:

4.62e-03, 4.63e-03, 4.65e-03, 4.69e-03, 4.74e-03

According to Lado, Fred. (1971) Numerical Fourier transforms in one, two, and three dimensions for liquid state calculations, the analytic solution to the problem is: (2pi )^(3/2)*e^(-k^(2)/2)

And the first values of the analytic solution with the same values of R and N are:

14.99, 12.92, 10.10, 7.15, 4.58

I also created a DFT program using a formula provided in the previous article which gives the expected results, but I haven't been able to replicate the analytic results in any of my attempts using the NumPy or SciPy fft libraries. Here's my program for the analytic and DFT results:

import math
import numpy as np

def F(r):   
    x=math.exp((-1/2)*(r**2))
    return x

def FT(r):
    x=((2*math.pi)**(3/2))*(math.exp((-1/2)*(r**2)))
    return x

R = float(10)
N = int(100)

ft = np.zeros(N)
fta = np.zeros(N)
dr = R/N
dk = math.pi/R

print ("\tk \t\t\t Discrete \t\t\t Analytic")

for j in range (1, N):
    kj = j*dk
    
#Discrete Transform
    sum = 0  
    for i in range(1, N):
        ri = i*dr
        sum = sum + (dr*ri*(F(ri))*(math.sin(kj*ri)))
        
    ft[j] = ((4*math.pi)/kj)*sum

#Analytic Transform
    fta[j] = FT(kj)
    
    #Print results
    print(kj, f" \t\t{ft[j]:.10E} \t\t{fta[j]:.10E}")

And these are the first few results:

 k                             Discrete                        Analytic
0.3141592653589793              1.4991263193E+01                1.4991263193E+01
0.6283185307179586              1.2928362116E+01                1.2928362116E+01
0.9424777960769379              1.0101494686E+01                1.0101494686E+01
1.2566370614359172              7.1509645344E+00                7.1509645344E+00
1.5707963267948966              4.5864901093E+00                4.5864901093E+00
Pedro Gutiérrez
  • 1
  • 1
  • ran your code but `print (y_fft)` can't reproduce the `1.74e+03, 1.41e+03, 7.86e+02, 3.27e+02, 1.22e+02` – manaclan Nov 29 '21 at 06:03
  • Thanks, you're right, I confused the results with those of another test I had previously run. I edited my question. – Pedro Gutiérrez Nov 29 '21 at 06:24
  • also, where is R and N in your gaussian function ? Is `14.99, 12.92, 10.10, 7.15, 4.58` the result when I replace k with `1,2,3,4,5` in the function `(2pi )^(3/2)*e^(-k^(2)/2)` – manaclan Nov 29 '21 at 06:38
  • `14.99, 12.92, 10.10, 7.15, 4.58` is the result when you replace k with `1*dk, 2*dk, 3*dk, ...` and so on with `dk = pi/R` in the analytic solution, I left my source code for the analytic and DFT results so it's clearer. – Pedro Gutiérrez Nov 29 '21 at 09:22
  • Your FFT is over half a Gaussian. You need to make sure that your “signal” contains the full bell curve. Here are some other Q&A that you might find useful: https://stackoverflow.com/q/52437002/7328782 , https://stackoverflow.com/q/49317834/7328782 , https://stackoverflow.com/q/61777788/7328782 – Cris Luengo Nov 29 '21 at 14:59

0 Answers0