To my knowledge, white Gaussian noise (WGN) is defined as a process with a correlation function:
$$ R[k]=\sigma^2 \delta[k] $$
and whose symbols are distributed according to $N(0,\sigma^2)$.
Naturally, when simulating such noise, the mean won't be zero for a specific realization but it would be close. However, if I try to force a zero mean (simply subtract the mean), I then get the undesirable characteristic in the frequency domain. I get that my noise process does not have a DC term (instead of having a constant spectrum).
Can I simulate WGN with zero mean and flat spectrum? Is the problem numeric, or is it my interpretation of the results which is flawed?