How do I generate a diffraction pattern using FFT from a single slit?

[py_engineer]

Hi,

I am a bit confused regarding using FFT to calculate the 1D single slit diffraction pattern.. Let's take the simple setup shown here:
http://electron9.phys.utk.edu/optics421/modules/m6/images/sslit.gif"

If I wanted to do a quick calculation using Matlab, I would create a 1xN matrix of zeros everywhere, and 1 in the middle wherever the slit is located. Then I can simply use the fft function from Matlab to get the diffraction pattern.

But what I really want to know is the diffraction pattern as a function of x (vertical x-axis) and determine the actual spot size of the central beam for example.

Let's say that I decide that the slit matrix at index 1 corresponds to x=-1 mm, slit matrix at index N corresponds to x=1 mm, and the slit is actually located between -0.1 mm and 0.1 mm (whatever indices these correspond to). Is there a correspondence in the x values for the fft of the slit? Does the value of the fft of the slit at index 1 corresponds to its value at x=-1 mm? Or how does it work?

Sorry for the very confusing message.. I hope someone can help.

Thanks!

 

[sophiecentaur]

The FFT is a discrete transform. You need to characterise your finite slit as a number of point sources, instead of a continuum.
Your model is of a single, infinitely narrow, slit and the pattern for that would be omnidirectional. As you increase the number of sources, the pattern will get nearer and nearer to the sinx/x function that you'd expect. Once your point sources are significantly less than a quarter wavelength apart, you get a picture which is more or less the same as the continuous slit apart from the far out sidelobes.

Edit: Of course, the FFT assumes a repeated pattern so you need a long array of zeros on either side of the row of ones.

 

[py_engineer]

Thanks for the reply. The problem is that I think I am not getting what I should get according to the theory.. Here is my Matlab code:

---
N=2^13;
lambda=525e-6; % Wavelength defined in mm
d_x0=2/N; % I go from -1 mm to 1 mm so d_x0 is the step between each points..
D=60;

for j=1:N
x0(j)=-1+(j-1)*d_x0;
if (x0(j)>=-2.5e-3&&x0(j)<=2.5e-3) % Slit location, total slit width is 50 microns
slit(j)=exp(j*2*pi/lambda*x0(j)^2/(2*D));
else slit(j)=0;
end;
end;

diffr=fft(slit);
diffr1=fftshift(diffr);

plot(x0,diffr1)
---

I feel that the pattern is not what I expect from theory.. The first zero should be for x0=lambda*D/(slit width).. but that is not what I am getting. Is it because N is too small?? Is there something clearly wrong with my code?

 

[sophiecentaur]

Tell me about this line:
"slit(j)=exp(j*2*pi/lambda*x0(j)^2/(2*D))"

Wouldn't the non zero values just be 1? Assuming the slit is uniformly illuminated. Sorry is that's a daft suggestion.
I think the scaling factor to get lambda to relate to the slit width should appear in the number of non zero values that you put into your fft. i.e. 50 microns' worth of ones.

 

[jfsalazars]

Hi,
The results of the FFT is the electric field, a complex magnitude, but what you see ( and of course what the photographic cameras register) is the intensity, (electric field)^2,
Try this program, is your program , buy I add some lines at the end
Hope this help
Jorge
---------------------------------------
N=2^13;
lambda=525e-6; % Wavelength defined in mm
d_x0=2/N; % I go from -1 mm to 1 mm so d_x0 is the step between each points..
D=60;

for j=1:N
x0(j)=-1+(j-1)*d_x0;
if (x0(j)>=-2.5e-3&&x0(j)<=2.5e-3) % Slit location, total slit width is 50 microns
slit(j)=exp(j*2*pi/lambda*x0(j)^2/(2*D));
else slit(j)=0;
end;
end;

diffr=fft(slit);
diffr1=fftshift(diffr);
%plot(x0,diffr1)
%%%
figure;
colormap 'gray'
imagesc(abs(diffr1));
---------------------------------------

 

[uart]

Tell me about this line:
"slit(j)=exp(j*2*pi/lambda*x0(j)^2/(2*D))"

Wouldn't the non zero values just be 1? Assuming the slit is uniformly illuminated. Sorry is that's a daft suggestion.

Not a daft suggestion at all sophiecentaur, I agree with it.
There's more than one problem with the OP's code but one glaring issue if that fact that he's using "j" as both a integer index variable and apparently as the sqrt(-1) complex variable. Matlab does pre-define both "i" and "j" as sqrt(-1), but if you use either of these variables for any other purpose then that pre-define is over-ridden.

The other issue is that the complex exponential in the OP's code is not needed, as the FFT itself takes care of that. Here's some extensively commented code I wrote to demonstrate how to generate the diffraction pattern, with correctly scaled power density, from the FFT of the aperture.

I=1000				% 1000 Watts per square meter = 1 Watt/mm (per meter of slit length).
N=2^12;				% Number of sample points.
L=0.66e-6			% Wavelength (red light).
D=1000e-6;			% Total x distance sampled (1mm). 
W=20e-6;			% Slit width (20um).
dx=D/N;				% The x-domain sample interval (approx).
dq=1/D;				% The q-domain sample interval  (q = theta/lambda).
dtheta=L*dq;    		% Angular displacement interval.
x=zeros(1,N);			% Construct x vector proportional to the field amplitude.
W2=round(W/2/dx);
x(N/2-W2 : N/2+W2)=sqrt(I);	% Set field strength constant over the slit width.
xf=dx*fftshift(fft(x));		% FFT scaled and shifted - Transform domain is q = theta/lambda).
x_domain_total_power=dx*sum(abs(x).^2)		% Compare the total power in x-domain with that
q_domain_total_power=dq*sum(abs(xf).^2)		% of the q-domain to confirm scaling is correct.
subplot(2,1,1)
plot([-N/2:N/2-1]*dx*1000,abs(x).^2/1000)	% Plot the x-domain intensity (converting to per mm).
xlabel('Slit Position (mm)')
ylabel('Watts/mm (per m length)')
title ('Slit power density versus position : Width=20um : Lambda=0.66um.')
subplot(212)
plot([-N/2 : N/2-1]*dtheta,abs(xf).^2/L)	% Plot q-domain intensity (converting to per radian).
xlabel('Angular Position (radians)')
ylabel('Watts/rad (per m length)')
title (Diffraction pattern, power density versus angular position.')

The above code generates the following graphs. (See attachment)