[Numpy-discussion] Horizontal lines in diffraction image (NumPy FFT)

Matthieu Brucher matthieu.brucher@gmail....
Wed Aug 6 04:32:20 CDT 2008


Exactly. Using FFT to do a convolution should be done after some
signal processing readings ;) (That's why I hate FFT to do signal
processing as well).

Matthieu

2008/8/6 Nadav Horesh <nadavh@visionsense.com>:
>
> I did about the same thing 9 year ago (in python of course). If I can recall correctly, you need to double the arrays size (with padding of 0) in order to avoid this artifact. I think that its origin is that the DFT is equivalent to periodic boundary conditions.
>
>  Nadav.
>
> -----הודעה מקורית-----
> מאת: numpy-discussion-bounces@scipy.org בשם Matthias Hillenbrand
> נשלח: ד 06-אוגוסט-08 05:26
> אל: numpy-discussion@scipy.org
> נושא: [Numpy-discussion] Horizontal lines in diffraction image (NumPy FFT)
>
> Hello,
>
> I am trying to calculate the propagation of a Gaussian beam with an
> angular spectrum propagation algorithm.
> My program does the following steps:
> 1. Generation and multiplication of the Gaussian beam, an aperture,
> and a lens -> u
> 2. FFT(u)
> 3. Multiplication of FFT(u) with the angular spectrum kernel H
> 4. IFFT(FU*H)
>
> Steps 2-4 are repeated 1000 times to show the propagation of the beam
> in z-direction
>
> Unfortunately the calculation results in a lot of horizontal lines
> that should not be there. The program produces reasonable results
> besides this problem.
>
> It is especially interesting that an equivalent calculation with
> Matlab or Octave has the same results without these artifacts. Both
> calculations take approximately 90 seconds.
>
> I am not sure whether I am doing something wrong or there is a
> difference in the FFT codes. I assume that the problem has something
> to do with the propagation kernel H but I cannot see a difference
> between both programs.
>
> Thank you very much for your help!
>
> Matthias
>
>
> ---------------------------------------------------------------------------------------------------
> Python code with comments (version without comments below)
> ---------------------------------------------------------------------------------------------------
>
> from pylab import *
> from numpy.lib import scimath  #Complex roots
> import os
>
> #########################################
> ##########  Defintion of the input variables  #########
> #########################################
>
>
> ##General variables
> lam     = 532.e-9                   #wavelength in m
> k       = 2*pi/lam                  #wave number
>
> ##Computational array
> arr_size    = 80.e-3                #size of the computation array in m
> N       = 2**17                     #points of the 1D array
>
> ##Aperture
> ap      = 40.e-3                    #aperture size of the 1D array in m
>
> ##Gaussian beam
> w0      =10.e-3                     #waist radius of the gaussian beam in m
>
> ##Lenses (definition of the focal length)
> n       = 1.56                      #refractive index of the glass
> f       = .1                        #focal length in m
> Phi     = 1/f                       #focal power of the lens
> r2      = 1.e18                     #radius of the second lens surface
> r1      = 1 / ( Phi/(n-1) + 1/r2 )  #computation of the radius of the
> first lens surface
>
> ##z distances
> zmin    = 0.99*f                    #initial z position
> zmax    = 1.01*f                    #final z position
> zsteps  = 1000                      #number of computated positions
> ROI     = 1000                      #Region of interest in the diffraction image
>
> ##############################################
> ###############  Program execution 1D  ################
> ###############################################
>
> x       = linspace(-arr_size/2,arr_size/2,N)
> A       = ones(N)
> A[where(ap/2<=abs(x))]  = 0         #Definition of the aperture
> G       = exp(- x**2/w0**2)         #Generation of the Gaussian beam
>
> delta = -r1*(1 - scimath.sqrt(1 - x**2/r1**2))
> +  r2*(1 - scimath.sqrt(1 - x**2/r2**2))
> m           = (r1**2 <= x**2 )
> delta[m]    = 0                     #correction in case of negative roots
> m           = (r2**2 <= x**2 )
> delta[m]    = 0                     #correction in case of negative roots
> lens        = exp(1j * 2 * pi / lam * (n-1) * delta)    #Calculation
> of the lens phase function
>
> u       = A*G*lens                  #Complex amplitude before beam propagation
>
> ############################################
> ###########  Start angular spectrum method  ###########
>
> deltaf  = 1/arr_size                            #spacing in frequency
> domain
> fx      = r_[-N/2*deltaf:(N/2)*deltaf:deltaf]   #whole frequency domain
> FU      = fftshift(fft(u))                      #1D FFT
> U       = zeros((zsteps,ROI))                   #Generation of the image array
> z       = linspace(zmin,zmax,zsteps)            #Calculation of the
> different z positions
>
> for i in range(zsteps):
>    H = exp(1j * 2 * pi / lam * z[i] * scimath.sqrt(1-(lam*fx)**2))
>    U[i]   = (ifft(fftshift(FU*H)))[N/2 - ROI/2 : N/2 + ROI/2]
>    if i%10 == 0:
>        t = 'z position: %4i'     % i
>        print t
>
> ############  End angular spectrum method  ############
> #############################################
>
> res         = abs(U)
>
> imshow(res)
> show()
>
>
>
>
>
>
> ----------------------------------------------------------
> Python code without comments
> ----------------------------------------------------------
>
> from pylab import *
> from numpy.lib import scimath  #Complex roots
> import os
>
> lam     = 532.e-9
> k       = 2*pi/lam
>
> arr_size    = 80.e-3
> N       = 2**17
>
> ap      = 40.e-3
>
> w0      =10.e-3
>
> n       = 1.56
> f       = .1
> Phi     = 1/f
> r2      = 1.e18
> r1      = 1 / ( Phi/(n-1) + 1/r2 )
>
> zmin    = 0.99*f
> zmax    = 1.01*f
> zsteps  = 1000
> ROI     = 1000
>
> x       = linspace(-arr_size/2,arr_size/2,N)
> A       = ones(N)
> A[where(ap/2<=abs(x))]  = 0
> G       = exp(- x**2/w0**2)
>
> delta = -r1*(1 - scimath.sqrt(1 - x**2/r1**2))
> +  r2*(1 - scimath.sqrt(1 - x**2/r2**2))
> m           = (r1**2 <= x**2 )
> delta[m]    = 0
> m           = (r2**2 <= x**2 )
> delta[m]    = 0
> lens        = exp(1j * 2 * pi / lam * (n-1) * delta)
>
> u       = A*G*lens
>
> deltaf  = 1/arr_size
> fx      = r_[-N/2*deltaf:(N/2)*deltaf:deltaf]
> FU      = fftshift(fft(u))
> U       = zeros((zsteps,ROI))
> z       = linspace(zmin,zmax,zsteps)
>
> for i in range(zsteps):
>    H = exp(1j * 2 * pi / lam * z[i] * scimath.sqrt(1-(lam*fx)**2))
>    U[i]   = (ifft(fftshift(FU*H)))[N/2 - ROI/2 : N/2 + ROI/2]
>    if i%10 == 0:
>        t = 'z position: %4i'     % i
>        print t
>
> res         = abs(U)
>
> imshow(res)
> show()
>
>
>
>
>
>
> ----------------------------------------------------------
> Matlab, Octave code
> ----------------------------------------------------------
>
> lam     = 532.e-9
> k       = 2*pi/lam
>
> arr_size    = 80.e-3
> N       = 2^17
>
> ap      = 40.e-3
>
> w0      =10.e-3
>
> n       = 1.56
> f       = .1
> Phi     = 1/f
> r2      = 1.e18
> r1      = 1 / ( Phi/(n-1) + 1/r2 )
>
> zmin    = 0.99*f
> zmax    = 1.01*f
> zsteps  = 1000
> ROI     = 1000
>
> x=linspace(-arr_size/2,arr_size/2,N);
> A       = ones(1,N);
> A(find(ap/2<=abs(x)))=0;
> G       = exp(- x.^2/w0^2);
>
> delta = -r1*(1 - sqrt(1 - x.^2/r1^2))  +  r2*(1 - sqrt(1 - x.^2/r2^2));
> delta(find(r1.^2 <= x.^2 ))=0;
> delta(find(r2.^2 <= x.^2 ))=0;
> lens = exp(1j * 2 * pi / lam * (n-1) * delta);
> u       = A.*G.*lens;
>
> deltaf  = 1/arr_size;
> fx      = [-N/2*deltaf:deltaf:(N/2-1)*deltaf];
> FU      = fftshift(fft(u));
> U       = zeros(zsteps,ROI);
> z       = linspace(zmin,zmax,zsteps);
>
> for i=1:zsteps
> H    = exp(1j * 2 * pi / lam * z(i) * sqrt(1-lam.^2*fx.^2));
> U(i,:) = ifft(fftshift(FU.*H))(N/2 - ROI/2 : N/2 + ROI/2-1);
> end
>
> imagesc(abs(U))
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
>



-- 
French PhD student
Website : http://matthieu-brucher.developpez.com/
Blogs : http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn : http://www.linkedin.com/in/matthieubrucher


More information about the Numpy-discussion mailing list