如何从2D FFT计算出波数域坐标

Question

我有一个由复数组成的二维数组，表示在现实空间中沿平面测量的潜力场。假设该数组为128个单元格乘以128个单元格，平面的总面积为500米乘以500米。该数组中的每个单元格都代表空间域中具有x和y坐标的点。

当我对这个二维数组使用scipy.fftpack中的2d FFT时，我得到了相同信息在波域中的表示。如何计算输出数组中的点的波域坐标kx和ky？

- Stu

Answer 1

下面是一些代码，完整演示了问题和我能够找到的解决方案。

from numpy import linspace , arange , reshape ,zeros
from scipy.fftpack import fft2 , fftfreq
from cmath import pi

# create some arbitrary data
some_data = arange(0.0 , 16384.0 , dtype = complex)

# reshape it to be a 128x128 2d grid
some_data_grid = reshape(some_data , (128 , 128) )

# assign some real spatial co-ordinates to the grid points   
# first define the edge values
x_min = -250.0
x_max = 250.0
y_min = -250.0
y_max = 250

# then create some empty 2d arrays to hold the individual cell values
x_array = zeros( (128,128) , dtype = float )
y_array = zeros( (128,128) , dtype = float )

# now fill the arrays with the associated values
for row , y_value in enumerate(linspace (y_min , y_max , num = 128) ):

  for column , x_value in enumerate(linspace (x_min , x_max , num = 128) ):

    x_array[row][column] = x_value
    y_array[row][column] = y_value

# now for any row,column pair the x_array and y_array hold the spatial domain
# co-ordinates of the associated point in some_data_grid

# now use the fft to transform the data to the wavenumber domain
some_data_wavedomain = fft2(some_data_grid)

# now we can use fftfreq to give us a base for the wavenumber co-ords
# this returns [0.0 , 1.0 , 2.0 , ... , 62.0 , 63.0 , -64.0 , -63.0 , ... , -2.0 , -1.0 ]
n_value = fftfreq( 128 , (1.0 / 128.0 ) )

# now we can initialize some arrays to hold the wavenumber co-ordinates of each cell
kx_array = zeros( (128,128) , dtype = float )
ky_array = zeros( (128,128) , dtype = float )

# before we can calculate the wavenumbers we need to know the total length of the spatial
# domain data in x and y. This assumes that the spatial domain units are metres and
# will result in wavenumber domain units of radians / metre.
x_length = x_max - x_min
y_length = y_max - y_min

# now the loops to calculate the wavenumbers
for row in xrange(128):

  for column in xrange(128):

    kx_array[row][column] = ( 2.0 * pi * n_value[column] ) / x_length
    ky_array[row][column] = ( 2.0 * pi * n_value[row] ) / y_length

# now for any row,column pair kx_array , and ky_array will hold the wavedomain coordinates
# of the correspoing point in some_data_wavedomain

我知道这可能不是最有效的方法，但希望它易于理解。我希望这能帮助某人避免一些挫败感。

- Stu

---

当我发现fftfreq()函数时，感觉非常幸福！ - PhilMacKay

Answer 2

嗯，看起来如果我使用FFT函数，那么直流信号会出现在零元素位置上，还有在你的情况下频率之间的间距是1/500m。因此，以下不太简洁的代码片段将给你获取频率轴的方法：

meter = 1.0
L     = 500.0 * meter
N     = 128

dF    = 1.0 / L
freqs = arange(0, N/L, dF)  # array of spatial frequencies.

自然地，这些频率以每米循环为单位，而不是弧度/米。如果我想让kx和ky成为弧度/米的空间频率数组，我只需要这样说：

 kx = 2*pi*freqs
 ky = 2*pi*freqs

（假设我已经导入了arange和pi等内容）。

编辑

Stu提出了一个关于超过奈奎斯特频率的频率的好观点，你可能更喜欢将其视为负数（我通常这样做，但代码中不是这样）。你可以随时执行以下操作：

freqs[freqs > 0.5*N/(2*L)] -= N/L

但是如果你真的想要负频率，你可能也想尝试一下fftshift -- 另一个棘手的问题。

- Adrian Ratnapala

---

在我上面的回答中，kx和ky在[0][0]处的波数都为0，表示该区域的直流平均值。你的频率没有显示在N/2处切换到负频率。 - Stu

我认为我们在这里没有分歧。输出数组的DC分量位于[0][0]。似乎主要的区别在于我使用一维频率轴，而你生成二维网格。如果我正确理解了你的代码。 - Adrian Ratnapala

当输入网格在x和y方向上的大小相等时，我将值计算为数组，但这并不总是适用于我的数据。但是，fft的频率输出从[0：n/2-1]范围内的0增加到2pi(n/2-1)，然后在[n/2:n]范围内继续从-2pi(n/2)到-0。 - Stu

如果你的输入空间域数据表示为复数（例如电磁场的情况），那么你需要将负频率视为实数据。 - Stu

你说得很对，如果输入的代码不是正方形的，你必须计算单独的kx和ky值（虽然你不需要网格，但它们可能很方便）。至于负频率与正频率，整个频率轴都是模采样率的，所以这是一种品味/方便/惯例问题。我更喜欢保持我的轴连续，所以没有负频率，除非我也做了fftshift；但那只是我的个人喜好。我猜你所说的EM场是指我们使用DFT来近似连续的傅里叶积分（其中-iv频率确实是绝对的）。没错。 - Adrian Ratnapala

我在波数域中对数据应用的变换取决于正确计算kx和ky。这包括它们在应该为负时为负。我知道当你试图以使人们理解的方式表示数据时，避免负频率会更容易，但是这样方程就无法工作。 - Stu

Answer 3

我在下面的FORTRAN代码中找到了子程序，并尝试将其实现为Matlab函数。请给我评论，告诉我我的翻译是否正确。 这是FORTRAN子程序：

subroutine kvalue(i,j,nx,ny,dkx,dky,kx,ky)
  c Subroutine KVALUE finds the wavenumber coordinates of one
  c element of a rectangular grid from subroutine FOURN.
  c
  c Input parameters:
  c i - index in the ky direction,
  c j - index in the kx direction.
  c nx - dimension of grid in ky direction (a power of two).
  c ny - dimension of grid in kx direction (a power of two).
  c dkx - sample interval in the kx direction,
  c dky - sample interval in the ky direction,
  c
  c Output parameters:
  c kx - the wavenumber coordinate in the kx direction,
  c ky - the wavenumber coordinate in the ky direction,
  c
  real kx,ky
  nyqx=nx/2+l
  nyqy=ny/2+l
  if(j.le.nyqx)then
      kx=(j-l)*dkx
  else
      kx=(j-nx-l)*dkx
  end if
  if(i.le.nyqy)then
          ky=(i-l)*dky
  else
          ky=(i-ny-l)*dky
  end if
  return
  end

我将其翻译为MATLAB函数：

function [kx,ky]=kvalue(gz,nx,ny,dkx,dky)
nyq_x=nx/2+1;
nyq_y=ny/2+1;
for i=1:length(gz)
    for j=1:length(gz)
        if j <= nyq_x
            kx(j)=(j-1)*dkx;
        else
            kx(j)=(j-nx-1)*dkx;
        end
        if i <= nyq_y
            ky(i)=(i-1)*dky;
        else
            ky(i)=(i-ny-1)*dky;
        end
    end
end

谢谢。

- Ferry Aries