python - 二维 FFT 显示高于奈奎西极限的意外频率

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，我只能看到低于奈奎西特极限 Ny 的信号，即 Ny = 采样频率/2。在我的情况下，我实际间距为 0,2，导致采样频率为 1/0,2 = 5，因此我的奈奎西特极限为 Ny = 5/2 = 2,5。

如果我的信号的频率确实高于奈奎西特限制，它们将被“折叠”回奈奎西特域，从而导致错误结果(混叠)。但即使我可能以太低的频率进行采样，理论上应该不可能看到任何高于 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))

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