[Numpy-discussion] Horizontal lines in diffraction image (NumPy FFT)
Nadav Horesh
nadavh at visionsense.com
Wed Aug 6 05:09:07 EDT 2008
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 at scipy.org בשם Matthias Hillenbrand
נשלח: ד 06-אוגוסט-08 05:26
אל: numpy-discussion at 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 at scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
-------------- next part --------------
A non-text attachment was scrubbed...
Name: winmail.dat
Type: application/ms-tnef
Size: 5713 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080806/1c9e95f4/attachment.bin>
More information about the NumPy-Discussion
mailing list