4

simple problem:

I plot out a 2D Gaussian function with a certain resolution in Matlab. I test with variance or sigma = 1.0. I want to compare it to the result of FFT(Gaussian), which should result in another Gaussian with a variance of (1./sigma). Since I am testing with sigma = 1.0, I would think that I should get two equivalent, 2D kernels.

i.e.

g1FFT = buildKernel(rows, cols, mu, sigma)    % uses normpdf over arbitrary resolution (rows, cols, 3) with the peak in the center

buildKernel:

function result = buildKernel(rows, cols, mu, sigma)  

result = zeros(rows, cols, 3);

center_w = floor(cols / 2);
center_h = floor(rows / 2);

for i = 1:rows
    for j = 1:cols
        distance = sqrt((center_w - j).^2 + (center_h - i).^2);
        g_val = normpdf(distance, mu, sigma);
        result(i, j, :) = g_val;
    end
end

% normalize so that kernel sums to 1
sumKernel = sum(result(:));
result = result ./ sumKernel;    

end

I am testing with mu = 0.0 (always), and variance or sigma = 1.0. I want to compare it to the result of FFT(Gaussian), which should result in another Gaussian with a variance of (1./sigma).

i.e.

g1FFT = circshift(g1FFT, [rows/2, cols/2, 0]);       % fft2 expects center to be in corners 
freq_G1 = fft2(g1FFT);
freq_G1 = circshift(freq_G1, [-rows/2, -cols/2, 0]); % shift back to center, for comparison's sake

Since I am testing with sigma = 1.0, I would think that I should get two equivalent, 2D kernels, because if sigma = 1.0, then 1.0/sigma = 1.0. So, g1FFT would equal freq_G1.

However, I do not. They have different magnitudes, even after normalization. Is there something I am missing?

8
  • Please provide all of your code to reproduce your error. You are correct in that the FFT of a Gaussian is another Gaussian, but with no code to show us what you did, we can't tell whether you're doing it properly.
    – rayryeng
    Jun 19, 2015 at 20:57
  • What does your buildKernel function look like?
    – rayryeng
    Jun 20, 2015 at 6:05
  • Really straightforward, but code added.
    – kdottiemo
    Jun 20, 2015 at 17:01
  • Why is your kernel a 3D matrix? You said it was a 2D Gaussian function, but there are three channels.
    – rayryeng
    Jun 20, 2015 at 17:32
  • 1
    You've probably found this already but you can't straightforwardly apply a continuous-time Fourier transform pair to discrete time/DFT. users.ece.gatech.edu/mrichard/… shows the equivalence for the discrete case depends on standard deviation being in terms of samples. Jun 23, 2015 at 14:37

1 Answer 1

0

To keep things simple, I will first cover the case for one-dimensional signals. Similar observations can be made for multi-dimensional cases.

The Fourier Transform of a continuous time Gaussian signal is itself a Gaussian function as indicated in this table. One can note that the wider the Gaussian in the time domain, the narrower the transformed Gaussian in the frequency domain and that for mu=0 and sigma=1/sqrt(2π) (which corresponds to a=1/(2*sigma^2)=π in the above transform table), the Fourier Transform of the continuous time function

enter image description here

would be the similar function (where only a change of variables occurred):

enter image description here

That's all good, but this is for a continuous time signal and we are really interested in discreet time signals. Unfortunately, and as also indicated on wikipedia, the Discrete Fourier Transform of a kernel obtained by sampling the continuous time Gaussian function, is not itself a sampled Gaussian function. Fortunately, this relationship is still often approximately true (without going into too much details, it requires the time-domain kernel to be wide enough but not too wide such that the frequency-domain approximation is also wide enough for the relationship to also be approximately true for the inverse transform). In this case, the Discrete Fourier Transform of the periodic extension (with period N) of the discrete time signal

enter image description here

where mu=0 and sigma=sqrt(N/2π) could be approximated by the similar function (up to a scaling factor and a change of variables):

enter image description here

You could then modify buildKernel to support different standard deviations sqrt(rows/2π) and sqrt(cols/2π) along the rows and columns respectively:

function result = buildKernel(rows, cols, mu, sigma)  

  if (length(mu)>1)
    mu_h = mu(1);
    mu_w = mu(2);
  else
    mu_h = mu;
    mu_w = mu;
  endif
  if (length(sigma)>1)
    sigma_h = sigma(1);
    sigma_w = sigma(2);
  else
    sigma_h = sigma;
    sigma_w = sigma;
  endif

  center_w = mu_w + floor(cols / 2);
  center_h = mu_h + floor(rows / 2);

  r = transpose(normpdf([0:rows-1],center_h,sigma_h));
  c = normpdf([0:cols-1],center_w,sigma_w);
  result = repmat(r * c, [1 1 3]);

  % normalize so that kernel sums to 1
  sumKernel = sum(result(:));
  result = result ./ sumKernel;    

end

which you could use to get a kernel whose FFT is a scaled version of itself. In other words a kernel obtained using

g1FFTin = buildKernel(rows, cols, mu, [sqrt(rows/2/pi) sqrt(cols/2/pi)]);

would be such that freq_G1 (as computed in your posted code) is nearly equal to g1FFTin * sqrt(rows*cols).

Finally given that your intention is really only to test that the kernel's FFT is also (approximately) Gaussian, you may wish to compare the FFT of a more arbitrary kernel with standard deviation sigma against another appropriately scaled Gaussian kernel computed directly in the frequency domain. In other words, assuming a spatial domain kernel obtained with:

g1FFTin = buildKernel(rows, cols, mu, sigma);

with corresponding frequency-domain representation obtained with:

g1FFT   = circshift(g1FFTin, [rows/2, cols/2, 0]);
freq_G1 = fft2(g1FFT);
freq_G1 = circshift(freq_G1, [-rows/2, -cols/2, 0]);

Then freq_G1 can be compared against another appropriately scaled Gaussian kernel computed directly in the frequency domain:

freq_G1_approx = (rows*cols/(2*pi*sigma^2))*buildKernel(rows, cols, mu, [rows cols]/(2*pi*sigma));
6
  • Thanks for your detailed explanation @SleuthEye. I think I mostly understand. I tried your first example (with sigma being [sqrt(rows/2/pi) sqrt(cols/2/pi)]) and found that scaling this kernel by sqrt(rows*cols) was approximately equal to the result calculated with fft. However, the second example didn't work for me. Scaling this next kernel by (rows*cols/(2*pi*sigma.^2)) did not (by far) equal the result calculated with fft. Any suggestions? I would like this to work for an arbitrary sigma, as you said.
    – kdottiemo
    Jun 30, 2015 at 21:55
  • @kdottiemo Just to clarify that it should work for more arbitrary standard deviations greater than ~0.5, and less than ~min(rows,cols)/3. Also, in case it was not clear in my post, freq_G1_approx should be compared against the result of the fft on buildKernel(rows,cols,mu,sigma) (followed by the circshift). So the second example would use two different calls to buildKernel: one constructed in the space domain and another directly in the frequency domain.
    – SleuthEye
    Jul 1, 2015 at 1:14
  • Thanks again. But I think I'm doing what you're saying, please let me know if I am still misunderstanding... %First kernel sigma = 2.0; g1 = buildKernel(rows, cols, 0.0, [rows, cols]/(2*pi*sigma)); scale = (rows*cols/(2*pi*sigma.^2)); scale = ones(size(g1)).*scale; new_G1 = g1.*scale; %Next kernel g1FFT = circshift(g1, [rows/2, cols/2, 0]); freq_G1 = fft2(g1FFT); freq_G1 = circshift(freq_G1, [-rows/2, -cols/2, 0]); %So, in this case, freq_G1 != new_G1
    – kdottiemo
    Jul 1, 2015 at 2:10
  • see update. the kernel width in spatial & freq domain only match for the very specific case of the 1st example. For 2nd example, the width won't match (so kernel is not its own transform). Following your comment, keeping the first kernel, but changing %Next kernel g1prime = buildKernel(rows,cols,mu,sigma); g1primeFFT = circshift(g1prime,...; freq_G1 = fft2(g1primeFFT); freq_G1 = circshift.... Then freq_G1 == new_G1.
    – SleuthEye
    Jul 1, 2015 at 2:50
  • Thanks; that cleared it up. This is so counter-intuitive to me - so to achieve the shape that one would expect with an input sigma, sigma itself actually becomes [rows cols]/2*pi*(sigma). Am I over-simplifying things by saying that? And the relation to its frequency-space representation given that same sigma is a scaling of ((rows * cols)/2*pi*(sigma^2))?
    – kdottiemo
    Jul 2, 2015 at 1:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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