二维 FFT 显示超出 Nyquisit 极限的意外频率

Question

注意：这个问题是基于我的另一个问题： Two dimensional FFT using python results in slightly shifted frequency

我有一些数据，基本上是一个函数 E(x,y)，其中 (x,y) 是 R^2 的（离散）子集，映射到实数。对于 (x,y) 平面，我在 x- 和 y 方向 (0,2) 上的数据点之间具有固定距离。我想使用 python 使用二维快速傅立叶变换 (FFT) 分析我的 E(x,y) 信号的频谱。

据我所知，无论我的信号中实际包含哪些频率，使用 FFT，我只能看到低于 Nyquisit 限制 Ny 的信号，即 Ny = 采样频率 / 2。就我而言，我实际间距为 0,2，导致采样频率为 1 / 0,2 = 5，因此我的 Nyquisit 限制为 Ny = 5 / 2 = 2,5。

如果我的信号的频率确实高于 Nyquisit 限制，它们将被“折叠”回 Nyquisit 域，从而导致错误结果（混叠）。但即使我可能会以太低的频率进行采样，但理论上应该不可能看到任何高于 Niquisit 限制的频率，对吗？

所以这是我的问题：分析我的信号只能导致最大 2.5 的频率，但我显然会得到更高的频率。鉴于我对这里的理论非常确定，我的代码中肯定有一些错误。我将提供一个缩短的代码版本，仅提供此问题的必要信息：

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

这样的结果是：

这怎么可能？当我对非常简单的信号使用相同的代码时，它工作得很好（例如，具有特定频率的 x 或 y 方向的正弦波）。

Answer 1

好的，我们开始吧！这里有几个简单的函数和一个完整的例子，你可以使用：它有一点额外的与绘图和数据生成相关的麻烦，但第一个函数 makeSpectrum 向你展示了如何使用 fftfreq 和 fftshift加上fft2 来实现你想要的。如果您有任何问题，请告诉我。

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

结果如下图，并打印出spectral peak=(1+7.660797103157986e-16j)，这正是纯复指数频率处频谱的正确值。