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))