Understanding a simple case of fft2 in python

Question

I have 2D arrays representing surfaces. I have created them by scanning some theoretically flat objects and I ended up with what can be called "imperfections", or else, deviation from the perfectly flat surface.

I hope to get information on the dominant wavelengths of imperfections and their amplitude, so I employed numpy's fft2. Since I am a first timer, I did some digging and tried to play around with some basic examples to get a grip.

Here is one:

import numpy as np
import matplotlib.pyplot as plt

"""
fft2 playground.
"""

# Initialise an empty array
field = np.empty([29, 481])

# Set the wave amplitude
amp = 2.5

# Create the synthetic sinusoidal field
for n, i in enumerate(np.linspace(0, 2 * np.pi, field.shape[0])):
    for m, j in enumerate(np.linspace(0, 2 * np.pi, field.shape[1])):
        field[n, m] = amp * np.sin(j)

# Perform 2D fourier and shift the result to centre
f = np.fft.fft2(field)
fshift = np.fft.fftshift(f)

# Calculate the magnitude and phase spectra
magnitude_spectrum = 20*np.log(np.abs(fshift))
phase_spectrum = np.angle(fshift)

# Reconstruct the initial field
f_ishift = np.fft.ifftshift(fshift)
re_field = np.abs(np.fft.ifft2(f_ishift))

# Plot
fig = plt.figure()

fig.add_subplot(411)
plt.imshow(field, cmap='gray')
plt.title('Field'), plt.xticks([]), plt.yticks([])
plt.colorbar()

fig.add_subplot(412)
plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude spectrum'), plt.xticks([]), plt.yticks([])

fig.add_subplot(413)
plt.imshow(phase_spectrum, cmap='gray')
plt.title('Phase spectrum'), plt.xticks([]), plt.yticks([])

fig.add_subplot(414)
plt.imshow(re_field, cmap='gray')
plt.title('Reconstructed field'), plt.xticks([]), plt.yticks([])
plt.colorbar()

plt.show()

Which gives:

Obviously the returned field is not the same as the input. The amplitude is correct but it is mirrored and 90 degrees off-phase. I feel I am missing something but I cannot nail it.

Also, I would expect that the magnitude spectrum would be something like 2 dots, one on either side of the central pixel and nothing else. What I understand is that the central pixel is the zero-frequency and moving to the edges are the higher frequencies, up to the Nyquist.

Worth noting, if I change the dimensions of the input field I get different results on the magnitude and phase spectra.

Could someone point out what am I missing please? How would I get back the initial field from the fft and why the magnitude spectrum is the way it is?

Thanx!

Answer 1

When you construct re_field, ditch the np.abs. It turns all negative values of the sine wave positive. Instead, do a np.realon the result of the ifft2 to get rid of the imaginary parts caused by rounding errors. Then you will get a result perfectly matching the original.

In the magnitude spectrum you cannot see the sine wave you inserted, because its frequency is to low. The vertical black line is the DC part in x direction, which is zero, the horizontal white line is the DC part in y direction, which has full ampitude because there is only DC in y direction.

Try some higher frequencies, then you will get a better idea.

Answer 2

The amplitude is correct but it is mirrored and 90 degrees off-phase. I feel I am missing something but I cannot nail it.

re_field has np.abs in it, so it's not the same thing. You've rectified the signal. To remove the negligible imaginary parts, use np.fft.ifft2(f_ishift).real instead.

Also, I would expect that the magnitude spectrum would be something like 2 dots, one on either side of the central pixel and nothing else.

It is, if you don't log-scale it first. Your frequency is 1 cycle per frame, though, which is very low, and will be right to the DC line. If you use a higher frequency, you can see the wider-spaced dots corresponding to that frequency.