Noise-to-signal ratio of single-trajectory spectral densities in centered Gaussian processes

We discuss the statistical properties of a single-trajectory power spectral density $S(\omega,\mathcal{T})$ of an arbitrary real-valued centered Gaussian process $X(t)$, where $\omega$ is the angular frequency and $\mathcal{T}$ the observation time. We derive a double-sided inequality for its noise-to-signal ratio and obtain the full probability density function of $S(\omega,\mathcal{T})$. Our findings imply that the fluctuations of $S(\omega,\mathcal{T})$ exceed its average value $\mu(\omega,\mathcal{T})$. This implies that using $\mu(\omega,\mathcal{T})$ to describe the behavior of these processes can be problematic. We finally evaluate the typical behavior of $S(\omega,\mathcal{T})$ and find that it deviates markedly from the average $\mu(\omega,\mathcal{T})$ in most cases.


Introduction
The power spectral density (PSD) of a deterministic or stochastic process X(t) encodes important information about its properties and is widely used in experimental, numerical and theoretical analyses (see, e.g., refs. [1][2][3][4] and references therein). The singletrajectory PSD S(ω, T ) is defined by where ω is the (angular) frequency and T the observation time. Usually, one averages S(ω, T ) over an ensemble of trajectories and eventually takes the T → ∞ limit to get where the overbar here and henceforth denotes the ensemble averaging. We emphasise that µ(ω) in eq. (2)-an ensemble-averaged property taken in the limit T → ∞-is conventionally referred to as the PSD and its calculation is the usual target of the standard analyses of spectral properties of random processes.
Let us remark that the T → ∞ can be formally taken in mathematical expressions but not in experimental or numerical analyses, and that therefore caution is required when making comparisons with theoretical predictions. Moreover, for many nonstationary stochastic processes the infinite-T limit of the expressions in eqs. (1) and (2) does not exist, what requires either using some alternative definitions of power spectral densities (see, e.g., [2,5]) or to confine oneself to the finite-T behaviour. Lastly, it is not always possible to have a large enough statistical sample in order to reliably perform the averaging.
Motivated by the latter circumstance, recent works [6][7][8][9][10] have concentrated on stochastic properties of random variable S(ω, T ) defined in eq. (1). Indeed, it is important to know how S(ω, T ) fluctuates from sample to sample, in order to estimate how large a statistical sample should be to allow a reliable evaluation of µ(ω, T ) from experimental or numerical analyses, especially in view of taking its large T limit.
It turns out that, for several centered Gaussian processes-standard Brownian motion [6], fractional Brownian motion with Hurst index H [7], scaled Brownian motion [8], diffusing diffusivity processes [9] and the Brownian gyrator model [10]the probability density function P (s) = P (S(ω, T ) = s) of the random variable S(ω, T ) is explicitly given, for arbitrary values of ω and T , by the universal form where I 0 (z) is the modified Bessel function and γ is the "noise-to-signal" ratio, defined by This latter parameter, which is also called the coefficient of variation of the distribution (3), is a measure of the relative weight of fluctuations of the finite-T PSD around its mean value. Note that the form of eq. (3) requires γ to satisfy the doublesided inequality which was directly verified in [6][7][8][9][10] for each particular case under study.
Because the inequality in eq. (6) and the distribution (3) appear to be valid for rather diverse Gaussian stochastic processes, one can conjecture that they hold in general for an arbitrary centered Gaussian process. In this paper we focus on this question and present a formal proof that this is indeed the case.
In Section 2 we prove the crucial inequality (6) for an arbitrary Gaussian centered process with arbitrary ω and T . Our proof is solely based on Wick's theorem. In Section 3 we take advantage of the Karhunen-Loeve decomposition of an arbitrary Gaussian process to calculate the moment-generating function of S(ω, T ), from which the general result in eq. (3) follows by the inversion of the Laplace transform. Section 4 is devoted to the analysis of the functional forms of P (s) in two limiting cases, and also to a discussion of the typical behavior of S(ω, T ). We conclude with a brief summary of our results in Sec. 5.

The noise-to-signal inequality
The second moment S 2 (ω, T ) of a single-trajectory PSD of a real-valued process X(t) can be formally written down as According to Wick's theorem, the four-time correlation function of the form for an arbitrary centered Gaussian process naturally decomposes as follows: This implies that the expression (7) can be formally rewritten as where we have used the shorthand notations Expression (9) implies that the variance var S(ω, Our first goal is to prove that var S(ω, T ) ≥ µ 2 (ω, T ), i.e., that γ ≥ 1. To this end, we notice that and rewrite formally eq. (11) as where The lower bound follows by merely noticing that ∆(ω, T ) > 0 for any ω and T . It may be instructive to note that despite the simplicity of its derivation, this bound has strong implications: namely, it shows that fluctuations of S(ω, T ) for arbitrary centered Gaussian processes exceed generically the mean value µ(ω, T ). This implies that the knowledge of µ(ω, T ) alone may not be sufficient to fully characterise the behaviour of this random variable. It also implies that estimating µ(ω, T ) from numerical or experimental data reliably well may require quite large statistical samples.
In order to prove the upper bound γ ≤ √ 2, we have to show that ∆(ω, T ) ≤ µ 2 (ω, T ) or, equivalently, that We rewrite next the latter inequality in the explicit form and take advantage of Mercer's theorem [11,12]. This theorem asserts that for an arbitrary symmetric, continuous non-negative kernel function X(t 1 )X(t 2 ) there exists an orthonormal set of eigenfunctions e k (t) defined in [0, T ] with positive eigenvalues λ k , such that the covariance function X(t 1 )X(t 2 ) can be expressed by and that this series converges absolutely and uniformly. The eigenvalues λ k and eigenfunctions e k (t) are generally found by solving the homogeneous Fredholm integral equation of the second kind: We substitute the expansion (17) in the inequality (6) and introduce the notation We can then cast the upper bound in inequality (6) where the existence of the limits on the right-and left-hand-sides is ensured by the absolute convergence of the expansion (17). Equation (20) is the standard Cauchy-Schwarz inequality, which thus proves the upper bound on the coefficient of variation γ. Note that the inequality becomes an identity for ω = 0 when both sides vanish.

The moment-generating function
In order to evaluate the probability distribution function of S(ω, T ), we consider the moment-generating function Φ κ of the single-trajectory PSD of an arbitrary real-valued, centered Gaussian process X(t). This function is defined by (21) ‡ In particular, when X(t) is the standard Brownian motion, e k (t) ∝ sin((k − 1/2)πt/T ) and λ k = 1/(π 2 (k − 1/2) 2 ). This corresponds to the celebrated Wiener representation of the BM. For the FBM the eigenfunctions e k (t) are combinations of the sine and cosine functions as above, but the eigenvalues λ k are not simply multiples of π 2 and are expressed through the zeros of the Bessel functions J H (x) and J 1−H (x) (see [13]).

We rewrite eq. (1) in the form
The expression in the last line ensures that S(ω, T ) ≥ 0 for any ω and T .
We take advantage of the Karhunen-Loeve decomposition (see, e.g., [14]), according to which any zero-mean square-integrable Gaussian stochastic process X(t), defined on the interval [0, T ], admits the following representation: Here e k (t) are the above defined orthonormal eigenfunctions and Z k are independent, normally distributed random variables, with zero mean and variance λ k . Substituting (23) in eq. (22), we obtain where the c k and s k are defined in eqs. (19).
Further on, we use the identity which permits us to write down Φ κ as the following two-fold integral: The averaging can now be straightforwardly performed to give where the coefficient in front of κ 2 is evidently positive, by virtue of eq. (20).
The last step consists in identifying the coefficients in front of κ and κ 2 in the last line in eq. (27). Using eqs. (11) and (12), we readily obtain that By inverting the Laplace transform in eq. (27) we obtain eq. (3). We note moreover that, similarly to the parental Gaussian process X(t), the probability density function in eq. (3) is entirely defined by the first two moments of S(ω, T ).

Limiting cases and typical behaviour of the PSD
We consider here two limiting situations in which the functional form of the probability density function in eq. (3) simplifies; namely, when γ → 1 (i.e., var S(ω, T ) → µ 2 (ω, T )) or when γ → √ 2 (i.e., var S(ω, T ) → 2µ 2 (ω, T )). Note that, in general, the coefficient of variation is an oscillating function of the frequency for fixed T , but attains a constant value in the T → ∞ limit. In particular, we have γ → 1 in the limit T → ∞ for, e.g., the sub-diffusive fractional Brownian motion [7] or for the Brownian gyrator model [10]. In turn, γ → √ 2 holds in the same limit for the super-diffusive fractional Brownian motion [7]. Moreover, γ = √ 2 holds as an identity for any T in any centered Gaussian process when ω = 0, a relation that follows immediately from eqs. (13) and (14). Indeed, S(ω = 0, T ) is equal to the squared area under the random Gaussian curve X(t), divided by the observation time (see eq. (1)). The limiting forms of the probability density function can be conveniently studied using the expression (27).
When γ = 1, the expression in the last line in eq. (27) becomes a full square, and thus one has implying that the probability density function is a simple exponential. This form has been experimentally verified for the sub-diffusive fractional Brownian motion [7] and for the Brownian gyrator model [10].
When γ = √ 2, the coefficient in front of κ 2 in eq. (27) vanishes and the momentgenerating function becomes which signifies that the probability density function converges to the χ 2 -distribution with one degree of freedom, i.e., Note that for ω = 0, when the single-trajectory PSD defines the squared area under random curve X(t), this result simply states that the area itself has a Gaussian distribution, which is an a priori known result. The form in (31) has also been verified in experimental analyses of fractional Brownian motion processes [7] and for the Brownian gyrator model [10].
Since γ ≥ 1, the magnitude of fluctuations of S(ω, T ) exceeds generically its mean value. In other words, the fluctuations of S(ω, T ) over different realizations of X(t) are significant and µ(ω, T ) is most likely dominated by some atypical realisations of the random process X(t). This means, in turn, that in order to correctly reproduce the analytical predictions from numerics or experiments, the number of realizations of the process has to be large enough in order to "catch" such rare trajectories. The function in eq. (3) attains its maximal value at s = 0, so that we can expect that typical trajectories most often yield a value of S(ω, T ) smaller than the average. In our case, we estimate the typical value of S(ω, T ) as the exponential of the average of the logarithm of S(ω, T ): For the probability density function given in eq. (3) we can perform the corresponding integral exactly, and obtain ln S(ω, T ) µ(ω, T ) = ln where C ≈ 0.577 is the Euler-Mascheroni constant. (Notice that the convergence of the integral is guaranteed by the double-sided inequality (6).) Consequently, the typical value of a single-trajectory PSD is given by As intuitively expected, µ typ (ω, T ) is smaller than µ(ω, T ) and the difference between them is more pronounced for γ close to √ 2 than for γ close to 1. Therefore, in order to obtain reliable estimates of µ(ω, T ), statistical samples need to be larger in the former than in the latter case.

Conclusions
Summarizing, we have analysed here the statistical properties of the fluctuations of a single-trajectory power spectral density S(ω, T ) of centered Gaussian processes X(t), going beyond its first moment, which is the main focus of most analyses. We have presented a formal proof of the statement that for an arbitrary Gaussian process the noise-to-signal ratio γ, defined as the ratio of the standard deviation and the mean value of S(ω, T ), obeys a double-sided inequality 1 ≤ γ ≤ √ 2 for any ω and any T . The bound γ > 1 implies that the magnitude of fluctuations generically exceeds the mean value, meaning that the realisation-to-realisation fluctuations are very significant. Using the Karhunen-Loeve decomposition of X(t), we evaluated the full probability density function of S(ω, T ), which holds for arbitrary centered Gaussian processes and for any value of the frequency and of the observation time. Finally, we have discussed the typical behaviour of a single-trajectory power spectral density, which is most likely to be observed for small statistical ensembles of trajectories.