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

Posted

技术标签:

【中文标题】如何从 2d FFT 计算波数域坐标【英文标题】:How to calculate wavenumber domain coordinates from a 2d FFT 【发布时间】:2011-11-01 22:47:15 【问题描述】:

我有一个 2d 复数数组,表示在真实空间中沿平面测量的势场。假设阵列是 128 个单元 x 128 个单元,平面的总面积是 500m x 500m。该数组中的每个单元格代表空间域中的一个点,坐标为 x 和 y。

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

【问题讨论】:

【参考方案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

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

【讨论】:

当我发现fftfreq()函数时真是太幸福了!【参考方案2】:

嗯,如果我使用 FFT 函数,DC 似乎会在零元素处结束,在你的情况下,频率之间的间距也是 1/500m。因此,以下(不是很简洁) sn-p 将获得您的频率轴:

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 成为以 rads/meter 为单位的空间频率数组,我会说

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

(假设我已经导入了 arange 和 pi 之类的东西)。

编辑

Stu 对高于奈奎斯特的频率提出了一个很好的观点,您可能更愿意将其视为负数(我通常这样做,但在代码中没有)。你总是可以这样做:

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

但是,如果您真的想要负频率,您可能还想玩 fftshift,以及 - 另一罐蠕虫。

【讨论】:

在我上面的回答中,[0][0] 处 kx 和 ky 的波数均为 0,表示该区域的 DC 平均值。您的频率未显示在 N/2 处切换到负频率。 我不认为我们在这里不同意。输出阵列的直流分量位于 [0][0]。 as 之间的主要区别似乎是我使用一维频率轴,而您生成二维网格网格。如果我正确理解您的代码。 假设输入网格在 x 和 y 上的大小相等,我将值计算为一个数组,假设输入网格在 x 和 y 中的大小相等,这对于我的数据并不总是正确的。但是 fft 的频率输出在 [0:n/2 -1] 的范围内从 0 增加到 2*pi*(n/2 -1) 但随后从 -2*pi*(n/2) 继续增加到-0 在 [n/2:n] 的范围内 您说的很对,如果输入代码不是正方形,您必须计算单独的 kx 和 ky 值(虽然您不需要网格,但它们可以很方便)。至于负频率与正频率,整个频率轴是采样率的模数,所以这是一个品味/方便/惯例的问题。我更喜欢保持轴连续,所以除非我也做了 fftshift,否则没有 -ive 频率;但这只是我。我猜你所说的 EM 场的意思是我们正在使用 DFT 来近似连续傅立叶积分(其中 -iv 频率确实是绝对的)。是的。 我对波数域中的数据应用的变换取决于正确计算 kx 和 ky。这包括当他们应该是消极的时候他们是消极的。我知道,当您尝试以对人们有意义的方式表示数据时,避免负频率会更容易,但这样方程式就不起作用了。【参考方案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

谢谢。

【讨论】:

以上是关于如何从 2d FFT 计算波数域坐标的主要内容,如果未能解决你的问题,请参考以下文章

如何从(0,0)开始以长度和角度获取2D坐标

如何从CoreLocation中的两个坐标获得对称坐标

如何从相机投影矩阵计算 ext 和 int 参数

如何有效地优化2D路径中的坐标数?

CoreML:iOS:如何获取检测模型中心的 2D 坐标

如何使用麦克风从麦克风/线路输入进行FFT