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

Pleura

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.

Maximilien

Thank you very much! Sure enough, the problem was with the fftshift. Now it works – Pleura

Answer 2

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.

Wilonah

The programming language for this question is MATLAB; not Python – Autocorrelation

Answer 3

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

Gothurd

Not agree with your answer: a call of fftshift before fft2 is nonsense. Maybe you want to write ifft2? – Maximilien

Not sure I agree - I meant what I said, but maybe I am wrong. – Gothurd

What is the mean of fftshift before a fft? It swaps the spatial parts of your data. I don't think it is what you want to do. If you do do this also after ifft2, well, yes you will have the good result ... but the two "external" fftshift are useless. fftshift is intended to be used on frequential data, as the doc says it. – Maximilien

If you don't fftshift the input data, then you'll have a constant phase shift across the frequency plane. Added this to my answer since I can't add code here. – Gothurd

After a few articles read: you are totally right!! I don't have learnt such a thing during my school, completely missed that point. Thank you so much for pointing that! (and all my apologies for my stubborn comments). For readers: please read this article to understand precisely why you have to do fftshift(fft2(ifftshift(X))) to have a correct DFT with centered DC frequency. – Maximilien

No worries. Turns out my PhD in optics wasn't a waste of time then ;) – Gothurd

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

Mignon