1. Foreword
In Fresnel Propagation 1 (Fresnel Propagation) and Fresnel Propagation 2 (Fresnel Propagation), we talked about the integral formula of Fresnel propagation, and introduced a method of simulating Fresnel propagation according to the integral formula. For this simulation method, we also verify the correctness of our simulation through the analytical solution of Fresnel propagation. In this article, we introduce a second method for simulating Fresnel propagation, the angular spectrum method.
2. Angular Spectral Method
2.1 Fresnel integral formula
In Fresnel Propagation 1 (Fresnel Propagation), we introduced the Fresnel integral formula, here we write it again,
u
2
(
x
2
,
the y
2
)
=
e
j
k
z
j
lambda
z
∫
?
∞
∞
∫
?
∞
∞
u
1
(
x
1
,
the y
1
)
e
j
k
2
z
[
(
x
2
?
x
1
)
2
+
(
the y
2
?
the y
1
)
2
]
d
x
1
d
the y
1
,
(1)
U_2\left(x_{2}, y_{2}\right) =\frac{e^{j k z}}{j \lambda z} \int_{-\infty}^{\infty} \int_{-\infty }^{\infty} U_1\left(x_{1}, y_{1}\right) e^{j \frac{k}{2 z}\left[\left(x_{2}-x_{1} \right)^{2} + \left(y_{2}-y_{1}\right)^{2}\right]} d x_{1} d y_{1},\tag{1}
U2?(x2?,y2?)=jλzejkz?∫?∞∞?∫?∞∞?U1?(x1?,y1?)ej2zk?[(x2x1?)2 + (y2y1?) 2] dx1?dy1?,(1)
In formula (1), if we look at this formula in the form of convolution, then the impulse response function is
h
(
x
1
,
the y
1
)
=
e
j
k
z
j
lambda
z
exp
?
[
j
k
2
z
(
x
1
2
+
the y
1
2
)
]
,
(2)
h\left ( x_1,y_1 \right ) =\frac{e^{jkz}}{j\lambda z} \exp \left [ \frac{jk}{2z}\left ( x_1^{2} + y_1^ {2}\right ) \right ] ,\tag{2}
h(x1?,y1?)=jλzejkz?exp[2zjk?(x12? + y12?)],(2)
With this impulse response function, we can further write the Fresnel integral in formula (1) as the following convolution form,
u
2
(
x
2
,
the y
2
)
=
?
u
1
(
x
1
,
the y
1
)
h
(
x
2
?
x
1
,
the y
2
?
the y
1
)
d
x
1
d
the y
1
,
(3)
U_2\left ( x_2,y_2 \right ) =\iint U_1\left ( x_1,y_1 \right )h\left ( x_2 – x_1, y_2 – y_1\right )dx_1dy_1 ,\tag{3}
U2?(x2?,y2?)=?U1?(x1?,y1?)h(x2x1?,y2y1?)dx1?dy1?,(3)
In convolution theory, the convolution of two functions can be obtained by the product of their respective Fourier transforms, then formula (3) can be further written as
u
2
(
x
,
the y
)
=
f
?
1
{
f
{
u
1
(
x
,
the y
)
}
f
{
h
(
x
,
the y
)
}
}
,
(4)
U_2(x,y)=F^{-1}\left \{ F\left \{ U_1(x,y) \right \} F\left \{ h(x,y) \right \} \right \ } ,\tag{4}
U2?(x,y)=F?1{F{U1?(x,y)}F{h(x,y )}},(4)
In formula (4),
f
f
F stands for Fourier transform,
f
?
1
F^{-1}
F?1 stands for inverse Fourier transform.
The expression of the impulse function (that is, formula (2)) in the Fourier frequency domain, we can also directly give
h
(
f
x
,
f
the y
)
=
e
j
k
z
exp
?
[
?
j
π
lambda
z
(
f
x
2
+
f
the y
2
)
]
,
(5)
H\left ( f_x,f_y \right ) =e^{jkz}\exp\left [ -j\pi \lambda z\left ( f_{x}^{2} + f_{y}^{2}\right ) \right ] ,\tag{5}
H(fx?,fy?)=ejkzexp[?jπλz(fx2? + fy2?)],(5)
Equation (5) is often called the transfer function. Thus, we have two simulation methods of the angular spectrum, one based on formula (2) and the other based on formula (5).
2.2 Fresnel angle spectral simulation of square hole based on formula (2)
First, we define the following function according to formula (2):
function [U2] = angular_spectrum_fresnel_prop_impulse(U1, lambda, delta1,z) % U1: source plane % lambda: wavelength % delta1: sampling interval in the source plane, delta_x1 and delta_y1, % delta_x1 = delta_y1 % z:distance of the central point between source plane and observation % plane, i.e., the propagation distance % U2: output in the observation plane [M, N] = size(U1); % obtain the number of points along each side k = 2 * pi / lambda; % wave number % source plane coordinates [x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-M/2 : 1: M/2 -1) * delta1); % impulse function h = exp(1j * k * z) / (1j * lambda * z) * exp((1j*k)/(2*z) * (x1.^2 + y1.^2)); % Fourier transform of U1 and h F_U1 = fftshift(fft2(fftshift(U1))); H = fftshift(fft2(fftshift(h))) * delta1^2; % product between F_U1 and H F_U1_H = F_U1 .* H; % Fresnel propagation integral U2 = fftshift(ifft2(fftshift(F_U1_H))); end
Square hole Fresnel diffraction based on formula (2)
clear;close all;clc
%% parameters
a = 1e-3; % width of square shape unit: m
N = 1024; % number of points along each side (row and column)
L = 1e-2; % physical size of the source plane, unit: m
delta1 = L / N; % sampling interval in the source plane
lambda = 532e-9; % wavelength,unti: m
z = 0.1; % propagation distance, unit: m
%% square shape
% cooridates in the source plane
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-N/2 : 1: N/2 -1) * delta1);
rect_shape = rect(x1, a).* rect(y1, a);
figure(1)
imagesc(x1(1,:),y1(:,1),rect_shape); %image display
colormap('gray'); %linear gray display map
axis square; % square figure
axis xy % y positive up
xlim([-0.005,0.005]) % x axis range
ylim([-0.005,0.005]) % y axis range
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);
%% Fresnel propagation of square shape
[Fres_rect] = angular_spectrum_fresnel_prop_impulse(rect_shape, lambda, delta1, z);
%%
figure(2),
imagesc(x1(1,:),y1(:,1),abs(Fres_rect)); %image display
colormap('gray'); %linear gray display map
axis square; % square figure
axis xy % y positive up
xlim([-0.005,0.005]) % x axis range
ylim([-0.005,0.005]) % y axis range
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);
The result of running the above code is shown in Figure 1.
(a) square holes;
(b) Fresnel diffraction spot based on formula (2).
Figure 1. (a) Square hole; (b) Fresnel diffraction spot based on formula (2).
———————————————
The comparison code with the analytical solution is as follows
U2 = analytical_square_fresnel_propagation(x1,y1,a,lambda,z);
%%
figure(3),
imagesc(x1(1,:),y1(:,1),abs(U2)); %image display
colormap('gray'); %linear gray display map
axis square; % square figure
axis xy % y positive up
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);
%% compare the 1D result
% compare the simulation result with the analytical result based on the 513 row
sim = abs(Fres_rect(513,:));
ana = abs(U2(513,:));
figure(4),
plot(x1(1,:),sim,'--r'), hold on
plot(x1(1,:),ana,'b')
% xlim([min(x1(1,:)),max(x1(1,:))]);
xlim([-0.005,0.005]) % x axis range
ylim([0,1.5]);
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('Magnitude', 'FontSize', 16);
legend('Simulation Result','Analytical Result')
The result of running the above code is as follows
(a) analytical solution;
(b) The comparison result of simulated solution and analytical solution at line 513;
Figure 2. (a) Analytical solution; (b) Comparison result of simulated solution and analytical solution at line 513;
2.3 Fresnel angle spectral simulation of square hole based on formula (5)
First, we define the following function according to formula (5):
function [U2] = angular_spectrum_fresnel_prop_transfer(U1, lambda, delta1,z) % U1: source plane % lambda: wavelength % delta1: sampling interval in the source plane, delta_x1 and delta_y1, % delta_x1 = delta_y1 % z:distance of the central point between source plane and observation % plane, i.e., the propagation distance % U2: output in the observation plane [M, N] = size(U1); % obtain the number of points along each side k = 2 * pi / lambda; % wave number % source plane coordinates [fx, fy] = meshgrid((-N/2 : 1: N/2 -1) * (1 / (N *delta1)), (-M/2 : 1: M/2 -1) * (1 / (M *delta1))); % transfer function H = exp(1j * k * z) * exp(-1j * pi * lambda * z * (fx.^2 + fy.^2)); % Fourier transform of U1 F_U1 = fftshift(fft2(fftshift(U1))); % product between F_U1 and H F_U1_H = F_U1 .* H; % Fresnel propagation integral U2 = fftshift(ifft2(fftshift(F_U1_H))); end
Square hole Fresnel diffraction based on formula (5)
clear;close all;clc
%% parameters
a = 1e-3; % width of square shape unit: m
N = 1024; % number of points along each side (row and column)
L = 1e-2; % physical size of the source plane, unit: m
delta1 = L / N; % sampling interval in the source plane
lambda = 532e-9; % wavelength,unti: m
z = 0.1; % propagation distance, unit: m
%% square shape
% cooridates in the source plane
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-N/2 : 1: N/2 -1) * delta1);
rect_shape = rect(x1, a).* rect(y1, a);
figure(1)
imagesc(x1(1,:),y1(:,1),rect_shape); %image display
colormap('gray'); %linear gray display map
axis square; % square figure
axis xy % y positive up
xlim([-0.005,0.005]) % x axis range
ylim([-0.005,0.005]) % y axis range
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);
%% Fresnel propagation of square shape
[Fres_rect] = angular_spectrum_fresnel_prop_transfer(rect_shape, lambda, delta1, z);
%%
figure(2),
imagesc(x1(1,:),y1(:,1),abs(Fres_rect)); %image display
colormap('gray'); %linear gray display map
axis square; % square figure
axis xy % y positive up
xlim([-0.005,0.005]) % x axis range
ylim([-0.005,0.005]) % y axis range
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);
The result of running the above code is as follows
Fig. 3. Simulation results of Fresnel diffraction with square hole based on transfer function.
Comparison with Analytical Solution
Figure 4. The comparison results of the Fresnel diffraction simulation results of the square hole based on the transfer function and the analytical solution at line 513;
3. Conclusion
In this paper, we investigate the angular spectroscopic simulation method of Fresnel diffraction. In this method, we implemented the method in two ways, the impulse function way and the transfer function way. These two methods are mathematically equivalent, but in simulation, due to discrete sampling, the simulation results will be slightly different. From the results, both methods can simulate Fresnel propagation well. Under the given parameter conditions in this paper, the simulation results corresponding to the transfer function have a better matching degree.