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. -----הודעה מקורית----- מאת: [EMAIL PROTECTED] בשם Matthias Hillenbrand נשלח: ד 06-אוגוסט-08 05:26 אל: [email protected] נושא: [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 [email protected] http://projects.scipy.org/mailman/listinfo/numpy-discussion
<<winmail.dat>>
_______________________________________________ Numpy-discussion mailing list [email protected] http://projects.scipy.org/mailman/listinfo/numpy-discussion
