% Testing the fast Fourier transform
% in dimension 2 (for images)
%
% Samuli Siltanen November 2014

% Load image 
im = imread('Tokyo2048.tif','tif');
im = double(im);

% Normalize image: pixel values between zero and one
im = im-min(im(:));
im = im/max(im(:));

% Compute the FFT
Fim = fft2(im);

% Reorganize the FFT elements
Fim = fftshift(Fim);

% Put some of the frequency components to zero
%nvec        = [[-M/2:-1],[0:(M/2-1)]];
Fim_filtered = Fim;
Fim_filtered(1024+[-5:5],1024+[-5:5]) = 0;
im_filtered  = real(ifft(fftshift(Fim_filtered)));
im_filtered = im_filtered-min(im_filtered(:));
im_filtered = im_filtered/max(im_filtered(:));

% Prepare the FFT side for plotting
Fim_filtered_show = log(abs(Fg_filtered));
Fim_filtered_show = Fim_filtered_show-min(Fim_filtered_show(:));
Fim_filtered_show = Fim_filtered_show/max(Fim_filtered_show(:));

figure(1)
clf

% Plot the original image and the FFT
%imagesc([im,Fim_filtered_show,im_filtered],[.4,.6])
imagesc([im,im_filtered])
colormap gray
axis equal
axis off


