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

8 次查看（过去 30 天）
Sohail Ahmed
Sohail Ahmed 2024-6-13
评论： Sohail Ahmed 2024-10-1
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 个评论
显示 1更早的评论
David Goodmanson
David Goodmanson 2024-6-17
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 2024-6-18
Yeah, there is a reason. The issue arises with the zero padding of the Fresnel inetgral while applying fft2.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Sandeep Mishra
Sandeep Mishra 2024-9-17
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.

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by