% Create a white box on a black background image
M= 256; N = 256;
image = zeros(M,N);
box = ones(64,64);
%box at centre
image(97:160,97:160) = box;

% Show Image

figure(1)
imshow(image)
title('Original Image');

% compute fft and display its spectra

F=fft2(double(image));
figure(2)
imshow(abs(fftshift(F)));
title('FFT Image Specta');

%compute Butterworth Low Pass Filter

u0 = 20; % set cut off frequency


u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);

for i = 1: M
    for j = 1:N
      %Apply a 2nd order Butterworth  
      UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
      
      H(i,j) = 1/(1 + UVw*UVw);
    end
end    


% display
figure(3)
imshow(fftshift(H))
title('Butterworth Low Pass Filter in Fourier Space');

% Apply filter and do inverse FFT
G=H.*F;
g=real(ifft2(double(G)));

% Show Result

figure(4)
imshow(g)
title('Butterworth Low Pass Filtered Image');