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 <[EMAIL PROTECTED]>: > > 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 > > > _______________________________________________ > Numpy-discussion mailing list > [email protected] > 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 _______________________________________________ Numpy-discussion mailing list [email protected] http://projects.scipy.org/mailman/listinfo/numpy-discussion
