NumPy FFT 产生偏离中心的输出

Question

TL;DR: NumPy FFT 在想要输出统一时创建非统一输出。我希望输出是均匀的电晕。

我正在尝试最终运行 Gerchberg-Saxton 相位检索算法。我一直在努力确保我了解 FFT 在 NumPy 中的工作原理。我使用fftshift 来创建外观正确的输出，但之后图像的强度不均匀。

我的输入图像是一个圆形，输出应该是一个从圆形光圈看日冕仪的东西。我正在尝试重现https://www.osapublishing.org/optica/fulltext.cfm?uri=optica-2-2-147&id=311836#articleSupplMat 中详述的结果。

我产生错误的算法：

1. 初始图像，f
2. FT(f)
3. x exp (i phase_mask)
4. IFT(FT(f)x exp(i phase_mask)

很高兴清理一切。

import numpy as np
import matplotlib.pyplot as plt
#Create 'pixels' for circle
pixels = 400
edge = np.linspace(-10, 10, pixels)
xv, yv = np.meshgrid(edge, edge)

def circle(x, y, r):
    '''
    x, y : dimensions of grid to place circle on
    r : radius
    Function defines aperture
    '''
    x0 = 0
    y0 = 0
    return np.select([((x-x0)**2+(y-y0)**2)>=r**2,
                      ((x-x0)**2+(y-y0)**2)<r**2],
                     [0,
                      1.])

#Create input and output images
radius = 4
input_img = circle(xv, yv, radius)
constraint_img = xcircle(xv, yv, radius)

img = input_img
constraint = 1 - img
max_iter = 10
re,im = np.mgrid[-1:1:400j, -1:1:400j] #Creates grid of values, 400=pixels
mask = 2*np.angle(re + 1j*im) #Gets angle from centre of grid
mask_i = mask


#Initial focal plane field, F. Initial image f.
f = np.sqrt(img)
F = np.fft.fftshift(np.fft.fft2(f)) * np.exp(mask * 1j)  #Focal plane field
F_1 = F 
am_f = np.abs(F_1) #Initial amplitude

g = np.fft.ifft2(F)
mask = np.angle(F/(F_1+1e-18))         #Final phase mask

recovery = (np.fft.ifft2(F*np.exp(-1j * mask)))
im3 = plt.imshow(np.abs(g)**2, cmap='gray')
plt.title('Recovered image')
plt.tight_layout()
plt.show()
plt.imshow(mask_i)
plt.colorbar()
plt.show()

Answer 1

您的问题在于这段代码：

pixels = 400
edge = np.linspace(-10, 10, pixels)

还有这个：

re,im = np.mgrid[-1:1:400j, -1:1:400j]

因为您使用fftshift*，所以您需要原点位于pixels//2。但是，您根本不对原点进行采样，它位于两个样本之间。

^{* 你应该真正使用ifftshift，它将原点从pixels//2 移动到0。fftshift 将原点从0 移动到pixels//2。对于偶数个样本，这两个做同样的事情。}

要正确采样原点，请按如下方式创建edge：

edge = np.linspace(-10, 10, pixels, endpoint=False)

我们现在看到edge[pixels//2] 等于 0。

对于np.mgrid，没有等效选项。您必须手动创建一个样本，然后删除最后一个样本：

re,im = np.mgrid[-1:1:401j, -1:1:401j] #Creates grid of values, 400=pixels
mask = 2*np.angle(re + 1j*im) #Gets angle from centre of grid
mask = mask[:-1, :-1]

通过这两个更改，您将看到对称输出。