function fftillustrator
% MATLAB script for one and two-dimensional FFTs
% R.F.Murphy, April 25, 2002

% generate a sine wave covering 0 to 2pi using 64 points
x=[0:63];
xn=x*2*pi/64;
xs=sin(xn);
% plot the sampled sine wave versus the data point number
figure(1);
plot(xs);

% take the fft of the sampled sine wave
xf=fft(xs);
% plot the real part of the fft in blue and the imaginary part in red
figure(2);
plot(x,real(xf),'b:',x,imag(xf),'r-');

% take the inverse fft of the fft
ixf=ifft(xf);
% compare the results of taking forward and inverse transform to the
% original data (the maximum of the absolute value of the difference
% should be small)
max(abs(real(ixf)-xs))

% make a copy of the fft and "zero out" the specific frequencies 
% corresponding to the original sin wave
xf2=xf;
xf2(2)=0;
xf2(64)=0;
% take the inverse of the "filtered" fft
ixf=ifft(xf2);
% the result should essentially be zero for all pixels
max(abs(real(ixf)))

% now generate an image which is a sin wave in both x and y
im=xs'*xs;
% and plot it with black corresponding to -1 and white corresponding to 1
figure(3);
imshow(im,[-1 1]);
     
% take the two-dimensional fft of this image and plot it 
imf=fft2(im);
% show that the imaginary part (corresponding to cosines) is essentially
% all zeros
min(min(imag(imf)))
max(max(imag(imf)))

% show the min and max of the real part (should be -1024 and 1024)
minr=min(min(real(imf)))
maxr=max(max(real(imf)))

% plot the fft as an image with -1024 mapped to black and 1024 mapped to white
figure(4);
imshow(real(imf),[minr maxr]);

% show the indices of the image where the min and max values are found
[i,j,v]=find(abs(imf)>1000)

% set the first non-zero term to zero, inverse transform and plot
imf(2,2)=0.;
imif=ifft2(imf);
figure(5);
imshow(real(imif),[-1 1]);

% set the second non-zero term to zero, inverse transform and plot
% (this should look like the original sin image but at half intensity
% since two of the four non-zero terms have been nulled out)
imf(2,64)=0.;
imif=ifft2(imf);
figure(6);
imshow(real(imif),[-1 1]);

% set the third non-zero term to zero, inverse transform and plot
% (this should look like a sin going from one corner diagonally to the other)
imf(64,2)=0.;
imif=ifft2(imf);
figure(7);
imshow(real(imif),[-1 1]);