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
would similar function (where change of variables occurred):
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
where mu=0 , sigma=sqrt(n/2π) approximated similar function (up scaling factor , change of variables):
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
Post a Comment