Fourier Transform of Gaussian Kernel in Matlab

Hi everyone,

I need your help!

I was reading a document on Discrete Fourier Transform and found the following example:

Here's the mentioned gaussian kernel:

g_k = (1/256)*[1 4 6 4 1;...
               4 16 24 16 4;...
               6 24 36 24 6;...
               4 16 24 16 4;...
               1 4 6 4 1];

As can be seen, the size of g_k or g(x,y) is 5 x 5 - while the size of G(u,v) is around 380 x 450

Can you please show me the Matlab code that can generate the above G(u,v) image or result?

Answer 1

[Matt J]

If you download gaussfitn from

https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn

then you can do,

g0=[1 4 6 4 1]';
params=gaussfitn((-2:2)',g0,[],{0,[],0},{0,[],0});

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

[A,mu,sig2]=deal(params{2:4});
sig=1/2/pi/sqrt(sig2);

fun=@(x) exp(-(x/sig).^2/2);

Gu=fun(linspace(-6*sig,+6*sig,450));


Guv=imresize(Gu'*Gu,[380,450]);

imshow(Guv);

Answer 2

[Matt J]

One could also do as below. This gives only an approximately Gaussian spectrum, however,

g_k = (1/256)*[1 4 6 4 1;...
               4 16 24 16 4;...
               6 24 36 24 6;...
               4 16 24 16 4;...
               1 4 6 4 1];

Guv=fftshift( abs(fft2(g_k,380,450)) );

imshow(Guv)

[Gobert]

Thank you! This is easier to understand.

[Gobert]

@Matt J - I have a question. When I follow instructions given in books, I always end up with an image that does not look like the expected better quality or improved version of the input image. See below. Can you help or tell me where I made a mistake?

img_in = rgb2gray(imread("peppers.png"));
figure, imshow(img_in,[]), title('input')

g_k = (1/256)*[1 4 6 4 1;...
    4 16 24 16 4;...
    6 24 36 24 6;...
    4 16 24 16 4;...
    1 4 6 4 1];

% Compute the FT of the input image
ft_in = fft2(img_in);

% Compute the FT of the gaussian kernel 
[n,m] = size(img_in);
ft_g_k = fftshift(abs(fft2(g_k,n,m)));

% Multiplication 
ft_mult = ft_g_k.*ft_in ;

% Inverse FT
inv_out = ifft2(ft_mult);
figure, imshow(inv_out,[]), title('output')

[Matt J]

Watch your zero-padding.

img_in = double(rgb2gray(imread("peppers.png")));
[m,n]=size(img_in);

g_k = (1/256)*[1 4 6 4 1;...
    4 16 24 16 4;...
    6 24 36 24 6;...
    4 16 24 16 4;...
    1 4 6 4 1];

% Compute the FT of the input image
ft_in = fft2(img_in,2*m,2*n);

% Compute the FT of the gaussian kernel
g_k(2*m,2*n)=0;%zero-pad
ft_g_k = fft2(circshift(g_k,[-2,-2]));

% Multiplication 
ft_mult = ft_g_k.*ft_in ;

% Inverse FT
inv_out = ifft2(ft_mult,'symmetric');
inv_out=inv_out(1:m,1:n); %unpad

figure, imshow(inv_out,[]), title('output')

[Gobert]

Thank you so much @Matt J

[Paul]

@Matt J

Why does the the zero-padded g_k need to be circshifted after zero-padding? Is it because the spatial indices of the upper left corner of g_k are (-2,-2) and the spatial indices of the upper left corner of img_in are (0,0)?

[Matt J]

@Paul Yes, that is why. If you don't do it, you will see a 2-pixel shift in the output.

[Gobert]

@Paul - These simple examples can give an idea of what "circshift" does to that "zero-padded matrix": See "ans"

% When m~=n
A = [16 2 3 13; 5 11 10 8; 9 7 6 12]
[n,m]=size(A)
A(2*n,2*m)=0
circshift(A,[-2,-2])

A =

16    2     3    13
5    11    10     8
9     7     6    12


n =

3


m =

4


A =

16    2     3    13     0     0     0     0
5    11    10     8     0     0     0     0
9     7     6    12     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0


ans =

6    12     0     0     0     0     9     7
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
3    13     0     0     0     0    16     2
10    8     0     0     0     0     5    11


% When m=n
B = [16 2 3 13; 5 11 10 8; 9 7 6 12;4 14 15 1]
[n,m]=size(B)
B(2*n,2*m)=0
circshift(B,[-2,-2])

B =

16    2     3    13
5    11    10     8
9     7     6    12
4    14    15     1


n =

4


m =

4


B =

16    2     3    13     0     0     0     0
5    11    10     8     0     0     0     0
9     7     6    12     0     0     0     0
4    14    15     1     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0


ans =

6    12     0     0     0     0     9     7
15    1     0     0     0     0     4    14
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
3    13     0     0     0     0    16     2
10    8     0     0     0     0     5    11