% FFT demo for image compression
% 
% for Octave: you will need to modify this file 

% load and display cat
I=imread('cheezburger.jpg');
X=im2double(rgb2gray(I)); % make grayscale and floating point

figure(1);
imagesc(X); 
colorbar; colormap gray;

% analyze frequencies with FFT in 2D
F = fft2(X);

% display spectrum with 0 frequency at origin
figure(2);
imagesc(log10(abs(fftshift(F))));
colorbar; colormap gray;

% drop all coefficients below a threshold
threshold = 1e2;
ii=find(abs(F)<threshold);
Fbis = F;
Fbis(ii)=0;
figure(3);
imagesc(log10(abs(fftshift(Fbis))));
colorbar; colormap gray;

% resynthetize image and compare
Xbis = ifft2(Fbis);
figure(4);
imagesc(Xbis);
colorbar; colormap gray;

% percentage reduction in size if we only needed to store the Fourier
% coefficients above the threshold. 
%
% Note: this ratio assumes that the storage for a Fourier coefficient is the
% same as that of an image pixel. This is not the case in this code but there
% are ways of storing the coefficients efficiently (for example: we don't need
% double precision to encode an image)
fprintf('compression ratio = %4.1f%%\n',...
 100*length(find(abs(F)>=threshold))/length(F(:)))