Python：如何对'unstructured'二维傅立叶变换数据进行插值
Python: How to interpolate 'unstructured' 2D Fourier transform data
我的目标是对函数的离散化连续二维傅立叶变换进行插值。问题似乎是每个维度的频率不是严格按升序输出的（参见here）。
fft.fft2 函数接受二维数组，在我的例子中，数组（我们称之为 A）的结构使得 A[i][j] = fun(x[i], y[j])、fun 是要转换的函数。将fft.fft2应用到A后，输出为与原数组维度相同的数组F，使得F[i][j]对应的频率坐标为(w_x[i], w_y[j]) ，其中 w_x = fft.fftfreq(F.shape[0]) 和 w_y = fft.fftfreq(F.shape[1])，它们都是非升序的一维数组。
超过 wx 和 wy 我想插值 F（比如函数 finterp），以便在调用 [=28 时返回插值=]、w_x 和 w_y 在 wx 的域和 wy 的范围内，但在其他方面是任意的。我研究了通过 scipy.interpolate 可用的各种插值，但在我看来，它们中的任何一个都不能处理这种类型的数据结构（坐标轴被定义为乱序一维数组和函数值在二维数组中）。
这有点抽象，所以我在这里做了一个简单的例子，结构与上面类似。假设我们希望在区域 x = [-1, 1] 和 y = [-1, 1] 上构造一个连续函数 f(x, y) = x + y，给定以下数据：
import numpy as np
# note that below z[i][j] corresponds to what we want f(x[i], y[j]) to be
x = np.array([0, 1, -1])
y = np.array([0, 1, -1])
z = np.array([0, 1, -1],[1, 2, 0],[-1, 0, -2])
z[i][j] 我们知道对应于在 x[i], y[j] 计算的函数。如何 (a) 直接插入此数据，给定其原始结构，或 (b) 重新排列数据，使 x 和 y 以升序排列，并且排列的 z是 z[i][j] 等于在重新排列的 x[i], y[j]?
处评估的函数
下面的代码展示了如何使用fftshift改变fft2和fftfreq的输出，使频率轴单调递增。应用 fftshift 后，您可以使用数组进行插值。我添加了数组显示，以便您可以验证数据本身是否未更改。原点从 top-left 角移到数组的中间，将负频率从右侧移到左侧。
import numpy as np
import matplotlib.pyplot as pp
x = np.array([0, 1, -1])
y = np.array([0, 1, -1])
z = np.array([[0, 1, -1],[1, 2, 0],[-1, 0, -2]])
f = np.fft.fft2(z)
w_x = np.fft.fftfreq(f.shape[0])
w_y = np.fft.fftfreq(f.shape[1])
pp.figure()
pp.imshow(np.abs(f))
pp.xticks(np.arange(0,len(w_x)), np.round(w_x,2))
pp.yticks(np.arange(0,len(w_y)), np.round(w_y,2))
f = np.fft.fftshift(f)
w_x = np.fft.fftshift(w_x)
w_y = np.fft.fftshift(w_y)
pp.figure()
pp.imshow(np.abs(f))
pp.xticks(np.arange(0,len(w_x)), np.round(w_x,2))
pp.yticks(np.arange(0,len(w_y)), np.round(w_y,2))
pp.show()
另一种方法是不使用 fftfreq 来确定您的频率，而是手动计算它们。默认情况下，FFT 计算 k=[0..N-1] 的 DFT。由于周期性，k 处的 DFT 等于 k+N 和 k-N 处的 DFT，其输出通常被解释为具有 k=[N//2...(N-1)//2]（但排列不同以匹配k=[0..N-1]);这是 fftfreq returns 的 k （它 returns k/N）。
因此，您可以改为说
N = f.shape[0]
w_x = np.linspace(0, N, N, endpoint=False) / N
现在您没有任何负频率，而是频率在 [0,N-1]/N.
范围内
我的目标是对函数的离散化连续二维傅立叶变换进行插值。问题似乎是每个维度的频率不是严格按升序输出的（参见here）。
fft.fft2 函数接受二维数组，在我的例子中，数组（我们称之为 A）的结构使得 A[i][j] = fun(x[i], y[j])、fun 是要转换的函数。将fft.fft2应用到A后，输出为与原数组维度相同的数组F，使得F[i][j]对应的频率坐标为(w_x[i], w_y[j]) ，其中 w_x = fft.fftfreq(F.shape[0]) 和 w_y = fft.fftfreq(F.shape[1])，它们都是非升序的一维数组。
超过 wx 和 wy 我想插值 F（比如函数 finterp），以便在调用 [=28 时返回插值=]、w_x 和 w_y 在 wx 的域和 wy 的范围内，但在其他方面是任意的。我研究了通过 scipy.interpolate 可用的各种插值，但在我看来，它们中的任何一个都不能处理这种类型的数据结构（坐标轴被定义为乱序一维数组和函数值在二维数组中）。
这有点抽象，所以我在这里做了一个简单的例子，结构与上面类似。假设我们希望在区域 x = [-1, 1] 和 y = [-1, 1] 上构造一个连续函数 f(x, y) = x + y，给定以下数据：
import numpy as np
# note that below z[i][j] corresponds to what we want f(x[i], y[j]) to be
x = np.array([0, 1, -1])
y = np.array([0, 1, -1])
z = np.array([0, 1, -1],[1, 2, 0],[-1, 0, -2])
z[i][j] 我们知道对应于在 x[i], y[j] 计算的函数。如何 (a) 直接插入此数据，给定其原始结构，或 (b) 重新排列数据，使 x 和 y 以升序排列，并且排列的 z是 z[i][j] 等于在重新排列的 x[i], y[j]?
下面的代码展示了如何使用fftshift改变fft2和fftfreq的输出，使频率轴单调递增。应用 fftshift 后，您可以使用数组进行插值。我添加了数组显示，以便您可以验证数据本身是否未更改。原点从 top-left 角移到数组的中间，将负频率从右侧移到左侧。
import numpy as np
import matplotlib.pyplot as pp
x = np.array([0, 1, -1])
y = np.array([0, 1, -1])
z = np.array([[0, 1, -1],[1, 2, 0],[-1, 0, -2]])
f = np.fft.fft2(z)
w_x = np.fft.fftfreq(f.shape[0])
w_y = np.fft.fftfreq(f.shape[1])
pp.figure()
pp.imshow(np.abs(f))
pp.xticks(np.arange(0,len(w_x)), np.round(w_x,2))
pp.yticks(np.arange(0,len(w_y)), np.round(w_y,2))
f = np.fft.fftshift(f)
w_x = np.fft.fftshift(w_x)
w_y = np.fft.fftshift(w_y)
pp.figure()
pp.imshow(np.abs(f))
pp.xticks(np.arange(0,len(w_x)), np.round(w_x,2))
pp.yticks(np.arange(0,len(w_y)), np.round(w_y,2))
pp.show()
另一种方法是不使用 fftfreq 来确定您的频率，而是手动计算它们。默认情况下，FFT 计算 k=[0..N-1] 的 DFT。由于周期性，k 处的 DFT 等于 k+N 和 k-N 处的 DFT，其输出通常被解释为具有 k=[N//2...(N-1)//2]（但排列不同以匹配k=[0..N-1]);这是 fftfreq returns 的 k （它 returns k/N）。
因此，您可以改为说
N = f.shape[0]
w_x = np.linspace(0, N, N, endpoint=False) / N
现在您没有任何负频率，而是频率在 [0,N-1]/N.