function A=fouriermat(N)
% fouriermat(N): return N*N real fourier basis
%
% fouriermat :: N:natural -> [[N,N]].

if rem(N,2)
	% for odd N 
	A=ones(n,1);

	t=(0:n-1)';
	for k=1:floor((n-1)/2)
		a=exp(2*i*pi*k*t/n);
		A=[A real(a) imag(a)];
	end
	if mod(n,2)==0,
		A=[A repmat([1; -1],[n/2 1])];
	end

	A=A*diag(1./sqrt(sum(A.^2)));
else
	% for even N. ought to fix this for odd N
	A=zeros(N,N);
	A(1,1)=1;
	A(N/2+1,N)=1;
	sq2=sqrt(2)/2;
	for j=2:N/2
        A(j,2*j-2)=sq2;
        A(2+N-j,2*j-2)=sq2;
        A(j,2*j-1)=sq2*i;
        A(2+N-j,2*j-1)=-sq2*i;
	end
	A=real(fft(A)/sqrt(N));
end


