Fourier Analysis using fourier function for bessel distribution.

Hello,

I am trying to write a code for fourier analysis from spatial domain to frequency domain.

I am using "fourier" and "ifourier" function but I am not getting the desired result.

I am sharing the code here. Please help me where I am getting wrong.

The desired result is;

Phase plot: is a spiral distribution

Amplitude Plot: An intensity NULL at the center with ring distribution around it for the amplitude distribution.

Thanks,

%% Initialize MATLAB
clc;
clear;
close all;
%% Define wave parameters(lambda is taken as per paper)
lambda = 632.8e-9;
k = 2*pi/lambda;
fx = 1.9048e-04;
fy = 1.9048e-04;
%% Define Symbolic variable
syms x y z

% Define r
r = sqrt(x^2 + y^2);


%% Bessel function of order 1, 1st kind.
% Define argument of the bessel function
arg = k*r;

% Define Bessel function
B = besselj(1,arg);

% Calculate Fourier transform of the above bessel function
FT_B = fourier(B);
%% Define spatial frequency transfer function
H = exp(1i*k*z)*exp(-1i*pi*lambda*z*(fx^2 + fy^2));

%% Compute Inverse Fourier transform
IFT_B = ifourier(FT_B*H)

IFT_B = 

%% Discretize domain
% Define Total no. of discrete points
Np = 1051;

% Define co-ordinate axis
x_axis = linspace(-100e-3,100e-3,Np);
y_axis = linspace(-100e-3,100e-3,Np);

% Define distance z at which the bessel function is to be calculated.
z_fixed = 5e3;

% Defining mesh grid
[X,Y] = meshgrid(x_axis,y_axis);      
Y = (-1)*Y;


% Calculate the expression obtained form the inverse fourier transform
% operation with symbolic expression
I = exp((z_fixed*2665344047195405i)/268435456)*besselj(1, (2665344047195405*(X.^2 + Y.^2)^(1/2))/268435456);
%% Plot Bessel function
% Plot amplitude distribution at z = 5 km
figure
surf(x_axis,y_axis,abs(I));
view(0,90)
shading interp;
colorbar;
xlabel('x-axis in m');
ylabel('y-axis in m');
title('Amplitude Distribution of the Bessel beam at z =  5 km');

% Plot phase distribution at z = 5 km
figure
surf(x_axis,y_axis,angle(I));
view(0,90)
shading interp;
colorbar;
xlabel('x-axis in m');
ylabel('y-axis in m');
title('Phase Distribution of the Bessel beam at z =  5 km');

Answer 1

[Rangesh]

Hello Biplob,

I see that you are seeking to understand why the expected results for the phase and amplitude plot do not align with the Bessel function. Here are a few potential issues I found in the code:

• The Fourier transform of variable ‘B’ involves two symbolic variables. By default, function "fourier" determines the independent variable as x and provides a transformation variable w.
• To address this, I recommend using r as a symbolic variable, where “r = (x^2 + y^2)”, as r is treated as a unified value throughout the code and specifying the independent and transformation variable.

    syms r w; 
    FT_B=fourier(B,r,w) 

• Similarly, when using the function “ifourier”, I suggest specifying the independent and transformation variable.

     IFT_B = ifourier(FT_B * H, w, r); 

• Given that z is set to 5000, the transfer function H is a constant. As a result, the inverse Fourier transform would produce another Bessel function, with some scaling factor H.
• This suggests that the amplitude plot would exhibit a null value at the center, displaying a ring distribution, and the phase plot would show all values as null. Therefore, I suspect there might be an issue with the specified transfer function.
• To evaluate the symbolic expression, you can use the "subs" function:

     radius=X.^2+Y.^2  
     subs(IFT_B, r, radius); 

For better understanding, kindly refer to the following link:

• fourier: https://www.mathworks.com/help/symbolic/sym.fourier.html
• ifourier: https://www.mathworks.com/help/symbolic/sym.ifourier.html
• subs: https://www.mathworks.com/help/symbolic/subs.html

By following the suggestions, I believe that you can obtain the expected results for the plots. I hope this resolves your query.

Thanks,

Rangesh.