Python：如何插入“非结构化”二维傅里叶变换数据

Question

我的目标是对函数的离散连续二维傅里叶变换进行插值。问题似乎是每个维度中的频率不是严格按升序输出的（见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），以便在调用finterp(w_x, w_y)、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]? 处评估的函数？

Answer 1

以下代码显示如何使用fftshift更改fft2和fftfreq，以便频率轴在单调上增加。应用fftshift时，可以使用阵列进行插值。我添加了阵列的显示，以便您可以验证数据本身是否不变。原点从左上角移动到阵列的中间，从右侧向左侧移动负频率。

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计算DFT for k=[0..N-1]。由于周期性，在k 987654329 @和k-N和k-N，它的输出通常被解释为具有k=[N//2...(N-1)//2]（但是以不同的方式匹配k=[0..N-1]）;这是k 987654334 @返回（它返回987654335 @）。

因此，您可以说

N = f.shape[0]
w_x = np.linspace(0, N, N, endpoint=False) / N

现在您没有任何负频率，而是在[0,N-1]/N的范围内具有频率。