# Code of the Angular Spectrum Method

• ecastro
In summary: either by trying to use spatial frequency coordinates (maybe someone can help with that) or by rewriting the code so that spatial frequency coordinates are not needed.
ecastro
I constructed my code of the Angular Spectrum Method. However, as the distance between the object and the plane of interest increases, the diffraction pattern never disappears; there is still some sort of a diffraction pattern, and I am expecting that it disappears as distance increases.

Here is the code:

Matlab:
clear; clc;

layer = zeros(499, 499);
diameter = 100*10^(-6);
wavelength = 500*10^(-9);
altitude = 2*10^(-3);
img_width = 1*10^(-3);
img_length = 1*10^(-3);

lambda = wavelength;                          % wavelength (meters)
z = altitude;                                 % altitude (meters)
k = 2*pi/lambda;                              % wavenumber
phy_x = img_width;                            % physical width (meters)
phy_y = img_length;                           % physical length (meters)
obj_size = size(layer);

% alpha and beta (wavenumber components)
dx = linspace(-phy_x/2, phy_x/2, obj_size(2));
dx = dx(1:length(dx));
dy = linspace(-phy_y/2, phy_y/2, obj_size(1));
dy = dy(1:length(dy));

for i = 1:obj_size(2);
for j = 1:obj_size(1)
if sqrt(dx(i)^2 + dy(j)^2) <= radius;
layer(j, i) = 1;
end;
end;
end;

U0 = fftshift(fft2(layer));

Fs_x = obj_size(2)/img_width;
Fs_y = obj_size(1)/img_length;

dx2 = Fs_x^(-1);
dy2 = Fs_y^(-1);

x2 = dx2*(0:(obj_size(2) - 1))';
y2 = dy2*(0:(obj_size(1) - 1))';

dFx = Fs_x/obj_size(2);
dFy = Fs_y/obj_size(1);

Fx = (-Fs_x/2:dFx:(Fs_x/2 - dFx));
Fy = (-Fs_y/2:dFy:(Fs_y/2 - dFy));

alpha = lambda.*Fx; beta = lambda.*Fy;

% gamma
gamma = zeros(length(beta), length(alpha));
for j = 1:length(beta);
for i = 1:length(alpha);
if (alpha(i)^2 + beta(j)^2) > 1;
gamma(j, i) = 0;
else
gamma(j, i) = sqrt(1 - alpha(i)^2 - beta(j)^2);
end;
end;
end;

U1 = ifft2(fftshift(fft2(layer)).*exp(1i*k.*gamma.*z));
I1 = (1/(16*pi)).*(U1.*conj(U1));

I think my error came from the calculation of the alpha and beta components. I just don't know where...

I adjusted your post as it looked like MATLAB code and I split some lines where you did assignements like somex = some_expression; somey=some_expression;

Have you tried stepping through the computation to see where things go astray?

Everything seems to be working fine at very small distances, where diffraction is expected to occur. However, as I increase the distance ##z##, the diffraction never disappears.

Aside: gamma is a built-in MATLAB function, so you probably shouldn't name your variable with the same name as it will interfere with that functionality for you and anyone that runs your program.

The variable z only shows up in one place after it is defined: as an exponent to exp in the calculation of U1. I assume this U1 is the diffraction pattern that you want? (You might want to comment your code better so that it's clear what all the variables are.)

The only other quantity of exp that could change is gamma, which depends on alpha and beta. Upon looking at gamma, I realized that your section of code that sets elements of gamma to zero is "dead" (seemingly not used):

Code:
if (alpha(i)^2 + beta(j)^2) > 1;
gamma(j, i) = 0;

When I ran your program, gamma has all nonzero values:

Code:
nnz(gamma)/numel(gamma)

ans =

1

So there might be something up there, if you were expecting some elements of gamma to be zero. But ultimately, if gamma is unchanged and you vary z, then the exp term gets bigger and bigger. I would expect the exponent to be negative here if you wanted the exp term to disappear as z increases.

kreil said:
Aside: gamma is a built-in MATLAB function, so you probably shouldn't name your variable with the same name as it will interfere with that functionality for you and anyone that runs your program.

I forgot that ##\texttt{gamma}## is a built-in function!

kreil said:
The variable z only shows up in one place after it is defined: as an exponent to exp in the calculation of U1. I assume this U1 is the diffraction pattern that you want? (You might want to comment your code better so that it's clear what all the variables are.)

Yes, it is the diffraction pattern that I want. Sorry for not putting comments on my code.

kreil said:
The only other quantity of exp that could change is gamma, which depends on alpha and beta. Upon looking at gamma, I realized that your section of code that sets elements of gamma to zero is "dead" (seemingly not used):

Code:
if (alpha(i)^2 + beta(j)^2) > 1;
gamma(j, i) = 0;

When I ran your program, gamma has all nonzero values:

Code:
nnz(gamma)/numel(gamma)

ans =

1

So there might be something up there, if you were expecting some elements of gamma to be zero. But ultimately, if gamma is unchanged and you vary z, then the exp term gets bigger and bigger. I would expect the exponent to be negative here if you wanted the exp term to disappear as z increases.

The condition is to disregard evanescent waves as they do not contribute to the diffraction pattern. I guess these show-up at certain values of ##z##.
Anyway, I just resolved the issue by not using spatial frequency coordinates, thank you!

Ok, glad you were able to resolve the problem! Consider posting the solution here for the benefit of others.

kreil said:
Ok, glad you were able to resolve the problem! Consider posting the solution here for the benefit of others.

Instead of guessing the spatial frequency values of the Fourier transformed electric field, which is done after ##\texttt{U0}## and before the definition of ##\texttt{gamma}##, I consider this another form of the angular spectrum method:

##\mathbf{\hat{E}} \left(k_x, k_y; z\right) = \mathbf{\hat{E}} \left(k_x, k_y; 0\right) e^{\pm i k_z z}##

where ##k_z = \frac{k z}{\sqrt{x^2 + y^2 + z^2}}##, and ##k = \frac{2\pi}{\lambda}##. I then used the inverse Fourier transform to plot the electric field at ##x, y, z##.

Hi ecastro,

I modified your code, delete some unnecessary parts and add a figure show to visualize the results. And note that instead of
Code:
U1 = ifft2(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z));

you should use
Code:
U1 = ifft2(ifftshift(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z)));
instead, since you fftshift the spectrem first, you need to first iffshift it back then can you impose the fft2 on the scalar field. Although the results seem to be the same because the captured intensity is equals abs(U1.*conj(U1))

The modified code is following:
Code:
clc; clear; close all;

% parameters
layer = zeros(499, 499);
diameter = 100e-6;
lambda = 500e-9;            % wavelength (meters)
k = 2*pi/lambda;            % wavenumber
z = 2e-4;                      % altitude (meters)
phy_x = 1e-3;               % physical width (meters)
phy_y = 1e-3;               % physical length (meters)
obj_size = size(layer);

% generate meshgrid
dx = linspace(-phy_x/2, phy_x/2, obj_size(2));
dx = dx(1:length(dx));
dy = linspace(-phy_y/2, phy_y/2, obj_size(1));
dy = dy(1:length(dy));

% generate a circle aperture
for i = 1:obj_size(2)
for j = 1:obj_size(1)
if sqrt(dx(i)^2 + dy(j)^2) <= radius;
layer(j, i) = 1;
end
end
end

Fs_x = obj_size(2)/phy_x;
Fs_y = obj_size(1)/phy_y;
dx2 = Fs_x^(-1);
dy2 = Fs_y^(-1);
dFx = Fs_x/obj_size(2);
dFy = Fs_y/obj_size(1);
Fx = (-Fs_x/2:dFx:(Fs_x/2 - dFx));
Fy = (-Fs_y/2:dFy:(Fs_y/2 - dFy));

% alpha and beta (wavenumber components)
alpha = lambda.*Fx;
beta = lambda.*Fy;

% gamma_cust
gamma_cust = zeros(length(beta), length(alpha));
for j = 1:length(beta)
for i = 1:length(alpha)
if (alpha(i)^2 + beta(j)^2) > 1
gamma_cust(j, i) = 0;
else
gamma_cust(j, i) = sqrt(1 - alpha(i)^2 - beta(j)^2);
end
end
end

% angular spectrm based formula
U1 = ifft2(ifftshift(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z)));
I1 = (1/(16*pi)).*(U1.*conj(U1));

% show results
figure;     imshow(layer, []);      title('circular aperture');
figure;     imshow(I1, []);         title('angular spectrum based diffraction');

## What is the Code of the Angular Spectrum Method?

The Code of the Angular Spectrum Method is a computational algorithm used in signal processing to reconstruct a wavefield from its measured values. It is based on the principle of the angular spectrum representation of the wavefield, where the wavefield is represented as a superposition of plane waves with different propagation angles and wavenumbers.

## How does the Angular Spectrum Method work?

The Angular Spectrum Method works by first taking the measured values of a wavefield and applying a Fourier transform to obtain its spatial frequency components. These components are then multiplied by a phase term to take into account the propagation angle of each plane wave. The result is then inverse Fourier transformed to obtain the reconstructed wavefield.

## What are the advantages of using the Angular Spectrum Method?

One of the main advantages of the Angular Spectrum Method is that it can handle arbitrary boundary conditions, making it suitable for a wide range of applications. It is also computationally efficient, making it a popular choice for reconstructing wavefields in real-time applications.

## What are the limitations of the Angular Spectrum Method?

While the Angular Spectrum Method has many advantages, it also has some limitations. One limitation is that it assumes the wavefield is bandlimited, which may not always be the case in practical applications. It also requires a large amount of computational resources, making it less suitable for low-power or embedded systems.

## What are some applications of the Angular Spectrum Method?

The Angular Spectrum Method has a wide range of applications in fields such as optics, acoustics, and medical imaging. It is commonly used for wavefront sensing, imaging through turbid media, and holographic reconstruction. It also has potential applications in wireless communication and radar systems.

Replies
1
Views
5K
Replies
5
Views
2K
Replies
1
Views
2K
Replies
7
Views
2K
Replies
1
Views
2K
Replies
1
Views
3K
Replies
1
Views
2K
Replies
1
Views
18K
Replies
2
Views
2K
Replies
5
Views
7K