Matlab - FFT of Gaussian - Equivalency

1k views Asked by At

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?

1

There are 1 answers

6
SleuthEye On BEST ANSWER

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