% This tutorial assumes basic comfort level with complex
% numbers (Z = A + Bi, where A and B are real numbers, and i is 
% the square root of -1), spatial frequency (measured in cycles 
% per length, e.g. cm-1).  
% It is aimed at explaining how Matlab deals with FFT's.

%% A test image to play with, saved as zebras.mat
load zebras

%% Starting with the one-dimensional case
line = zebras(128,:);  % just one line through the middle

%% The 1D FFT algorithm returns a vector with as many complex-
%% valued Fourier coefficients as there are points in the input.
f = fft(line);
whos line f   % just get Matlab to show how big and what
              % type of variable

%% Each Fourier coefficient tells the amplitude and the
%% phase of one spatial frequency component.  The FFT algorithm
%% returns real and imaginary components; amplitude is the
%% square root of the sum of squares of the real and imaginary
%% component; and phase is the angle between the two

%% Picking just one spatial frequency component (the one that
%% represents 2 cycles per window):  
f(3)
real(f(3))%  the real part of the complex number
imag(f(3))%  the imaginary part of the complex number
%% (The first component is the DC value, or mean intensity.)
%% We rarely care about scaling everything, and every algorith
%% does it differently, but Matlab scales the FT by N, the number
%% of points in the original function
mean(line(:))  % the average of the orignal function
abs(f(1))/length(line)  % the magnitude of the 0th FT component,
                        % divided by N, the number of points
                        
%% Show that the amplitude (abs) is sqrt(Real + Imag)
abs(f(3))  % I use amplitude and magnitude interchangeably to 
            % to describe the length of a complex number when
            % it is plotted as a vector in the complex plane 
            % (which is a coordinate system where the x-axis 
            % is the real part of the number and the y-axis
            % is the imaginary part of the number
sqrt(real(f(3)).^2 + imag(f(3)).^2)

%% And show that the phase is the angle between the two
phase(f(3))
atan(imag(f(3))/real(f(3)))

%% Or the imaginary part of the complex logarithm
imag(log(f(3)))  % just for fun ...

%% To look at the entire spectrum
subplot(4,1,1); plot(line)  % just to see what it looks like
xlabel('pixel number'); ylabel('intensity');
title('line through zebra image')

%% Plotting just the amplitude spectrum
subplot(4,2,3)
plot(abs(f)); xlabel('cycles per line'); ylabel('amplitude')

%% And the phase spectrum
subplot(4,2,4)
plot(angle(f)); xlabel('cycles per line'); ylabel('phase')
%% Except in rare cases, the "phase spectrum" looks like
%% noise, uniformly distributed between -pi and pi.  But the
%% actual values are important and are what generate edges
%% and features in the image.

%% Usually, the DC component has a lot of power, relative to
%% the rest of the spectrum, so it's easier to visualize if
%% we subtract the mean out of the line before doing this.
subplot(4,1,3)
f_dc0 = fft(line - mean(line(:)));
plot(abs(f_dc0)); xlabel('Fourier component');
ylabel('amplitude')
axis tight; title('Fourier transform of original function')

%% The second half of the Fourier vector doesn't represent
%% high spatial frequencies.  For Fourier transforms calculated
%% from real-valued images or vectors, there's a symmetry in 
%% the spectrum - the second half is equal (but opposite in phase)
%% to the first half. One of the weird things about the FFT 
%% algorithm is the order in which Fourier components are 
%% written into the resulting matrix.  The fftshift function is 
%% used to center the FFT:
f = fftshift(fft(fftshift(line-mean(line(:)))));
%% Now the symmetry is apparent, and we can plot and label 
%% our axes accordingly. 
subplot(4,1,4)
plot(-128:127,abs(f)); xlabel('cycles per image'); ylabel('amplitude')
axis tight;

%% A negative spatial frequency component can be interpreted as a statement
%% about relative phase.  A sine wave has a certain phase; a cosine
%% function has a 90 degree phase shift.  In the Fourier transform, these
%% are represented as the sum of positive and negative frequency components
%% with different phases (different real and imaginary components).
clf
x = -64:63;
sinExample = sin(2*pi*x/12.8); % this function will oscillate every 12.8 points
subplot(3,3,1); plot(x,sinExample); axis tight; title('odd symmetry')
cosExample = cos(2*pi*x/12.8);
subplot(3,3,2); plot(x,cosExample); axis tight; title('even symmetry')
mixExample = cos(2*pi*x/12.8 - 2);  % add a random phase shift so it's neither
subplot(3,3,3); plot(x,mixExample); axis tight; title('shifted cos')
  
%% The highest spatial frequency that can be represented in a
%% digitized image or vector is one cycle every two points.  This
%% is why the highest Fourier component of a 128 element line
%% is 64 cycles per line (window). The spacing of spatial frequency 
%% values (spatial frequency resolution) is determined
%% by the size of the image or vector: 1/FOV.  That's why the 
%% finest increment is one cycle per image; larger spatial frequencies
%% can't be represented accurately in an image.
freq = -64:63;
subplot(3,3,4); plot(freq,real(fftshift(fft(fftshift(sinExample)))),'r'); hold on
 plot(freq,imag(fftshift(fft(fftshift(sinExample)))),'g'); hold off
  xlabel('cycles per window'); title('FT of odd-symmetric function')
  axis tight; legend('real','imaginary')
subplot(3,3,5); plot(freq,real(fftshift(fft(fftshift(cosExample)))),'r'); hold on
 plot(freq,imag(fftshift(fft(fftshift(cosExample)))),'g'); hold off; axis tight
  xlabel('cycles per window'); title('FT of even-symmetric function')
subplot(3,3,6); plot(freq,real(fftshift(fft(fftshift(mixExample)))),'r'); hold on
 plot(freq,imag(fftshift(fft(fftshift(mixExample)))),'g'); hold off; axis tight
  xlabel('cycles per window'); title('FT of shifted function')  
%% The ringing artifact is due to the truncation and undersampling - 
%% artifacts of discrete Fourier transforms.  But the key point is that the
%% even-symmetric (cosine) function has no imaginary component; the
%% odd-symmetric function has no real component.  Something in between has
%% real and imaginary parts to it's FT.  All of the functions have the same
%% amplitude spectrum, but different phase spectra.  

subplot(3,3,7); plot(freq,10*abs(fftshift(fft(fftshift(sinExample))))/length(x),'k'); hold on
  plot(freq,angle(fftshift(fft(fftshift(sinExample)))),'b'); hold off
  xlabel('cycles per window'); title('Amp. and phase of sin function')  
  axis tight; legend('magnitude','phase')
subplot(3,3,8); plot(freq,10*abs(fftshift(fft(fftshift(cosExample))))/length(x),'k'); hold on
  plot(freq,angle(fftshift(fft(fftshift(cosExample)))),'b'); hold off;
  xlabel('cycles per window'); title('Amp. and phase of cos function'); axis tight
subplot(3,3,9); plot(freq,10*abs(fftshift(fft(fftshift(mixExample))))/length(x),'k'); hold on
  plot(freq,angle(fftshift(fft(fftshift(mixExample)))),'b'); hold off; axis tight
  xlabel('cycles per window'); title('Amp. and phase of shifted function')  

%% High spatial frequencies oscillate more rapidly; low spatial frequencies
%% oscillate more slowly.  In the (centered) FT, low spatial frequencies
%% are located toward the center of the spectrum; high spatial frequencies
%% are toward the edge.  Picking just the cos-phase example above, and
%% plotting only the real part of the Fourier transform, investigate this:

clf
dx = 1; xmax = 64;
x = -xmax:dx:(xmax-dx); % make a vector to represent x coordinates
cosExample1 = cos(2*pi*x/6.4);
cosExample2 = cos(2*pi*x/12.8);
cosExample3 = cos(2*pi*x/25.6);
subplot(2,3,1); plot(x,cosExample1); axis tight; title('20 cycles per window')
subplot(2,3,2); plot(x,cosExample2); axis tight; title('10 cycles per window')
subplot(2,3,3); plot(x,cosExample3); axis tight; title('5 cycles per window')

% now make a vector to represent to represent frequency - instead of
% "cycles per window", use "cycles per unit length"
fmin = -0.5/dx; fmax = 0.5/dx; df = 0.5/xmax; % a little fudging to deal with 
                                              % centering a function on an
                                              % even number of points
freq = fmin:df:(fmax-df);
subplot(2,3,4); plot(freq,real(fftshift(fft(fftshift(cosExample1)))),'r'); 
  xlabel('cycles per unit length'); title('FT of even-symmetric function')
  axis tight; legend('real')
subplot(2,3,5); plot(freq,real(fftshift(fft(fftshift(cosExample2)))),'r');
  xlabel('cycles per unit length'); axis tight
subplot(2,3,6); plot(freq,real(fftshift(fft(fftshift(cosExample3)))),'r');
  xlabel('cycles per unit length'); axis tight
  
%% change dx to .1 in the above example and change the denominator in
%% cosExample1 to 1, then execute those lines again.  This is the frequency
%% resolution limit of the function - if you zoom in on the plot, or look
%% at the numbers in cosExample1, you can see that the values are
%% alternating from -1 to 1 ... one cycle every unit of x.  

%% Now change the demoninator to something less than 1.  Where does the
%% peak appear now in the Fourier transform?  It should appear at a
%% relatively low spatial frequency (less than 1 cycle per unit length),
%% even though the actual oscillation frequency is much higher than one
%% cycle per unit length.  This is an example of aliasing - the peak is 
%% actually outside the window, but it "reflects" (is represented at the
%% wrong frequency) into the window. To fix this, change dx to something 
%% smaller (e.g. if the denominator in cosExample1 is 0.5, then dx has to 
%% be 0.25 or less).
%%
%% Note that if the denominator in cosExample is exactly 0.5, then there is
%% no oscillation - we're sampling the function only when it's high - this
%% is an extreme case of aliasing!  The only value in the FT should
%% therefore be in the center - at the DC component.
  
%% Going back to the line taken out of the zebra image, we can use the
%% above information to scale our spatial frequency axis into
%% a more useful metric, like cycles per mm

%% Example 1: Image is 10 cm wide

xmax = 10/2; % since we're going from -xmax to xmax
dx = 2*xmax/length(line); % width of a single pixel

x = -xmax:dx:(xmax-dx); % always that problem of including 0

deltaF = 0.5/xmax; fmax = 0.5/dx; % units are now cyc/cm
freq = -fmax:deltaF:(fmax-deltaF);

subplot(1,2,1)
line = line - mean(line(:)); % take the mean out
plot(x,line); % a scan through the zebra image
xlabel('distance (cm)'); ylabel('pixel intensity (dif. from mean)')

subplot(1,2,2);
plot(freq,abs(fftshift(fft(fftshift(line)))))
xlabel('cycles per mm'); ylabel('amplitude'); axis tight
title('if image is 10 mm wide')

%% There is a group of peaks on the left side of the line that, if you zoom
%% in on the plot, look like there are about 7 or 8 peaks per cm.  This
%% corresponds to the messy peak at ~7.5 cyc/cm in the FT.  The entire left
%% side of the line also seems to be modulated so about 3 clusters of peaks
%% appear in half-an-image - this corresponds to a messy peak at low
%% spatial frequencies - 3 cycles per 5cm, or 0.6 cyc/cm.  Where the messy
%% peaks around 4 cyc/cm come from is anybody's guess ...

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Two dimensional Fourier transforms %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f = fftshift(fft2(fftshift(zebras-mean(zebras(:)))));
subplot(1,2,1)
imagesc(abs(f))

%% I pick this image because the orientation bias is obvious -
%% diagonal stripes are dominant on the zebras, so the 
%% Fourier transform shows the same diagonal bias.
%% The DC value is in the center of the matrix (129,129).
%% The x coordinate of a given point tells about the projection
%% of the two-dimensional spatial frequency component (grating)
%% onto the x-axis; y coordinates show y projections. Low
%% spatial frequencies are toward the center, high spatial
%% frequencies are on the edge

%% Without using fftshift, the low spatial frequency values are
%% on the edges. 
%% The first row represents the perfectly vertical gratings 
%% (horizontal spatial frequency components); the first column
%% represents the perfectly horizontal gratings (vertical
%% spatial frequency components).  The very first value in the 
%% matrix, the first row and first column, is the DC component.
%% It should be zero if we subtract the mean before doing the FT. 
subplot(1,2,2)
F = fft2(zebras-mean(zebras(:)));
F(1,1)
imagesc(abs(F))

%% One concept we use frequently is the idea of projecting an image in one
%% dimension.  This is the same as A) taking the average of the image along
%% either the row or column dimentsion, or B) looking only at the perfectly
%% horizontal or perfectly vertical components of the Fourier transform.

%% Try killing everything but horizontal frequency components in the FT.
subplot(2,2,1)
temp = F; temp(2:end,:) = 0;
imagesc(real(ifft2(temp))); axis image
  %% the inverse Fourier transform function is the
  %% same as the forward Fourier transform function except for a scaling
  %% constant

%% Try killing everything but vertical
subplot(2,2,3)
temp = F; temp(:,2:end) = 0;imagesc(real(ifft2(temp))); axis image

%% These should look the same as spatial averaging along
%% each dimension
subplot(2,2,2)
imagesc(mean(zebras,1))

subplot(2,2,4)
imagesc(repmat(mean(zebras,2),[1 256])); axis image
  %% here I use Matlab's repmat function to turn the vector into a matrix
  %% so I can display it in a similar shape ...
