Fourier transform of a Gaussian is not a Gaussian, but thats wrong! - Python

np.fft.fft returns a result in so-called "standard order": (from the docs)

If A = fft(a, n), then A[0] contains the zero-frequency term (the mean of the signal), which is always purely real for real inputs. Then A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

The function np.fft.fftshift rearranges the result into the order most humans expect (and which is good for plotting):

The routine np.fft.fftshift(A) shifts transforms and their frequencies to put the zero-frequency components in the middle...

So using np.fft.fftshift:

import matplotlib.pyplot as plt
import numpy as np

N = 128
x = np.arange(-5, 5, 10./(2 * N))
y = np.exp(-x * x)
y_fft = np.fft.fftshift(np.abs(np.fft.fft(y))) / np.sqrt(len(y))
plt.plot(x,y)
plt.plot(x,y_fft)
plt.show()

A fourier transform implicitly repeats indefinitely, as it is a transform of a signal that implicitly repeats indefinitely. Note that when you pass y to be transformed, the x values are not supplied, so in fact the gaussian that is transformed is one centred on the median value between 0 and 256, so 128.

Remember also that translation of f(x) is phase change of F(x).

It is being displayed with the center (i.e. mean) at coefficient index zero. That is why it appears that the right half is on the left, and vice versa.

EDIT: Explore the following code:

import scipy
import scipy.signal as sig
import pylab
x = sig.gaussian(2048, 10)
X = scipy.absolute(scipy.fft(x))
pylab.plot(x)
pylab.plot(X)
pylab.plot(X[range(1024, 2048)+range(0, 1024)])

The last line will plot X starting from the center of the vector, then wrap around to the beginning.

Your result is not even close to a Gaussian, not even one split into two halves.

To get the result you expect, you will have to position your own Gaussian with the center at index 0, and the result will also be positioned that way. Try the following code:

from pylab import *
N = 128
x = r_[arange(0, 5, 5./N), arange(-5, 0, 5./N)]
y = exp(-x*x)
y_fft = fft(y) / sqrt(2 * N)
plot(r_[y[N:], y[:N]])
plot(r_[y_fft[N:], y_fft[:N]])
show()

The plot commands split the arrays in two halfs and swap them to get a nicer picture.