# Sampling Noise Given Power Spectral Density

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?

• 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 ? Jun 6 at 13:17
• @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). Jun 6 at 16:17
• Thanks for the paper you linked! I generate the pseudo-random noises like this. See the other answers there as well.
– 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). 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$$.