1
$\begingroup$

I have a question about my code which computes numerically the output field of a 4f setup with a pinhole in the middle which works as a spatial filter. My setup consists of two lenses with 50mm focal length (distance 2f), and a pinhole between both lenses. The input field has a Gaussian field distribution, where the Gaussian beam has a waist of 4mm and in addition I add to the input field some noise. My goal is to look how good I can filter out the noise for different pinhole diameters.

I wrote a program that simulates a spatial filter process in a 4f setup, and in general the spatial distribution of the output field has the shape that one would expect (Gaussian field goes in and Gaussian field comes out). The only problem that I have is that the field amplitudes of the output field are really high compared to the amplitudes of the input field (max values of fields 4 vs. 2000).

In my program I'm using the Fraunhofer approximation and solve the Fraunhofer diffraction integral (for a field starting in the front focal plane of a lens) numerically by using the FFT/IFFT functions from python. My code is inspired by the book: "Numerical simulation of optical wave propagation with examples in MATLAB" from Jason D. Schmidt (chapter 4). I guess there is something wrong with the scaling factors for the fft/ifft, but I have no idea where the error in my code is, since I'm using the exact same scaling factors like in the schmidt book.

(fft2 -> dl^2, ifft2 -> (N*dl_f)**2, where dl_f = 1/(N*dl))

The functions of my code are:

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 computational 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 computational 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

With these functions, I then calculate the spatial filtering process in the 4f setup as follows:

   # 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)
   

These are the plots of my results for a 0.1mm pinhole diameter:

Initial field

Filtered field in Fourier plane

Output field

(I can only give you these stack overflow links as I am not allowed to post images yet)

As you can see the filtering process works in my simulation, but the amplitudes of the output field are really high compared to the input field. Has someone an idea what I'm doing wrong? I'm looking for a solution of this problem for a view days now but without success. As I said before, I guess there is probably something wrong with my scaling factors for the fft.

I'm really grateful to anyone who wants to help me, because right now I have no idea what I'm doing wrong in my code.

$\endgroup$
  • $\begingroup$ Does the trivial case data=ifft(fft(data)) (without filtering) work? $\endgroup$ – ConvexHull Aug 29 at 19:59
  • $\begingroup$ You mean if I do the calculation without applying the filter function? In this case, the result is more or less the same. The field amplitudes are just a bit higher. $\endgroup$ – Nabla94 Aug 30 at 8:16
  • $\begingroup$ Input and output should be identical up to round off errors, also the amplitude. What factor do you get if you apply this to a constant function $factor=ifft(fft(c))/c$? $\endgroup$ – ConvexHull Aug 30 at 9:26
  • $\begingroup$ Yes, I agree, that's the behavior I would expect. So, if I understood you correctly, you mean that I should calculate the fft/ifft of a field value, and then divide the result by that initial field value, right? If I do that, I have a factor in the size of 609.66. $\endgroup$ – Nabla94 Aug 30 at 10:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.