FFT of Gaussian and Lorentzian functions in Python

Question

I need to calculate the Fourier Transform of both Gaussian and Lorentzian functions and plot the result. I know the result for both, but I can't seem to get them right and I don't really understand how the fft works in Python...

This is what I have for the lorentzian function so far

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt

a = 1.0
N = 500 #number of points
x = np.linspace(-5,5,N)
lorentz = (a/np.pi) * (1/(a**2 + x**2)) #lorentzian function

fourier = fftpack.rfft(lorentz)

fourier = (2 * np.pi**2 / N)* abs(fourier[0:N/4]) #supposed to normalize

plt.plot(fourier)
plt.show()

I'm really lost on this. I should be getting something like exp(-pi|k|), I don't know what I should do with this to make that look more like what I want.

Answer 1

There are a number of reasons why the fft() function doesn't give the result you were expecting:

• Firstly, the fft() function only computes a discrete Fourier transform, which is defined as a summation over a finite number of regularly spaced points, and implicitly requires a periodic function. The continuous-time Fourier transform is defined as an integral over an infinite extent of time, without making assumptions of the signal being periodic.
• Secondly, the way you're defining your lorentz time-series effectively makes this centred on a time that is shifted away from the origin. This will introduce a phase-rotation to all the Fourier coefficients (which may be why you took the absolute value before generating your plot?).

If you take sufficiently large numbers of sample points, and extend the time-extent of your input waveform sufficiently, I believe you should indeed converge to the "correct" answer, e.g. as obtained by contour integration. For example, if you use

x = numpy.linspace(-50, 50, 5000)

you get a discrete Fourier transform that looks more like: