2D FFT Gaussian filtering exercise doesn't return what is expected

Hi, for the purpouse of exercise I have been trying to program my own Gaussian filter in fourier domain. So far I have googled and talked to chatGPT with no luck. Lots of codes online I found just use FOR LOOP for simple frequency cut-off filters, which is insufficient in this case. The result is an image which seems to be overlayed with the original one and flipped 90 degrees each time. It does that even without the filter.

One of the Fourier characteristics is that Convolution in real space is equal to Multiplication in fourier space:

So I wanted to test this and created a filter same size as the original image (see the function "gauss_distribution_2" below), and then just multiplied it element by element with the absolute value of fourier transformed image. Could you please tell me what am I doing wrong?

%% -- function I programmed --
function gauss_distribution_2 = gauss_distribution_2(m,n,sigma)
% CREATES 2D MATRIX WITH GAUSS DISTRIBUTED DATA IN FOURIER DOMAIN NORMED FOR 1
[x, y] = meshgrid(-floor(n/2):floor((n-1)/2), -floor(m/2):floor((m-1)/2));
r = x.^2 + y.^2;
gauss_distribution_2 = sqrt(pi/sigma)*exp(-r.^2 /(4*sigma));
gauss_distribution_2 = gauss_distribution_2./max(gauss_distribution_2(:));
end

%% -- code --
im = imread("image.png");
im = rgb2gray(im);
m = size(im, 1); n = size(im, 2); sigma = 1e8; % filter values
% FT
G = gauss_distribution_2(m,n,sigma); 
im_freq = fft2(im);
im_freq = fftshift(im_freq);
im_freq = abs(im_freq);
filtered_spectr = G .* im_freq; % filtration
% spectrum seems to be filtered alright
figure; imagesc(double(log(abs(filtered_spectr)+1)));colormap gray; axis equal; grid off; xticks([]); yticks([]) 
filtered_spectr(filtered_spectr<1e-10) = 0; % perhaps numbers were too small?
% Inverse FT
filtered_spectr = ifftshift(filtered_spectr);
im_filtered = ifft2(filtered_spectr);
im_filtered = abs(im_filtered);
% image returned is nonsense
figure, imagesc(uint8(im_filtered)), colormap gray, title('filtered image'), axis equal;

Answer 1

[Bruno Luong]

You made many short cuts and reasoning that are not right.

The convolution formula you post is for continuous signal/surface, define on the entire R or R^2.

The convolution on image processing is discrete convolution, with truncated domain. There is an equivalent theorem, but it applies with periodic discrete signal. So the product of discrete DFT is actually the DFT of the wrap around signal. See Circular convolution theorem and cross-correlation theorem in this page.

Furthermore the DFT of the truncated discrete Gaussian signal is NOT a Gaussian in Fourier domain.

All that make your calculation returns an unexpected filtered image.

Take a look at this FEX to see how to perform MATLAB discrete convolution using DFT

[Josef]

I'm not orienting myself in the topic into such depth. Although thank you for the hints. Would you have more specific advice for a beginner? Is it just better to FT the Gaussian and then multiply or that won't work either?

[Bruno Luong]

You have first and foremost to pad the image and kernel with 0s to avoid wrapping around artefact. This can be done with the second argument of FFT without explicitly padding. After ifft you have to trim back to get the original image size.

Yes working on FFT of the kernel do not assume that the kernel of a Gaussian is the Gaussian due to border effect.

Look at the file exchange I did all that with code.

[Josef]

So I foundmy mistake. The abs value can not get filtered.

im_freq = abs(im_freq); 
filtered_spectr = G .* im_freq; % filtration

I rewrote the code into a new script and it worked just fine. So I guess my assumption was correct after all. I have compared the result to imgaussfilt(im,5). The results seem to be quite similar - my function yields a bit more contrast. The spectra is however different.

im = imread("image.jpg");
% Parameters
m = size(im, 1);   
n = size(im, 2);    
sigma = 1e6;   
% FFT
im_freq = fft2(im); 
im_freq = fftshift(im_freq); 
% Filtration
G2 = gauss_distribution_2(m,n,sigma);
imagesc(G2); colorbar;
filtered_spectr = im_freq.*G2;
% inverse FFT
im_filtered = ifftshift(filtered_spectr); % fftshift zpt
im_filtered = ifft2(im_filtered); % inverse fft2
im_filtered = abs(im_filtered);
figure, imshow(uint8(im_filtered)), colormap gray, title('filtered image my function'), axis equal;
B = imgaussfilt(im, 5);
figure, imshow(uint8(B)), colormap gray, title('filtered image imgaussfilt'), axis equal;