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

Hi,

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;
figure
imagesc(X,Y,abs(psi_init.^2)) % plots pulse before propagation
hold on
colorbar
title('before')
%
% create T matrix
BC = zeros(res);
T = abs(abs_psi);
T = relaxation(BC, abs_psi, T, res, dl);
figure
imagesc(X,Y,T) % plots pulse after propagation
hold on
colorbar
title('T')

%
for i = 1:100
    psi = bpm(T, psi_init, dz, k0, n0, H, beta, X, Y);
    psi_init = psi;
    disp(i)
end
figure
imagesc(X,Y,abs(psi_init.^2)) % plots pulse after propagation
hold on
colorbar
title('after')

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;
end

Answer 1

[Matt J]

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