1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ dsp.stackexchange.com might be a more appropriate place to ask this question. I haven't read it, but can you try taking $N$ larger, say $N = 128$? $\endgroup$ – snar Dec 4 '19 at 15:24
  • $\begingroup$ seems to me you compute an inverse FFT ifft2 $\endgroup$ – Ahmad Bazzi Dec 4 '19 at 15:29
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.