function F=fouriermat(n,fig)

% This function calculates the nxn Discrete Fourier Transform matrix F,
% i.e. if x is nx1, then F*x is the discrete fourier transform of x.   
%
% F is the symmetric nxn Vandermonde matrix 
%   F = [1 zeta zeta^2 ... zeta^(n-1)]
% where zeta is a column vector of the n complex roots of unity
% 
% If fig is set to a positive integer, then a figure window shows the real
% part of each column of F.

if ~exist('fig','var')
    fig=0;
end


% Create a discrete fourier transform matrix of size nxn

%tic
if 1
  F=ones(n);
  for I=1:(n-1)
    F(I+1,2:end)=exp(-2*(I:I:(I*(n-1)))*pi*i./n);
    
  end

elseif 0
  % Intuitive yet unstable and more expensive way to do the calculation, by 
  % explicitly taking powers
  zeta=[1; exp(-2*(1:(n-1))'.*pi*i./n)];  
  F=vandermonde(zeta);
end
%toc

if fig
  figure(fig)
  for j=1:n
    plot(real(F(j,:)))
    title(sprintf('Plot of real part of Fourier mode at frequency -2*pi*%d/%d',(j-1),n))
    fprintf('Hit any key to continue\n')
    pause
  end
end