二维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，我只能看到低于Nyquisit极限Ny的信号，即Ny =采样频率/ 2。在我的情况下我有一个0.2的实际间距，导致采样频率为1/0,2 = 5，因此我的Nyquisit限制是Ny = 5/2 = 2,5。

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

所以这里是我的问题：分析我的信号应该只会导致频率最高2.5，但我cleary得到更高的频率。鉴于我对这里的理论非常肯定，我的代码中必然存在一些错误。我公司将提供缩短代码版本，只为这个问题提供必要的信息：

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)，这正是用于在纯复指数的频率频谱中的正确的值。