[Numpy-discussion] Horizontal lines in diffraction image (NumPyFFT)

Nadav Horesh nadavh@visionsense....
Thu Aug 7 00:34:06 CDT 2008


The 0 padding is easy in numpp/pylab as in octave/matlab, by just adding one argument. In pylab it is the "a" keyword:

 y = fft(x, n=2*len(x))

if n is not given then n=len(x) --- normal fft
if n > len(x) then it pads x with 0.

 Nadav.


-----הודעה מקורית-----
מאת: numpy-discussion-bounces@scipy.org בשם Matthias Hillenbrand
נשלח: ה 07-אוגוסט-08 00:30
אל: numpy-discussion@scipy.org
נושא: Re: [Numpy-discussion] Horizontal lines in diffraction image (NumPyFFT)
 
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))
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 5644 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/numpy-discussion/attachments/20080807/d1c828fc/attachment-0001.bin 


More information about the Numpy-discussion mailing list