二維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)，這正是用於在純復指數的頻率頻譜中的正確的值。