Two dimensional FFT showing unexpected frequencies above Nyquisit limit

Question

Note: This question is building on another question of mine: Two dimensional FFT using python results in slightly shifted frequency

I have some data, basically a function E(x,y) with (x,y) being a (discrete) subset of R^2, mapping to real numbers. For the (x,y) plane i have a fixed distance between data points in x- as well as in y direction (0,2). I want to analyze the frequency spectrum of my E(x,y) signal using a two dimensional fast fourier transform (FFT) using python.

As far as i know, no matter which frequencies are actually contained in my signal, using FFT, i will only be able to see signals below the Nyquisit limit Ny, which is Ny = sampling frequency / 2. In my case i have a real spacing of 0,2, leading to a sampling frequency of 1 / 0,2 = 5 and therefore my Nyquisit limit is Ny = 5 / 2 = 2,5.

If my signal does have frequencies above the Nyquisit limit, they will be "folded" back into the Nyquisit domain, leading to false results (aliasing). But even though i might sample with a too low frequency, it should in theory not be possible to see any frequencies above the Niquisit limit, correct?

So here is my issue: Analyzing my signal should only lead to frequencies of 2,5 max., but i cleary get frequencies higher than that. Given that i am pretty sure about the theory here, there has to be some mistake in my code. I will provide a shortened code version, only providing necessary information for this issue:

simulationArea =...  # length of simulation area in x and y direction
x = np.linspace(0, simulationArea, numberOfGridPointsInX, endpoint=False)
y = x
xx, yy = np.meshgrid(x, y)
Ex = np.genfromtxt('E_field_x100.txt')  # this is the actual signal to be analyzed, which may have arbitrary frequencies
FTEx = np.fft.fft2(Ex)  # calculating fft coefficients of signal
dx = x[1] - x[0]  # calculating spacing of signals in real space. 'print(dx)' results in '0.2'

sampleFrequency = 1.0 / dx
nyquisitFrequency = sampleFrequency / 2.0
half = len(FTEx) / 2

fig, axarr = plt.subplots(2, 1)

im1 = axarr[0, 0].imshow(Ex,
                         origin='lower',
                         cmap='jet',
                         extent=(0, simulationArea, 0, simulationArea))
axarr[0, 0].set_xlabel('X', fontsize=14)
axarr[0, 0].set_ylabel('Y', fontsize=14)
axarr[0, 0].set_title('$E_x$', fontsize=14)
fig.colorbar(im1, ax=axarr[0, 0])

im2 = axarr[1, 0].matshow(2 * abs(FTEx[:half, :half]) / half,
                          aspect='equal',
                          origin='lower',
                          interpolation='nearest')
axarr[1, 0].set_xlabel('Frequency wx')
axarr[1, 0].set_ylabel('Frequency wy')
axarr[1, 0].xaxis.set_ticks_position('bottom')
axarr[1, 0].set_title('$FFT(E_x)$', fontsize=14)
fig.colorbar(im2, ax=axarr[1, 0])

The result of this is:

How is that possible? When i am using the same code for very simple signals, it works just fine (e.g. a sine wave in x or y direction with a specific frequency).

Answer 1

Ok here we go! Here’s a couple of simple functions and a complete example that you can use: it’s got a little bit of extra cruft related to plotting and for data generation but the first function, makeSpectrum shows you how to use fftfreq and fftshift plus fft2 to achieve what you want. Let me know if you have questions.

import numpy as np
import numpy.fft as fft
import matplotlib.pylab as plt


def makeSpectrum(E, dx, dy, upsample=10):
    """
    Convert a time-domain array `E` to the frequency domain via 2D FFT. `dx` and
    `dy` are sample spacing in x (left-right, 1st axis) and y (up-down, 0th
    axis) directions. An optional `upsample > 1` will zero-pad `E` to obtain an
    upsampled spectrum.

    Returns `(spectrum, xf, yf)` where `spectrum` contains the 2D FFT of `E`. If
    `Ny, Nx = spectrum.shape`, `xf` and `yf` will be vectors of length `Nx` and
    `Ny` respectively, containing the frequencies corresponding to each pixel of
    `spectrum`.

    The returned spectrum is zero-centered (via `fftshift`). The 2D FFT, and
    this function, assume your input `E` has its origin at the top-left of the
    array. If this is not the case, i.e., your input `E`'s origin is translated
    away from the first pixel, the returned `spectrum`'s phase will *not* match
    what you expect, since a translation in the time domain is a modulation of
    the frequency domain. (If you don't care about the spectrum's phase, i.e.,
    only magnitude, then you can ignore all these origin issues.)
    """
    zeropadded = np.array(E.shape) * upsample
    F = fft.fftshift(fft.fft2(E, zeropadded)) / E.size
    xf = fft.fftshift(fft.fftfreq(zeropadded[1], d=dx))
    yf = fft.fftshift(fft.fftfreq(zeropadded[0], d=dy))
    return (F, xf, yf)


def extents(f):
    "Convert a vector into the 2-element extents vector imshow needs"
    delta = f[1] - f[0]
    return [f[0] - delta / 2, f[-1] + delta / 2]


def plotSpectrum(F, xf, yf):
    "Plot a spectrum array and vectors of x and y frequency spacings"
    plt.figure()
    plt.imshow(abs(F),
               aspect="equal",
               interpolation="none",
               origin="lower",
               extent=extents(xf) + extents(yf))
    plt.colorbar()
    plt.xlabel('f_x (Hz)')
    plt.ylabel('f_y (Hz)')
    plt.title('|Spectrum|')
    plt.show()


if __name__ == '__main__':
    # In seconds
    x = np.linspace(0, 4, 20)
    y = np.linspace(0, 4, 30)
    # Uncomment the next two lines and notice that the spectral peak is no
    # longer equal to 1.0! That's because `makeSpectrum` expects its input's
    # origin to be at the top-left pixel, which isn't the case for the following
    # two lines.
    # x = np.linspace(.123 + 0, .123 + 4, 20)
    # y = np.linspace(.123 + 0, .123 + 4, 30)

    # Sinusoid frequency, in Hz
    x0 = 1.9
    y0 = -2.9

    # Generate data
    im = np.exp(2j * np.pi * (y[:, np.newaxis] * y0 + x[np.newaxis, :] * x0))

    # Generate spectrum and plot
    spectrum, xf, yf = makeSpectrum(im, x[1] - x[0], y[1] - y[0])
    plotSpectrum(spectrum, xf, yf)

    # Report peak
    peak = spectrum[:, np.isclose(xf, x0)][np.isclose(yf, y0)]
    peak = peak[0, 0]
    print('spectral peak={}'.format(peak))

Results in the following image, and prints out, spectral peak=(1+7.660797103157986e-16j), which is exactly the correct value for the spectrum at the frequency of a pure complex exponential.