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!