Fresnel diffraction in two steps

Question

I wrote a short matlab script file that is suppose to run Fresnel's propagation (diffraction), such that given a certain input field U0, it will tell you how the field looks after distance z0. I compared the result to textbook results, and it seems like my program works fine. The problem is if I try to take two propagation steps instead of just one. i.e, instead of taking a single iteration of the program to propagate a distance z0, I take two iterations of the program to propagate a distance z0/2 each. Then I get complete nonsense, and I can't figure out what's the problem. Any advice will be accepted with great gratitude. Here is the code:

function U = fresnel_advance (U0, dx, dy, z, lambda)
% The function receives a field U0 at wavelength lambda
% and returns the field U after distance z, using the Fresnel
% approximation. dx, dy, are spatial resolution.

k=2*pi/lambda;
[ny, nx] = size(U0); 

Lx = dx * nx;
Ly = dy * ny;

dfx = 1./Lx;
dfy = 1./Ly;

u = ones(nx,1)*((1:nx)-nx/2)*dfx;    
v = ((1:ny)-ny/2)'*ones(1,ny)*dfy;   

O = fftshift(fft2(U0));

H = exp(1i*k*z).*exp(-1i*pi*lambda*z*(u.^2+v.^2));  

U = ifft2(O.*H);  

Answer 1

After calling fft2, you call also fftshift to have the DC frequency at the middle.

But when you call ifft2, the function assumes that you have still the DC frequency at (1,1). So you have to come back to this format before doing the inverse FFT.

So changing the final line to U = ifft2(fftshift(O.*H)) may solve the problem.

EDIT

I've just seen that Matlab advises to use ifftshift after fftshift insteaf of twice ifftshift (can't find the version in which it is introduced). According to the documentation, the sequences of calls ifftshift(fftshift(X)) and ifftshift(fftshift(X)) are not equivalent in case of odd sizes.

So I think it is better to do: U = ifft2(ifftshift(O.*H)) in the last step of your code.

Answer 2

It would be helpful if you could post some example inputs to your program to demonstrate the problem.

I suspect the problem may be related to the fact that you aren't calling FFTSHIFT enough times. Typically one considers the centre of the optical field matrix to be on the 'origin', where as FFT2 considers the 'bottom-left' corner. Therefore, you should FFTSHIFT before FFT2 as well as after.

You need to do the same thing for the IFFT2 piece too.

EDIT Justification for adding two calls to FFTSHIFT: compare these two:

N = 512; [x,y] = meshgrid(-1:1/N:(N-1)/N);
mask = (x.*x + y.*y) < 0.001;
figure(1)
imagesc(angle(fftshift(fft2(fftshift(mask)))))
figure(2)
imagesc(angle(fftshift(fft2(mask)))

Answer 3

In fact the problem lays in the way you run your fft. It is well explained in Fourier Optics and Computational Imaging by Khedar Kare, Wiley 2015:

The appropriate sequence for 2D FFT in most programming platforms for the result to be meaningful from Physical standpoint (e.g. in describing phenomena like diffraction) is thus given by: fftshift(fft2(ifftshift(...))).

In your code you should:O = fftshift(fft2(ifftshift(U0)));

Should you be interested in software developed in Python there is a rapidly developing Python toolbox for optics including diffraction: PyOptica. In PyOptica, a wavefront can be propagated using:

import astropy.units as u
import numpy as np

import pyoptica as po


wavelength = 500 * u.nm 
pixel_scale = 22 * u.um
npix = 1024

w = 6 * u.mm
h = 3 * u.mm
axis_unit = u.mm

wf = po.Wavefront(wavelength, pixel_scale, npix)
ap = po.RectangularAperture(w, h)
wf = wf * ap
fig_1 = po.plotting.plot_wavefront(wf, 'intensity', axis_unit=axis_unit)

The second step is to propagate:

f = 50 * u.cm
fs_f = po.FreeSpace(f)
wf_forward = wf * fs_f
fig_2 = po.plotting.plot_wavefront(wf_forward, 'intensity',  axis_unit=axis_unit)

It is important to remember sampling conditions for Fresnel propagation that is: ( z <= N (dx)^2 / lambda) where:

• N - number of pixels (in one direction);
• dx - pixel size;
• lambda - wavelength.

That condition is based on Computational Fourier Optics: A MATLAB Tutorial by David Voelz, SPIE 2011

You should implement that condition in your code. In PyOptica it is always checked before propagation; if the requested distance violates the condition, the propagation distance is broken up into substeps.

Answer 4

I think there's an error the phase term. "H = exp(1ikz).exp(-1ipilambdaz*(u.^2+v.^2)); "
should be "H = exp(1ikz).exp(-1ipi/lambda/z*(u.^2+v.^2)); "