[Numpy-discussion] Horizontal lines in diffraction image (NumPy FFT)
Matthias Hillenbrand
mailanhilli at googlemail.com
Wed Aug 6 17:30:59 EDT 2008
Hello,
> When you convolve two signals, of lengths N and M, you need to pad the
> FFTs to length (N+M-1) before multiplication.
> You can take a look at my linear position-invariant filtering code at:
> http://mentat.za.net/hg/filter
I understand your comments that I need zero padding when doing FFT
filtering. I will need some time to understand Stéfan's program as I
am a NumPy/Python newbie and not used to object oriented programming.
What I don't understand is that Octave/Matlab and Numpy generate
different results with nearly the same code. Why are the results with
Octave/Matlab OK although I have not applied additional zero padding.
In my opinion I should get exactly the same results.
Regards,
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