function F=dftfmap(N,Fs)
% dftfmap - Frequency map for DFT of real signal 
%
% dftfmap :: N:natural, real ~'sampling rate' -> dmap(dftbins(N)).

	M=dftbins(N);
	F=dmap(M,@map,@revmap);

	function I=map(X)
		I=round(N*X/Fs);
		I(X<0)=-inf;
		I(X>Fs/2)=inf;
	end

	function X=revmap(I)
		I=shiftdim(shiftdim(I),-1);
		X1=Fs*(2*I-3)/(2*N);
		X2=Fs*(2*I-1)/(2*N);
		X1(I<=1)=0;
		X2(I>=M)=Fs/2;
		X=cat(1,X1,X2);
	end
end
		
