フィルターのクリア

FFT of gaussian if radial domain

4 ビュー (過去 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 performfft 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:
Hope this helps!
Thanks
Balaji

カテゴリ

Help Center および File ExchangeDiscrete Fourier and Cosine Transforms についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by