二维 FFT 显示高于奈奎西特极限的意外频率
Two dimensional FFT showing unexpected frequencies above Nyquisit limit
注意：这个问题是建立在我的另一个问题之上的：
我有一些数据，基本上是函数 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.
如果我的信号确实具有高于奈奎西特限制的频率，它们将 "folded" 回到奈奎西特域，导致错误结果（混叠）。但是即使我可能以太低的频率采样，理论上应该不可能看到任何超过 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 方向的正弦波）。
好的，我们开始吧！这里有几个简单的函数和一个您可以使用的完整示例：它有一些与绘图和数据生成相关的额外内容，但第一个函数 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)，这正是纯复指数频率下频谱的正确值。
注意：这个问题是建立在我的另一个问题之上的：
我有一些数据，基本上是函数 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.
如果我的信号确实具有高于奈奎西特限制的频率，它们将 "folded" 回到奈奎西特域，导致错误结果（混叠）。但是即使我可能以太低的频率采样，理论上应该不可能看到任何超过 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 方向的正弦波）。
好的，我们开始吧！这里有几个简单的函数和一个您可以使用的完整示例：它有一些与绘图和数据生成相关的额外内容，但第一个函数 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)，这正是纯复指数频率下频谱的正确值。