Generate time-domain random signal from PSD

Given an analytical description of the PSD, for example (MATLAB "pseudocode"):

N = 20400;
bw = 2 * 2040 / N; % double-sided spectrum, 20% occupied
PSD = zeros(N,1);
P0 = sigma * sqrt(N/bw); % sigma^2 = (P0^2/N) * bw
PSD(1:2040) = P0;
PSD(18362:end) = P0;

what is the "correct" way to generate a time-domain signal $x[n]$ which is a realization of this PSD, i.e. which has variance of $\sigma^2$ and a DFT where all bins $2041, \dots, 18361$ are zero?

I found this: http://www.mathworks.com/matlabcentral/newsreader/view_thread/264846. However, I am confused since there only the phase is taken randomly but not the magnitude which would result in a perfect "brick-wall" spectrum when taking the DFT of each realization (because the magnitude and magnitude-squared is not random). However, this is only the case if the number of samples approaches infinity.

Important: I do not want to generate iid. white Gaussian samples and filter them - I would like to generate the random signal in DFT domain and then use the inverse DFT to generate the time domain signal.