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

Matthias Hillenbrand mailanhilli@googlemail....
Tue Aug 5 21:26:08 CDT 2008

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

More information about the Numpy-discussion mailing list