I am trying to perform numerical Fourier transforms using numpy.fft.fft. As a test function, I used a rather simple function: a Gaussian with a complex width parameter
.
The Fourier transform of this function can easily be calculated analytically and gives
.
My problem is that I cannot reproduce the result of the analytically Fourier transform numerically using numpy.fft.fft. Here is my code
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.fftpack import fft, ifft, fftshift, fftfreq, fft2, ifft2
def fun(t, alpha, beta):
# a function of time and 2 parameters
norm = np.sqrt(alpha/np.pi)
return norm * np.exp( -t**2 * (alpha + 1j*beta) )
def ft_fun_analytic(om, alpha, beta):
# the analytical Fourier transform of fun(t)
gamma = 1. / (1. + 1j*beta/alpha)
norm = np.sqrt( gamma )
return norm * np.exp( -om**2 * gamma / (4.*alpha) )
t = np.linspace(-4., 4., 300000, endpoint=True)
om = 2.*np.pi * fftshift(fftfreq(len(t), d = t[1]-t[0]))
alpha, beta = [160., 0.07]
t_series = fun(t, alpha, beta)
grid_norm = (t[1]-t[0])
om_spec = grid_norm * fftshift(fft(t_series))
om_spec_analytic = ft_fun_analytic(om, alpha, beta)
fig1 = plt.figure(figsize=(10,5))
plt.subplot(121, xlabel=r'$t$', ylabel='Intensity', title=r'$f(t)$')
plt.plot(t, np.abs(t_series)**2)
plt.xlim([-0.4,0.4])
plt.subplot(122, xlabel=r'$t$', ylabel='Phase', title=r'$f(t)$')
plt.plot(t, np.angle(t_series))
plt.xlim([-0.4,0.4])
plt.show()
fig2 = plt.figure(figsize=(10,5))
plt.subplot(121, xlabel=r'$\omega$', ylabel='Intensity', title=r'$FT[f](\omega)$')
plt.plot(om, np.abs(om_spec)**2, label='numerical FT')
plt.plot(om, np.abs(om_spec_analytic)**2, '--', label='analytical FT')
plt.xlim([-120., 120.])
plt.legend()
plt.subplot(122, xlabel=r'$\omega$', ylabel='Phase', title=r'$FT[f](\omega)$')
plt.plot(om, np.angle(om_spec), label='numerical FT')
plt.plot(om, np.angle(om_spec_analytic), '--', label='analytical FT')
plt.xlim([-120., 120.])
plt.legend()
plt.show()
The function in the time domain looks like this:
In the frequency domain, the numerical Fourier transform gives the correct intensity. But the phase includes an offset and high frequency oscillation:
I am aware of effects like aliasing and problems with the sampling range for fast Fourier transforms. I tried to play with both the resolution and range. Even for extemely fine sampling that my machine can only just handle this problem occurs.
Is there a way to reproduce the analytical result using numpy.fft.fft? Maybe by choosing a smart grid instead of a really fine one?
I am not 100% sure if this question belongs on the scientific computing site or here. Since my (possibly unfounded) suspicion is that this is about numerical precision in the evaluation of the exponential function, I will try here first.
