% fourmat.m, the Fourier matrix.
clear

N=4;
F=dftmtx(N)/sqrt(N); % Matlab does not make F norm-preserving (unitary)

% construct F*F as a Hankel matrix
dum=zeros((N+1),1);
dum(2)=1;
F2=hankel(dum(2:N+1),dum(1:N));
sprintf('Test on F*F %0.3g',norm(F2-F*F,'fro'))
pause

% construct Ia
Ia=hankel(flipud(dum(2:N+1)));

% constuct a general ciculant Toeplitz matrix
dum=randn((N+1),1);
dum(N+1)=dum(1);
T=toeplitz(dum(1:N),flipud(dum(2:N+1)));

% and the Hankel
H=T*Ia;

% illustrate how to diagonalize these matrises
sprintf('Test if FT*hermit(F) is diagonal %0.3g',norm(F*T*F'-diag(diag(F*T*F')),'fro'))
sprintf('Test if FHF is diagonal %0.3g',norm(F*H*F-diag(diag(F*H*F)),'fro'))
pause

% compare the diagonal of F*T*F' to the eigenvalues of T
sprintf('Now follows eigenvals of T')
sort(diag(F*T*F'))'
sort(eig(T))'
norm(sort(diag(F*T*F'))-sort(eig(T)),'fro') % dangerous as sort is imperfect
pause

% compare the diagonal of F*H*F to the eigenvalues of H
sprintf('Now follows eigenvals of H')
sort(diag(F*H*F))'
sort(eig(H))'
% Now do it the proper way
dum=diag(F*H*F);
for k=2:ceil(N/2)
   dum1=sqrt(dum(k)*dum(N-k+2));
   dum(k)=dum1;
   dum(N-k+2)=-dum1;
end
sprintf('Now follows eigenvals of H the proper way')
sort(real(dum))' % lack of numerical precision leaves a small imaginary part
norm(sort(real(dum))-sort(eig(H)),'fro') % now it works as real-valued data
