How to compute diffraction integral using fast Fourier transform (fft2) ?

4 views (last 30 days)
Hi there !
The following code computes Fresnel diffraction integral using fft2. It works for (obj.E = ones(Ny, Nx);) and performs the fft2 analysis, but when the input is (obj.E = zeroes(Ny, Nx);), it plots nothing. I checked the same code on python with zeros and it worked but in MATLAB it is not performing fft2 for zeros (gives only zero values). Kindly help me resolve this issue.
The Sheet.m (code):
classdef Sheet
properties
dx
dy
x
y
xx
yy
Nx
Ny
E
end
methods
function obj = Sheet(extent_x, extent_y, Nx, Ny)
obj.dx = extent_x / Nx;
obj.dy = extent_y / Ny;
obj.x = obj.dx * ((0:Nx-1) - floor(Nx/2));
obj.y = obj.dy * ((0:Ny-1) - floor(Ny/2));
[obj.xx, obj.yy] = meshgrid(obj.x, obj.y);
obj.Nx = Nx;
obj.Ny = Ny;
obj.E = zeros(Ny, Nx);
end
function rectangular_slit(obj, x0, y0, lx, ly)
% Creates a slit centered at the point (x0, y0) with width lx and height ly
mask = ((obj.xx > (x0 - lx/2)) & (obj.xx < (x0 + lx/2))) & ((obj.yy > (y0 - ly/2)) & (obj.yy < (y0 + ly/2)));
obj.E = obj.E + mask;
end
end
end
The main code :
% Simulation input
Lx = 1.4;
Ly = 0.4;
Nx = 2500;
Ny = 1500;
sheet = Sheet(2*Lx, 2*Ly, Nx, Ny);
% Slit separation
mm = 1e-3;
D = 128 * mm;
sheet.rectangular_slit(-D/2, 0, 22 * mm, 88 * mm);
sheet.rectangular_slit(D/2, 0, 22 * mm, 88 * mm);
% Distance from slit to the screen (mm)
z = 5000;
% Wavelength (mm)
lambda = 18.5*1e-7;
k = 2*pi / lambda;
% Initialize complex array for phase information
fft_c = fft2(sheet.E .* exp(1i * k/(2*z) * (sheet.xx.^2 + sheet.yy.^2)));
c = fftshift(fft_c);
% Plot with MATLAB
abs_c = abs(c);
% Screen size (mm)
dx_screen = z*lambda/(2*Lx);
dy_screen = z*lambda/(2*Ly);
x_screen = dx_screen * ((1:Nx) - Nx/2);
y_screen = dy_screen * ((1:Ny) - Ny/2);
figure;
subplot(2, 1, 1);
imagesc(x_screen([1 end]), y_screen([1 end]), abs_c);
colormap(gray);
axis equal;
xlabel('x (mm)');
ylabel('y (mm)');
title('Probability Density |\psi|^2');
xlim([-2, 2]);
ylim([-1, 1]);
subplot(2, 1, 2);
plot(x_screen, abs(c(Ny/2, :)).^2, 'LineWidth', 1);
xlabel('x (mm)');
ylabel('Probability Density |\psi|^2');
title('Probability Density |\psi|^2 along y=0');
xlim([-2, 2]);
  3 Comments
David Goodmanson
David Goodmanson on 17 Jun 2024
Hi SK. if the input is all zeros, is there a reason that you expect that the output would not be all zeros? If anyting I might be questioning the python code that gives nonzero output for zero input.
Sohail Ahmed
Sohail Ahmed on 18 Jun 2024
Yeah, there is a reason. The issue arises with the zero padding of the Fresnel inetgral while applying fft2.

Sign in to comment.

Answers (1)

Sandeep Mishra
Sandeep Mishra on 17 Sep 2024
Hi Sohail,
Upon executing the code snippet in MATLAB R2024a, I encountered the same outcome. During the debugging process, I identified that the 'rectangular_slit' function within the 'Sheet' class is not updating the 'E' variable as expected.
To resolve this issue, you should modify the 'rectangular_slit' function to return the modified object. Additionally, ensure that the 'Sheet' object is updated accordingly after invoking the 'rectangular_slit' function in your main script.
Please update the "rectangular_slit" method of "Sheet" class as follows:
% Updated function for sheet.m file
function obj = rectangular_slit(obj, x0, y0, lx, ly)
% Creates a slit centered at the point (x0, y0) with width lx and height ly
mask = ((obj.xx > (x0 - lx/2)) & (obj.xx < (x0 + lx/2))) & ((obj.yy > (y0 - ly/2)) & (obj.yy < (y0 + ly/2)));
obj.E = obj.E + mask;
end
Please update the call to "rectangular_slit" in the main code as follows:
% Updation function call for main file
% Add two slits to the sheet
sheet = sheet.rectangular_slit(-D/2, 0, 22 * mm, 88 * mm);
sheet = sheet.rectangular_slit(D/2, 0, 22 * mm, 88 * mm);
For more information, refer to the following MathWorks documentation:
  1. ‘Method Syntax’: https://www.mathworks.com/help/releases/R2024a/matlab/matlab_oop/specifying-methods-and-functions.html
  2. Using a Class to Display Graphics’: https://www.mathworks.com/help/releases/R2024a/matlab/matlab_oop/a-class-code-listing.html#buq2wcy
I hope this helps.

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!