FFT of gaussian if radial domain

11 次查看（过去 30 天）
elis02
elis02 2023-5-20
回答： Balaji 2023-9-8
Hi
I defind a gaussian with cartesian dimain as written below.
How do I make the same but with radial? without the angle (theta).
In other words I want to transfer instead of x,y to 'r' :
exp(-(r(i_index)^2 )/w0^2)
and keep everything.
% Constants and Parameters
c = 299792458; % Speed of light [m/s]
twidth = 500e-15; % Time window [s]
%twidth = 5e-12;
n = 2^10; % Number of points
lambda0 = 795e-9; % Center wavelength [m]
W0 = 2*pi*c/lambda0; % Center frequency [rad/s]
lambda0_um = 1e6*lambda0;
FWHM = 34e-15; % Initial pulse duration [s]
Epsilon_0 = 8.8541878128*1e-12;
% Time domain
T = linspace(-twidth/2, twidth/2, n); % Time grid
dt = T(2) - T(1);
Dnu = 1/dt;
delta_T = twidth/n; % Time interval [s]
% Gaussian pulse in time domain
y = exp(1i*W0*T); % Carrier
sigma = FWHM/((2*log(2))^0.5);
gaussian = exp(-T.^2/sigma^2);
N = sqrt(sum(gaussian.^2)*dt);
gaussian = gaussian/N;
fin_t = y.*gaussian; % Electric field in time domain
% Spatial domain
x = linspace(-10e-3, 10e-3, 256); % at the focus, beam will be ~300um, so step size should be there approx 30um, not including self focusing
dx = x(2) - x(1);
w0 = 13e-3/2;
gaussian_space = zeros(length(x),length(x));
for i_index = 1:length(x)
for j =1:length(x)
gaussian_space(i_index,j) = exp(-(x(i_index).^2 + x(j).^2)/w0^2);
end
end
N = sqrt(sum(sum(gaussian_space.^2))*(dx)^2);
gaussian_space = gaussian_space / N;
Average_power = 14; % watt
Laser_operation_freq = 5e3; % Laser operation freq
Pulse_energy = Average_power / Laser_operation_freq;
Total_initial_intensity = 0.5*c*Epsilon_0*(sum(sum(gaussian_space.^2))*(dx)^2) * (sum(gaussian.^2)*dt);
Normalization = sqrt(Total_initial_intensity/Pulse_energy);
% Electric field in spatial and time domain
E_x_y_t = (1/Normalization)*gaussian_space.*reshape(fin_t,1,1,[]);
% Create frequency grids
frequencies = linspace(-Dnu/2, Dnu/2, n); % Frequency grid (time dimension)
n_x = length(x);
dfx = 1/(x(end)-x(1));
dfy = 1/(x(end)-x(1));
fx = linspace(-dfx*(n_x/2), dfx*(n_x/2), n_x); % Frequency grid (x dimension)
fy = linspace(-dfy*(n_x/2), dfy*(n_x/2), n_x); % Frequency grid (y dimension)
wavelength_freq_um = 1e6*c./frequencies; %wavelength in um
%% First Lens
first_focal_length = 4; %focal length in [m]
[i_index, j, k] = ndgrid(1:length(x), 1:length(x), 1:length(frequencies));
phase_first_focal_length = (2.*pi*frequencies(k)./c).*(x(i_index).^2+x(j).^2)/(2*first_focal_length);
%% now let's add the lens phase and return to the time domain
E_f = fftshift(fft(E_x_y_t, [], 3),3); % Frequency domain signal
E_f = E_f.*exp(-1i*phase_first_focal_length);
E_x_y_t = ifft(ifftshift(E_f,3), [], 3);
later i also add a different phase which correspond to the fourier of x,y:
phase_shift = sqrt(k_w_air(k).^2 - (2*pi*fx(i_index)).^2 - (2*pi*fy(j)).^2) * L_first_air_path;
delta_f = frequencies(k) - c/lambda0;
E_fx_fy_f_shifted = E_fx_fy_f_shifted .* exp(1i * phase_shift) .* exp(-1i*2*pi*delta_f *L_first_air_path/Vg_air);
Where fx , fy are the corresponding fourier of x, y so idealy should change to fr

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

回答（1 个）

Balaji
Balaji 2023-9-8
Hi Elis,
I understand that you want to perform ‘fft’ in radial domain.
Firstly you can change the Electric Field to radial domain as follows:
r = 0:dr:r_max;
n_r = length(r);
frequencies = linspace(-Dnu/2, Dnu/2, n);
% Electric field in radial and time domain
E_r_t = zeros(n_r, n);
for i = 1:n_r
radial_profile = exp(-r(i)^2/w0^2);
E_r_t(i, :) = radial_profile * fin_t;
end
Secondly for implementing fft in radial domain please refer to the following MATLAB answer:
How do you do a radial Fourier transform in MATLAB? - MATLAB Answers - MATLAB Central (mathworks.com)
Hope this helps!
Thanks
Balaji

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by