discrete Fourier transform of Gaussian

I'm trying to understand the following code, where we calculate the discrete Fourier transform of the Gaussian kernel:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n1=64;
n2=64;
n3=64;
i = sqrt(-1); 
wx=exp(i*2*pi/n1); % nth root of unity.
wy=exp(i*2*pi/n2); % nth root of unity.
wz=exp(i*2*pi/n3); % nth root of unity.
[x,y,z]=meshgrid([1:n1],[1:n2],[1:n3]); 
x=permute(x,[2 1 3]); 
y=permute(y,[2 1 3]); 
z=permute(z,[2 1 3]); 
KERNEL = n1*n1*(2-wx.^(x-1)-wx.^(1-x)) + n2*n2*(2-wy.^(y-1)-wy.^(1-y)) + ...
    n3*n3*(2-wz.^(z-1)-wz.^(1-z));
KERNEL = exp( -dt*KERNEL );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The Gaussian kernel is defined as follows: .

Let and and grid points .

The discrete Fourier transform (1D) of a grid function is the coefficient vector with .

But here in the code we compute the kernel in a different way. Does somebody know why the Fourier transform of the kernel can be computed this way?

What I mean is, why do we get constants (I write n for , ) because I thought it should be since we are in the 3D case. And I don't understand why we compute , because we then compute the kernel in a range since x takes the values instead of and this isn't the range from the discrete Fourier transform definition. I'd really appreciate if someone could give me a hint to understand this computation.