python - 模拟空间滤波器 4f 设置的问题（Python）

问题描述

我有一个关于我的代码的问题，它以数字方式计算 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 比例因子可能有问题。

我真的很感谢任何想帮助我的人，因为现在我不知道我的代码做错了什么。

标签： pythonfftsimulationspatialoptics-algorithm

python - 模拟空间滤波器 4f 设置的问题（Python）

问题描述

解决方案

推荐阅读