% Read the image
g = imread("image.png");
% Convert the image to grayscale and normalize into [0, 1]
g = double(im2gray(g)) / 255;

% Since h == w, get image size N
[h, w] = size(g);
N = h;

% Compute matrix U
U = zeros(h, w);
for n = 0:w-1
    for l = 0:h-1
        U(l+1, n+1) = exp(-2*pi*1i*(l*n/N)) / N;
    end
end

g_freq = U * g * U;

% Here we randomly mask some frequency components. 
% Later in this course, we will apply some other filters.
g_freq(2:10, :) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Complete the following lines to reconstruct the grayscale image 
% using the above frequency g_freq and U.
% You are not allowed to use the built-in MATLAB function "ifft2"

g_reconstructed = (N*U') * g_freq * (N*U');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Display original and reconstructed image
figure;
subplot(1, 2, 1);
imshow(g);
title("Original Image");
subplot(1, 2, 2);
imshow(real(g_reconstructed));
title("Reconstructed Image");

