CODEWITHSUNDEEP
pythonarraysnumpyfft
cstahl
cstahl

Reputation: 43

Numpy FFT fails on meshgrids?

I'm running into an odd problem using Numpy meshgrids with the FFT functions. Specifically, either the fft2 or the ifft2 function seems to fail when used on an array built using meshgrids. 

x = np.arange(-4, 4, .08)
y = np.arange(-4, 4, .08)
X, Y = np.meshgrid(x, y)
field = (X + i*Y)*np.exp(X**2 + Y**2)

As a check before I proceeded with my project, I did 

fieldCheck1 = np.fft.fft2(field)
fieldCheck2 = np.fft.ifft2(fieldCheck1)

which should yield back my original array, but in fact eliminates the real part ( a plot of abs(fieldCheck2)**2 is flat zero, where originally it was a gaussian) and completely scrambles the phase information (a phase plot of fieldCheck2 looks like static instead of a phase ramp)

I've checked the documentation, but I don't see anything in there that would explain this. Any insight into the source of the problem would be helpful.

Upvotes: 4

Views: 543

Answers (2)

SleuthEye
SleuthEye

Reputation: 14577

To add to the problem described by @jakevdp, it looks like you go into this whole mess because your function (after replacing ^ with ** in your code)

field = (X + 1.0j*Y)*np.exp(X**2 + Y**2)

which you describe as being a Gaussian is in fact not quite so, as illustrated in this graph:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,np.abs(field));

enter image description here

To get an actual 2D Gaussian function, try:

x = np.arange(-4, 4, .08)
y = np.arange(-4, 4, .08)
X, Y = np.meshgrid(x, y)
field = np.exp(-(X**2 + Y**2))  # notice "-" sign in the exponent

Which would then give you the graph:

enter image description here

And corresponding error np.abs(fieldCheck2-field) for the round trip transform as (which is more in line with the expected numerical precision):

enter image description here

Upvotes: 0

jakevdp
jakevdp

Reputation: 86433

The problem (after replacing ^ with ** in your code) is that the contrast between your smallest and largest values is nearly 30 orders of magnitude:

>>> abs(field).max() / abs(field).min()
8.8904389513698014e+28

Floating point arithmetic only has finite precision, so identities that work for real numbers will not always work for floating point numbers. As a simple example:

>>> x = 1
>>> y = 1e30
>>> z = x + y

>>> x == z - y
False

The FFT is essentially a more complicated version of this: you're adding very small numbers to very big numbers, and when you subtract the very big numbers again you get zero rather than the small number you expect.

Upvotes: 3

Related Questions