我有一个关于我的代码的问题，它以数字方式计算 4f 设置的输出字段，中间有一个针孔，用作空间滤波器。我的设置包括两个焦距为 50mm（距离 2f）的镜头，以及两个镜头之间的针孔。输入场具有高斯场分布，其中高斯光束的腰部为 4mm，此外我在输入场中添加了一些噪声。我的目标是看看我能过滤掉不同针孔直径的噪音有多好。
我编写了一个程序，在 4f 设置中模拟空间滤波过程，通常输出场的空间分布具有人们期望的形状（高斯场进入，高斯场出现）。我遇到的唯一问题是输出场的场幅度与输入场的幅度相比非常高（场 4 的最大值与 2000 相比）。
在我的程序中，我使用弗劳恩霍夫近似并通过使用 python 的 FFT/IFFT 函数数值求解弗劳恩霍夫衍射积分（对于从镜头前焦平面开始的场）。我的代码受到 Jason D. Schmidt（第 4 章）一书的启发：“用 MATLAB 中的示例对光波传播进行数值模拟”。我猜 fft/ifft 的比例因子有问题，但我不知道我的代码中的错误在哪里，因为我使用的比例因子与施密特书中的完全相同。
(fft2 -> dl^2, ifft2 -> (N*dl_f)**2，其中 dl_f = 1/(N*dl))
我的代码的功能是：
import numpy as np
from matplotlib import pyplot as plt
def lens_FT(U_in, lambda_0, L, N, focal):
'''Computes the field dirstribution of a field starting in the front focal
plane of a lens,and is then focused into the back focal plane of the
lens using Fraunhofer diffraction.
We assume a square shaped computational domain in this function.
Arguments
---------
U_in : 2d-array
Initial_field dirstibution
lambda_0 : float
Vacuum wavelength of the initial field in mm
focal : float
Focal length of the lens
L : float
length of computational domain at one side
N : int
Number of grid points for one side of the computaional domain
Returns
-------
U_out : 2d array
Field distribution in the Fourier plane of the initial field after
interaction with lens
kx : 1d array
Coordinates in the output plane in x-directiom
ky : 1d array
Coordinates in the output plane in y-directiom
'''
pass
# vacuum wavenumber
k_0 = 2*np.pi/lambda_0
# grid spacing of the initial field
dl = L/N
# define grid for output plane
kx = np.fft.fftfreq(N, d=dl)
ky = np.fft.fftfreq(N,d=dl)
# sort arrays
kx = np.asarray(np.append( kx[int(N/2)::], kx[0:int(N/2)] ))*lambda_0*focal
ky = np.asarray(np.append( ky[int(N/2)::], ky[0:int(N/2)] ))*lambda_0*focal
# define mesh for output field
kxx, kyy = np.meshgrid(kx, ky)
# shift U_in in order to perform the fft correct
U_in = np.fft.fftshift(U_in)
# compute output field with Fraunhofer integral
U_out = (np.exp(1j * k_0 * 2*focal) / (1j * lambda_0 * focal) *
np.fft.fftshift(np.fft.fft2(U_in)) * (dl)**2)
return U_out, kx ,ky
def lens_inv_FT(U_in, lambda_0, L, N, focal):
'''Computes the inverse Fourier tranform of a field which interacts with a
thin lens using Fraunhofer diffraction.
We assume a quadratic computational domain in this function.
Arguments
---------
U_in : 2d-array
Initial_field dirstibution
lambda_0 : float
Vacuum wavelength of the initial field in mm
focal : float
Focal length of the lens
L : float
length of computational domain at one side
N : int
Number of grid points for one side of the computaional domain
Returns
-------
U_out : 2d array
Field distribution in the Fourier plane of the initial field after
interaction with lens
x : 1d array
Coordinates in the output plane in x-directiom
y : 1d array
Coordinates in the output plane in y-directiom
'''
pass
# vacuum wavenumber
k_0 = 2*np.pi/lambda_0
# grid spacing of the initial field in spatial domain
dl = L/N
# grid spacing in spatial frequency domain
dl_f = 1/(N*dl)
# define mesh for output field
x = np.linspace(-L/2, L/2, N, endpoint=False)
y = np.linspace(-L/2, L/2, N, endpoint=False)
xx, yy = np.meshgrid(x, y)
U_in = np.fft.ifftshift(U_in)
U_out = (np.exp( 1j * k_0 *2*focal) / (1j * lambda_0 * focal) *
np.fft.ifftshift(np.fft.ifft2(U_in))* (N*dl_f)**2)
return U_out, x ,y
def pinhole_filter_fct(cx, cy, r, N ):
'''Creates filter function with pinhole shape
Arguments
---------
cx : int
matrix x index where pinhole is centered
cy : int
matrix y index where pinhole is centered
r : int
radius of the pinhole in mm
N : int
number of grid points of one side of the computational domain
Returns
-------
filter_fct : 2d-array
filter function with circular pinhole shape (matrix with 0 and 1 elements)
'''
# convert radius length into pixel
pixel_pinhole_rad = round(radius*N/L)
# Define matrix for filter fct
x = np.arange(-N/2, N/2)
y = np.arange(-N/2, N/2)
filter_fct = np.zeros((N, N))
# creating pinhole in
mask = (x[np.newaxis,:]-cx)**2 + (y[:,np.newaxis]-cy)**2 < pixel_pinhole_rad**2
filter_fct[mask] = 1
return filter_fct
使用这些函数，然后我计算 4f 设置中的空间滤波过程，如下所示：
# define parameter for simulation and physical parameter
# number of grid points per side
N = 602
# length of one side of the computational domain in mm
L = 15
# Gauss width in mm
w = 4
# radius of pinhole in mm
radius = 0.1
# vacuum wavelength of initial field in mm
lambda_0 = 0.00081
# focal length of the lens in mm
focal = 50
# define initial field (gaussian beam field distribution)
# define grids
x = np.linspace(-L/2, L/2, N)
y = np.linspace(-L/2, L/2, N)
xx, yy = np.meshgrid(x, y)
# build 2D gaussian field distribution array
gaussian_beam_in = 3*np.exp(-(xx**2 + yy**2)/(2*w**2))
# add noise to input field
noise = np.random.uniform(0,1,(N,N))
gaussian_beam_in = gaussian_beam_in + noise
# Compute field in Fourier domain
gaussian_beam_out, kx, ky = lens_FT(gaussian_beam_in, lambda_0, L, N, focal)
# compute filter function
pinhole_filter = pinhole_filter_fct(0, 0, radius, N)
# apply filter function by elementwise matrix multiplication with field in Fourier plane
gaussian_beam_out_filtered = np.multiply(gaussian_beam_out, pinhole_filter)
# compute field in output plane at 4f
gaussian_beam_out_4f, x, y = lens_inv_FT(gaussian_beam_out_filtered, lambda_0, L, N, focal)
这些是我对 0.1 毫米针孔直径的结果图：
（我只能给你这些堆栈溢出链接，因为我还不允许发布图像）
如您所见，过滤过程在我的模拟中有效，但输出场的幅度与输入场相比非常高。有人知道我做错了什么吗？我现在正在寻找这个问题的解决方案，但没有成功。正如我之前所说，我猜我的 fft 比例因子可能有问题。
我真的很感谢任何想帮助我的人，因为现在我不知道我的代码做错了什么。
