Sampling Noise Given Power Spectral Density

Question

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?

Answer 1

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$.