放大 np.fft2 结果

Question

有没有办法从 np.fft2 选择 x/y 输出轴范围？ 我有一段代码计算光圈的衍射图案。孔径以 2k x 2k 像素阵列定义。衍射图案基本上是孔径 2D FT 的内部部分。 np.fft2 为我提供了一个与输入大小相同但具有一些预设 x/y 轴范围的输出数组。当然我可以使用图像查看器放大，但我已经丢失了细节。解决方法是什么？

谢谢， 格特

import numpy             as np
import matplotlib.pyplot as plt

r= 500
s= 1000

y,x = np.ogrid[-s:s+1, -s:s+1]
mask = x*x + y*y <= r*r
aperture = np.ones((2*s+1, 2*s+1))
aperture[mask] = 0

plt.imshow(aperture)
plt.show()

ffta= np.fft.fft2(aperture)

plt.imshow(np.log(np.abs(np.fft.fftshift(ffta))**2))
plt.show()

Answer 1

不幸的是，FFT 的大部分速度和精度来自与输入大小相同的输出。

增加输出傅立叶域中表观分辨率的传统方法是对输入进行零填充：

np.fft.fft2(aperture, [4 * (2*s+1), 4 * (2*s+1)])

告诉 FFT 将输入填充为

4 * (2*s+1)

像素高和宽，即，使输入大四倍（像素数的十六倍）。

从一开始我说“表观”分辨率是因为您拥有的实际数据量没有增加，但傅里叶变换会显得更平滑，因为输入域中的零填充会导致傅里叶变换对输出进行插值。在上面的示例中，可以用一个像素看到的任何特征都将用四个像素显示。为了使这个完全具体，这个例子表明零填充 FFT 的每四个像素在数值上与原始未填充 FFT 的每个像素相同：

# Generate your `ffta` as above, then
N = 2 * s + 1
Up = 4
fftup = np.fft.fft2(aperture, [Up * N, Up * N])

relerr = lambda dirt, gold: np.abs((dirt - gold) / gold)
print(np.max(relerr(fftup[::Up, ::Up] , ffta))) # ~6e-12.

（这

relerr

只是一个简单的相对误差，你想要接近机器精度，大约

2e-16

。零填充FFT和未填充FFT的每4个样本之间的最大误差是

6e-12

这非常接近机器精度，这意味着这两个数组在数值上几乎是等效的。）结束

零填充是解决问题的最直接方法。但它确实会花费你 lot 的内存。这令人沮丧，因为您可能只关心转换的一小部分。有一种称为 chirp z 变换（CZT，或通俗地称为“缩放 FFT”）的算法可以执行此操作。如果您的输入是

N

（对您

2*s+1

）并且您只想在任何地方评估FFT输出的

M

样本，它将计算三个大小为

N + M - 1

的傅立叶变换以获得所需的

M

样本输出。这也可以解决您的问题，因为您可以在感兴趣的区域中请求

M

样本，并且它不需要过多的内存，尽管它需要至少 3 倍的 CPU 时间。缺点是 CZT 的可靠实现还没有在 Numpy/Scipy 中：请参阅scipy issue 和它引用的code。 Matlab 的 CZT 似乎是可靠的，如果可以的话； Octave-forge 也有一个，Octave 的人通常会努力匹配/超越 Matlab。

但是如果你有记忆，零填充输入是要走的路。

Answer 2

Python 的 scipy 模块有一个 czt 函数。在您的情况下，您可以简单地使用 czt 作为快速 dft。但是它只是 1d。要执行 2d dft，我必须实现自己的方法。它基本上是沿着 x 和 y 顺序移动的。这是一个可运行的代码示例。

import numpy as np
import scipy
import matplotlib.pyplot as plt


def get_blob_pattern(size=100, sigma=30):
    # get an image of spherical aperture with smooth edge
    radius = 0.5 * size
    x_grid, y_grid = np.meshgrid(np.linspace(-radius, radius, size), np.linspace(-radius, radius, size))
    dist = (x_grid**2 + y_grid**2) ** 0.5
    img = 1 / (1 + np.exp(np.clip(3 * (dist - sigma), -20, 20)))
    return img


def get_dft2d(data, m, kmax):
    # perform 2d dft on image (using scipy 1d czt), get a 2d spectrum with frequency range of (-kmax, +kmax) on both axis
    # m: spectrum bin number
    w, h = data.shape
    aw = 2j * kmax / ((m - 1) * w)
    f1 = scipy.signal.czt(data, m=m, w=np.exp(aw), a=np.exp(0.5 * (m - 1) * aw), axis=0)
    data1 = f1 * np.exp(-0.5j * kmax * np.linspace(-1, 1, m))[:, np.newaxis]
    ah = 2j * kmax / ((m - 1) * h)
    f2 = scipy.signal.czt(data1, m=m, w=np.exp(ah), a=np.exp(0.5 * (m - 1) * ah), axis=1)
    data2 = f2 * np.exp(-0.5j * kmax * np.linspace(-1, 1, m))[np.newaxis, :]
    return data2


img = get_blob_pattern(size=100, sigma=20)
plt.imshow(img, vmin=0, vmax=1)
plt.show()

img_f = get_dft2d(img, 500, 100)
plt.imshow(np.abs(img_f) ** 2, vmin=0, vmax=10000)
plt.show()

我昨天遇到了这个确切的问题并且很难找到解决方案，最终不得不自己解决。遗憾的是 8 年没有答案。希望将来能对像我这样的物理系学生有所帮助。