The Fast Fourier Transform of this Kernel Blows Up

Question

Warning: This type of math is new to me, and I am self learning for my thesis, so I am coming at this question with a shakey understanding to begin with.

I am trying to take the Fast Fourier Transform of a 2D Gaussian Dispersal Kernel that looks like this:

$$ G(x,y) = \frac{1}{4\pi \mu} e^{-\frac{x^2 + y^2}{4\mu}} $$

The spatial domain that I am using is a 2-dimensional 7x7 grid in which each cell has length of 500m. I chose N = 8 because to my understanding using N = $ 2^{m} $ makes the FFT algorithm efficient, and choosing N = 8 divides my grid in 7x7 cells quite nicely - although this could change if it is needed. Below is the code that I use in my problem:

% Setting Parameters

xl = 1750; yl = 1750;            % Length of x-axis & y-axis
N = 8          ;                 % Number of intervals
dx = 2*xl/N; dy = 2*yl/N ;       % Width of each cell
x = linspace(-xl,xl-dx,N);       % Define the x-axis
y = linspace(-yl,yl-dy,N);       % Define the y-axis
[X,Y] = meshgrid(x,y);           % Create a spatial grid. 

% Dispersal parameters & Kernel

Dt = 0.667;                      % Estimate pulled from previous study
mu = Dt;
K2D=1/(4*pi*mu)*exp(-(X.^2+Y.^2)./(4*mu));

% Initialize a population in the (-1<=y<=1)x(-1<=x<=1) square
p0 = (abs(X)<=500 & abs(Y)<=500); 

% Perform the FFT on p0 and K2D and get p1.
fp0 = fft2(p0);
fK2D = dx * dy * fft2(K2D);
fp1 = fp0.*fK2D;
p1 = real( fftshift( ifft2(fp1) ) );

What seems to be of issue for me is that when I go to calculate fK2D, all of the values get scaled by a factor of 4, while the original kernel has a value that does not have this scaling. The problem is that when I go to calculate p1, which conceptually you can think of as the population of some animal after p0, the values retain the scaling and populations seem to blow up. From what I have seen, I am not sure whether or not this means that I need to "normalize" the kernel, and if it is the case that I do need to normalize the kernel, it's not clear to me what it should be normalized by.

Answer 1

One thing for sure is that it is not due to $\mu$. See the insights below. Note that MATLABs FFT performs a sum with no normalization, while the ifft does. Typing help fft on MATLAB tells us

For length N input vector x, the DFT is a length N vector X,
with elements
                 N
   X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                n=1
The inverse DFT (computed by IFFT) is given by
                 N
   x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
                k=1

Some insights on the FFT

It is worth noting the following FFT property: $$\mathcal{F}(e^{-at^2}) = \sqrt{\frac{\pi}{a}} \exp(- \pi^2 \frac{f^2}{a} )$$ The 2D FFT in your case will give a normalization factor (for $a = \frac{1}{4\mu}$) that looks like this $\sqrt{4\pi\mu}$ for one dimension. Square it you get $4\pi\mu$. This normalization factor cancels the one in your expression.