Matlab - FFT of Gaussian - Equivalency -


simple problem:

i plot out 2d gaussian function resolution in matlab. test variance or sigma = 1.0. want compare result of fft(gaussian), should result in gaussian variance of (1./sigma). since testing sigma = 1.0, think should 2 equivalent, 2d kernels.

i.e.

g1fft = buildkernel(rows, cols, mu, sigma)    % uses normpdf on arbitrary resolution (rows, cols, 3) peak in center 

buildkernel:

function result = buildkernel(rows, cols, mu, sigma)    result = zeros(rows, cols, 3);  center_w = floor(cols / 2); center_h = floor(rows / 2);  = 1:rows     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 kernel sums 1 sumkernel = sum(result(:)); result = result ./ sumkernel;      end 

i testing mu = 0.0 (always), , variance or sigma = 1.0. want compare result of fft(gaussian), should result in gaussian variance of (1./sigma).

i.e.

g1fft = circshift(g1fft, [rows/2, cols/2, 0]);       % fft2 expects center in corners  freq_g1 = fft2(g1fft); freq_g1 = circshift(freq_g1, [-rows/2, -cols/2, 0]); % shift center, comparison's sake 

since testing sigma = 1.0, think should 2 equivalent, 2d kernels, because if sigma = 1.0, 1.0/sigma = 1.0. so, g1fft equal freq_g1.

however, not. have different magnitudes, after normalization. there missing?

to keep things simple, first cover case one-dimensional signals. similar observations can made multi-dimensional cases.

the fourier transform of continuous time gaussian signal gaussian function indicated in this table. 1 can note wider gaussian in time domain, narrower transformed gaussian in frequency domain , mu=0 , sigma=1/sqrt(2π) (which corresponds a=1/(2*sigma^2)=π in above transform table), fourier transform of continuous time function

enter image description here

would similar function (where change of variables occurred):

enter image description here

that's good, continuous time signal , interested in discreet time signals. unfortunately, , indicated on wikipedia, discrete fourier transform of kernel obtained sampling continuous time gaussian function, not sampled gaussian function. fortunately, relationship still approximately true (without going details, requires time-domain kernel wide enough not wide such frequency-domain approximation wide enough relationship approximately true inverse transform). in case, discrete fourier transform of periodic extension (with period n) of discrete time signal

enter image description here

where mu=0 , sigma=sqrt(n/2π) approximated similar function (up scaling factor , change of variables):

enter image description here

you modify buildkernel support different standard deviations sqrt(rows/2π) , sqrt(cols/2π) along rows , 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 kernel sums 1   sumkernel = sum(result(:));   result = result ./ sumkernel;      end 

which use kernel fft scaled version of itself. in other words kernel obtained using

g1fftin = buildkernel(rows, cols, mu, [sqrt(rows/2/pi) sqrt(cols/2/pi)]); 

would such freq_g1 (as computed in posted code) equal g1fftin * sqrt(rows*cols).

finally given intention test kernel's fft (approximately) gaussian, may wish compare fft of more arbitrary kernel standard deviation sigma against another appropriately scaled gaussian kernel computed directly in frequency domain. in other words, assuming 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 compared against another appropriately scaled gaussian kernel computed directly in frequency domain:

freq_g1_approx = (rows*cols/(2*pi*sigma^2))*buildkernel(rows, cols, mu, [rows cols]/(2*pi*sigma)); 

Comments

Popular posts from this blog

powershell Start-Process exit code -1073741502 when used with Credential from a windows service environment -

twig - Using Twigbridge in a Laravel 5.1 Package -

c# - LINQ join Entities from HashSet's, Join vs Dictionary vs HashSet performance -