Numpy FFT fails on meshgrids?

Question

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.

Is i the imaginary unit in your environment somehow? It is 1.0j in mine. — norio
– norio, Commented Feb 18, 2016 at 20:43
yes, I set it that way at the beginning of all my code. Old bad habits are easier to work around than change. — cstahl
– cstahl, Commented Feb 18, 2016 at 21:01
Bitwise-or (^) will fail when used between a float and an int. Are you trying to exponentiate (**)? — jakevdp
– jakevdp, Commented Feb 18, 2016 at 23:58
Drat, I know, and I have it correct in my actual code--mistyped for here. I've corrected it in the original post. — cstahl
– cstahl, Commented Feb 19, 2016 at 14:51

jakevdp · Accepted Answer · 2016-02-19 00:08:35Z

3

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.

answered Feb 19, 2016 at 0:08

jakevdp

88.6k15 gold badges140 silver badges177 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

cstahl Over a year ago

This was the problem exactly. I was able to adjust some coefficients (not shown in my simplified code in the question) to reduce range and now everything works as expected.

SleuthEye · Accepted Answer · 2016-02-19 01:01:42Z

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

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:

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

Collectives™ on Stack Overflow

Numpy FFT fails on meshgrids?

2 Answers 2

1 Comment

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

Comments

Your Answer

Sign up or log in

Post as a guest

Related