N = 4;
M = N*2;

% 1) obtaining the Fourier matrix: 
%    for illustration only, since grossly inefficient, better use fft


w = cos(2*pi/M)  + sin(2*pi/M)*1i;     % root of unity
F = zeros(M,M);                       % matrix of a Fourier transform
f1 = zeros(M,1);
f1(1) = 1;
for i=2:M,
  f1(i) = w*f1(i-1);     
end    
F(:,1) = f1;

for i=2:M,
  F(:,i) = w*F(:,i-1);     
end  
  
imagesc([real(F) imag(F)])
 title 'Fourier matrix F, real and imaginary parts';
 colormap gray
  hold on
  plot([0 2*M],[N+0.5 N+0.5],'w', 'LineWidth',3);         % split into 4 blocks
    plot([N+0.5 N+0.5], [0 M+1],'w', 'LineWidth',3); 
    plot([N+M+0.5 N+M+0.5], [0 M+1],'w', 'LineWidth',3); 
    plot([M+0.5 M+0.5], [0 M+1],'w', 'LineWidth',3); 
  hold off
  
  % note that F11 = F22, and F12 = F21 = -F11 
  
  
% 2) illustration of fft speed and importance of M being a product of small
% primes 
    
n = 10;
M = 2^17-1       % 2^17 -1 happens to be a Mersenne prime
tic
 for i=1:n,
   Y = fft(1:M);       
 end
toc

% 3)
% illustration of spectral/fft method: for simulation on a regular grid:
% Gaussian covariance 

N = 16000;
M = N*2;

scale = 500;
c = zeros(M,1);
c(1:N+1) = exp(-((0:N)/scale).^2);

c(M:-1:N+2) = c(2:N);  % symmetrizing
plot(c)

f = fft(c);          
min(real(f)) 
max(imag(f))    %f should be real but might have a small imaginary part due to rounding 
f = max(real(f),0);           % f should be nonnegative for large enough N, 
                        % however make sure that rounding errors did not interfere 
sig = f(1:N+1)*N;
sig(1) = 2*sig(1);   sig(N+1) = 2*sig(N+1);
X = randn(N+1,1).*sqrt(sig); 
X(M:-1:N+2) = X(2:N); 
Y = randn(N+1,1).*sqrt(sig); 
Y(M:-1:N+2) = -Y(2:N);            
Y(1) =0;
Z = X + 1i*Y;
V = ifft(Z);      

 plot(real(V))        % V has a small imaginary part due to roundoff


 % 4) 
 % illustration of 2-dim FFT and its use in image compression
 
 X = double(imread('./data/face.bmp'));
 X = X(:,:,1);        % this file is originally an RGB bitmap, take only one component
 % imwrite(X,'face1.bmp');
 imagesc(X)
  colormap gray
 sz =  size(X)     
 nrow = sz(1);  ncol = sz(2);
  
 f = fft2(X);
 imagesc(-[real(f) imag(f)])
   colormap gray
 
  f(2,1:5)   % shows symmetries
  f(96,1:5)
   
  % for compressed image: delete some Fourier coefficients 

  f = fft2(X); 
  Kx = 40;   % how many coefficients to keep for each direction
  Ky = 30;
  f(Ky+1:nrow - Ky+1, :) = 0;    
  f(:, Kx + 1: ncol-Kx+1) = 0;      
  
  compratio = Kx*Ky/(nrow*ncol)  
  imagesc([X real(ifft2(f))]);
    title 'True and compressed images';
     colormap gray