필터 지우기
필터 지우기

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 CenterFile Exchange에서 Discrete 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