I am trying to write a Monte Carlo physics simulation which involves, given a power spectral density, sampling rate, and total number of samples, generating noise with such a power spectral density. I have some lecture notes which work all of this out in the continuum, but I am having issues translating this to the digital domain.

The continuum formulas are as follows: we have a stochastic noise signal $$A(t) = \int_{-\infty}^{\infty} f(\omega)e^{i\omega t} \frac{d\omega}{\sqrt{2\pi}},$$ where for all $\omega \in \mathbb{R}$, $f(\omega)$ is a Gaussian random variable that satisfies $\langle f(\omega)f(\omega') \rangle = 2\pi S(\omega) \delta(\omega + \omega')$. Here $S(\omega)$ is the power spectral density of the signal. My lecture notes go on to derive using stochastic calculus the average response of certain quantum systems subject to this noise.

I am trying to replicate these results via a Monte Carlo simulation. I've tried to discretize the formulas above as best as possible (with help from this paper), and seem to be getting qualitatively correct results, but am having issues with the scale of my signal. In particular, the power of my noise signal seems to depend strongly on the sampling rate and number of samples I choose, which is obviously not desirable.

So my question is basically: assuming I correctly discretized the above formulas, is there any reason I should still expect a dependence on the sampling rate and number of samples? I would expect there to be some sort of asymptotic behavior, where as the sampling rate gets large enough the signal approaches the continuum limit, but right now I'm getting a quadratic increase in the power of the signal with higher and higher sampling rates (holding number of samples constant).

More generally, I'm curious as to why there doesn't seem to be a lot of information on this problem on the internet. In my mind, sampling noise of a given power spectral density seems like a very natural problem that would be useful in a variety of settings, and yet I have been able to find almost no information on how to do it (and very little code as well). Why is this not a problem that gets much attention?

  • $\begingroup$ Not sure I can follow what your asking. Do you just want to create time domain noise signal with a given Power Spectral Density (which is easy to do) or is there something more to your question ? $\endgroup$
    – Hilmar
    Jun 6 at 13:17
  • $\begingroup$ @Hilmar I need to generate a set of time domain signals with a given power spectral density. And in order to get the correct results from my Monte Carlo simulation, I need to generate these signals with the correct probability densities from the set of all sample paths (according to the rule I gave in the above question). $\endgroup$ Jun 6 at 16:17
  • $\begingroup$ Thanks for the paper you linked! I generate the pseudo-random noises like this. See the other answers there as well. $\endgroup$
    – Ed V
    Jul 7 at 15:24

Sampling any signal including a signal of a given power spectral density will create a unique spectrum over the range of the sampling rate (so any $n(f: f + F_s)$ for n integer will be unique in the sampled spectrum where $F_s$ is the sampling rate and the signal is assumed to be complex. Typically we use ranges of $0:F_s$ or $-F_s/2: F_s/2$ and refer to this as the "First Nyquist zone". Any continuous-time spectrum that is outside this range will map in the sampling process to this range, which includes the aliasing effects when a signal extends beyond this unique range, or what would appear as frequency translation when the signal is not in the first Nyquist zone (under-sampling).


That said, to the extent the original power spectral density is completely contained in the unique spectral range (and to the extent we know which Nyquist zone the signal occupies in continuous-time if not in the first Nyquist zone directly), we can achieve a one to one mapping from the continuous to digital domain, which means the spectral shape and total power will be conserved. When the signal's power spectral density goes beyond a single Nyquist zone, aliasing will occur destroying our ability to assess the shape of the power spectral density, but the total noise power will be maintained: consider the simple case of sampling a Gaussian white noise process where each sample is completely independent of the next (regardless of sampling rate), given enough samples we will experimentally still achieve the same histogram with the same standard deviation meaning the power is completely independent of the sampling rate). For any process in general we can see that as the sampling rate gets smaller in proportion to the frequency distribution of the signal, the sampled signals due to that can be considered independent (white in frequency), due to the frequency spread of the power spectral density combined with the aliasing effect.

The total power should be conserved assuming the sampling clock and the noise process are completely independent, but in order for the spectral features to be preserved between the continuous and discrete time domains, the noise process must either be band-limited to less than the sampling rate (for the case of a complex signal and sampling process), or the power spectral density must be periodic in frequency at intervals of the sampling rate (essentially matching what the sampled spectrum will become).

Given the OP is seeing a total power dependence on the sampling rate that continues quadratically, I suspect that the required scaling by $\Delta T$ (the sample period) when mapping from continuous to discrete domains was excluded: as we see in the mapping techniques such as in the Euler and Tustin methods, but generally, integration operations must be multiplied by $\Delta T$ (as the $dt$ goes to $\Delta T$ when the integration goes to summation) and similarly derivatives are divided by $\Delta T$ as $dx/dt$ goes to $\Delta x/ \Delta T$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.