Why does my FFT highly depends on the resolution of the data?

1 view (last 30 days)
I am trying to implement fourier beam propagation method (FFT=BPM) in matlab. This algorithm requires making a fourier transform of an EM wave which is just a 2D matrix with different values.
When the input beam is a Gaussian exp(-r^2/2sigma^2) the transform returns another Gaussian. But, when I change the number of points in my matrix the transform gives very different values which are incorrect.
Basically, I am looking on the fft of my psi_init in the function that is in the bottom. When the resolution (1st line) is high, i.e res = 100, the result is almost a delta function - 10^8 in the origin and 0 almost everywhere else. When the resolution it is similar to a Gaussian.
Why do I get such different results? The fft of a guassian should return a guassian as well.
% define variables
res = 100; %resolution
z_res = 100; %time resolution
cm = 1e-02;
mili_m = 1e-03;
micro_m = 1e-06;
nano_m = 1e-09;
Lx = 4 * mili_m;
Ly = 4 * mili_m;
Lz = 20 * cm;
dl = Lx / res;
dz = Lz / z_res;
lambda = 445 * nano_m;
k0 = 2 * pi / lambda;
n0 = 1.82533;
alpha = 1;
kappa = 0.7;
beta = 14e-06;
X = linspace(0, Lx, res) - Lx / 2;
Y = linspace(0, Ly, res) - Ly / 2;
k_x = X;
k_y = Y;
% create H matrix
temp = (dz / (2 * k0 * n0)); %beta = k0 * n0
H_transfer = (exp(1i*temp*k_x.^2))'*(exp(1i*temp*k_y.^2));
H = fftshift(H_transfer);
% create initial wave function
fwhm = 5000 * micro_m;
w0 = fwhm/sqrt(2*log(2));
psi_init =(exp(-(1/w0^2) * X.^2))'*exp(-(1/w0^2) * Y.^2); %Gaussian pulse
abs_psi = (abs(psi_init)).^2;
imagesc(X,Y,abs(psi_init.^2)) % plots pulse before propagation
hold on
% create T matrix
BC = zeros(res);
T = abs(abs_psi);
T = relaxation(BC, abs_psi, T, res, dl);
imagesc(X,Y,T) % plots pulse after propagation
hold on
for i = 1:100
psi = bpm(T, psi_init, dz, k0, n0, H, beta, X, Y);
psi_init = psi;
imagesc(X,Y,abs(psi_init.^2)) % plots pulse after propagation
hold on
function [psi] = bpm(T, psi_init, dz, k, n0, H, beta, X, Y)
N = - beta .* T;
delta_n = 0.4;
s = -1i*k*N*dz;
S = exp(s);
z = fft2(psi_init);
zz = z .* H;
psi_prime = ifft2(zz);
psi = S .* psi_prime;

Answers (1)

Matt J
Matt J on 2 Oct 2022
Edited: Matt J on 3 Oct 2022
The fft of a guassian should return a guassian as well.
Remember, the FFT and the continuous Fourier transform are not the same thing. The FFT will approximate the continuous transform when,
(a) The sampling is sufficiently fine,
(b) The spatial extenf ot he sampling is sufficiently large to cover the signal, and
(c) the FFT is multiplied by the sampling interval lengths


Find more on Fourier Analysis and Filtering in Help Center and File Exchange


Community Treasure Hunt

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

Start Hunting!

Translated by