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

I am trying to utilize Numpy's fft function, however when I give the function a simple gausian function the fft of that gausian function is not a gausian, its close but its halved so that each half is at either end of the x axis.

The Gaussian functi开发者_开发知识库on I'm calculating is y = exp(-x^2)

Here is my code:

from cmath import *
from numpy import multiply
from numpy.fft import fft
from pylab import plot, show

""" Basically the standard range() function but with float support """
def frange (min_value, max_value, step):
    value = float(min_value)
    array = []
    while value < float(max_value):
        array.append(value)
        value += float(step)
    return array


N = 256.0 # number of steps
y = []
x = frange(-5, 5, 10/N)

# fill array y with values of the Gaussian function   
cache = -multiply(x, x)
for i in cache: y.append(exp(i))

Y = fft(y)

# plot the fft of the gausian function
plot(x, abs(Y))
show()

The result is not quite right, cause the FFT of a Gaussian function should be a Gaussian function itself...

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()

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.

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.

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).

Following on from Sven Marnach's answer, a simpler version would be this:

from pylab import *
N = 128

x = ifftshift(arange(-5,5,5./N))

y = exp(-x*x)
y_fft = fft(y) / sqrt(2 * N)

plot(fftshift(y))
plot(fftshift(y_fft))

show()

This yields a plot identical to the above one.

The key (and this seems strange to me) is that NumPy's assumed data ordering --- in both frequency and time domains --- is to have the "zero" value first. This is not what I'd expect from other implementations of FFT, such as the FFTW3 libraries in C.

This was slightly fudged in the answers from unutbu and Steve Tjoa above, because they're taking the absolute value of the FFT before plotting it, thus wiping away the phase issues resulting from not using the "standard order" in time.