samer@32: function A=fouriermat(N)
samer@32: % fouriermat(N): return N*N real fourier basis
samer@32: %
samer@32: % fouriermat :: N:natural -> [[N,N]].
samer@32: 
samer@32: if rem(N,2)
samer@32: 	% for odd N 
samer@32: 	A=ones(n,1);
samer@32: 
samer@32: 	t=(0:n-1)';
samer@32: 	for k=1:floor((n-1)/2)
samer@32: 		a=exp(2*i*pi*k*t/n);
samer@32: 		A=[A real(a) imag(a)];
samer@32: 	end
samer@32: 	if mod(n,2)==0,
samer@32: 		A=[A repmat([1; -1],[n/2 1])];
samer@32: 	end
samer@32: 
samer@32: 	A=A*diag(1./sqrt(sum(A.^2)));
samer@32: else
samer@32: 	% for even N. ought to fix this for odd N
samer@32: 	A=zeros(N,N);
samer@32: 	A(1,1)=1;
samer@32: 	A(N/2+1,N)=1;
samer@32: 	sq2=sqrt(2)/2;
samer@32: 	for j=2:N/2
samer@32:         A(j,2*j-2)=sq2;
samer@32:         A(2+N-j,2*j-2)=sq2;
samer@32:         A(j,2*j-1)=sq2*i;
samer@32:         A(2+N-j,2*j-1)=-sq2*i;
samer@32: 	end
samer@32: 	A=real(fft(A)/sqrt(N));
samer@32: end
samer@32: 
samer@32: