qNoise: A generator of non-Gaussian colored noise

TL;DR

qNoise is an open-source C++ library that generates customizable non-Gaussian, colored noise based on Tsallis' $q$-distribution, enabling diverse applications in modeling and analyzing stochastic systems with adjustable statistical properties.

Contribution

This paper introduces a novel C++ generator for non-Gaussian colored noise based on Tsallis' $q$-distribution, with adjustable parameters for broad applicability.

Findings

Noise statistics can be tuned via parameters $q$ and $ au$

The generator produces noise with compact support for $q<1$ and finite variance up to $q=5/3$

System responses are optimized at certain $q$ values different from 1

Abstract

We introduce a software generator for a class of colored (self-correlated) and non-Gaussian noise, whose statistics and spectrum depend on two parameters, $q$ and $\tau$. Inspired by Tsallis' nonextensive formulation of statistical physics, the so-called $q$-distribution is a handy source of self-correlated noise for a large range of applications. The $q$-noise -- which tends smoothly for $q=1$ to Ornstein--Uhlenbeck noise with autocorrelation $\tau$ -- is generated via a stochastic differential equation, using the Heun method (a second order Runge--Kutta type integration scheme). The algorithm is implemented as a stand-alone library in C++, and is made available as open source in the Github repository. Noise' statistics can be specified handily; by only varying parameter $q$: it has compact support for $q<1$ (sub-Gaussian regime) and finite variance up to $q=5/3$ (supra-Gaussian regime). Once $q$ is fixed, noise' autocorrelation can be tuned independently by means of parameter $\tau$. The presented qNoise generator provides a readily tool to modeling wide range of real-world noise types, and is suitable to study the effects of correlation and deviations from the normal distribution in systems of stochastic differential equations, key to understand system dynamics in numerous applications. The effect of noises' statistics on the response of a range of nonlinear systems is briefly discussed. In many of these examples, the systems' response turns optimal for some $q\neq1$. Hence, this paper aims to introduce qNoise generator for C++ at the class level and evaluate the kind of noise it generates, alongside their use in a range of applications.

Equations

$$\tau \dot{\eta }=-\eta (t)+\xi (t),$$

$$⟨\xi (t)⟩=0,⟨\xi (t)\xi (t^{\prime })⟩=2D\delta (t-t^{\prime }).$$

$$\tau \dot{\eta }=-\frac{d}{d\eta }{V}_q(\eta )+\xi (t)$$

$$V_q(\eta )=\frac{D}{\tau (q-1)}\ln \left[1+\frac{\tau (q-1)}{D}\frac{{\eta }^2}{2}\right],$$

$$P_q^{\mathrm{st}}(\eta )=\frac{1}{Z_q}{\left[1+\frac{\tau (q-1)}{D}\frac{{\eta }^2}{2}\right]}^{\frac{1}{1-q}},$$

$$⟨{\eta }^2⟩=\frac{2D}{\tau (5-3q)},$$

$$|\eta |<{\eta }_c\equiv \sqrt{\frac{2D}{\tau (1-q)}}.$$

$${\tau }_q\approx \frac{2\tau }{5-3q}.$$

$$⟨{\xi }_{OU}(t){\xi }_{OU}(t^{\prime })⟩=\frac{D}{\tau }\exp \left(-\frac{|t-t^{\prime }|}{\tau }\right),$$

$$P^{\mathrm{st}}({\xi }_{OU})=Z^{-1}\exp \left(-\frac{\tau }{D}\frac{{\xi }_{OU}^2}{2}\right).$$

$$m\frac{d^{2}x}{d{t}^2}=-\gamma \frac{dx}{dt}-V^{\prime }(x)-F+\xi (t)+\eta (t),$$

Full text

††journal: SoftwareX

1 Motivation

Most studies on noise-induced phenomena [1, 2] have assumed the noise source to have Gaussian distribution, either “white” (memoryless) or “colored” (red, pink, etc …). Although customarily accepted without criticism on the basis of the central limit theorem, the true rationale behind this assumption lies in the possibility of obtaining some analytical results, and avoiding many difficulties arising in generating and handling non-Gaussian noise. There is however experimental evidence that at least in some cases (particularly in sensory and biological systems) non-Gaussian noise sources may add desirable features to noise-induced phenomena (e.g. robustness, fault tolerance [3]). These findings add practical interest to the task of finding viable ways to deal with non-Gaussian noise.

Here we introduce a lightweight (generic C++ class) generator for non-Gaussian, colored stochastic processes. The expected applications of this algorithm are as diverse as the modeling of some types of vibration or fluctuation which are typically non-Gaussian, the generation of noise which is naturally confined to a domain, or the investigation of the response of many dynamical systems embedded in noise, as the latter deviates from being Gaussian.

The main features of noise obeying Tsallis’ statistics are summarized in Sec. 2; section 3 is devoted to the description of the software architecture and properties; section 4, provides statistical analysis of the generated noise in the qualitatively different cases, alongside a discussion on the $q$-dependence of the variance and effective self-correlation time. The measured self-correlation times of the obtained series are compared with a fitting expression [4]. In Sec. 5, a brief review is provided on related work, namely where non-Gaussian noise-induced phenomena have been studied.

2 $q$-noise with Tsallis’ statistics

The exponentially self-correlated Gaussian noise $\eta(t)$ named after Ornstein and Uhlenbeck (OU noise, or “colored” Gaussian noise) can be dynamically generated through the differential equation

[TABLE]

where $\xi(t)$ is centered Gaussian white noise with variance $D$, namely

[TABLE]

This way, $\eta$’s self-correlation time is $\tau$.

A straightforward generalization of Eq. (1) was proposed by Borland some time ago [5] as a model for correlated diffusion:

[TABLE]

where the potential $V_{q}$ is given by:

[TABLE]

As much as the OU noise allows to explore spectral effects within the class of exponentially correlated noise, this generalization provides moreover a device to explore statistics effects by varying just one parameter (namely $q$, at constant $\tau$ and $D$).

The stationary properties of $\eta$ (including its autocorrelation function) are thoroughly described elsewhere [4, 6, 7, 8, 9], we here summarize the main results. Using the Fokker–Planck formalism, one obtains the stationary probability distribution

[TABLE]

which can be normalized only for $q<3$ ($Z_{q}$ is a normalization factor). The first moment $\langle\eta\rangle$ always vanishes [4, 6, 7, 8, 9] and the second moment,

[TABLE]

is finite only for $q<5/3$.

Some properties of the noise are summarized in Fig. 1. The bold line depicts the Gaussian limit ($q=1$). Curves of weaker full lines show that for $q>1$, the second moment is larger than the Gaussian limit $D/\tau$. For $q<1$ (dotted lines) the distribution has a cut-off and is only defined for

[TABLE]

Same distributions are shown in linear and semilogarithmic scales (Fig. 1, left and right panels respectively).

The autocorrelation time $\tau_{q}$ of the process $\eta(t)$ in its stationary regime, also diverges for $\to 5/3\approx 1.66$. Far from its divergence point, it can be approximated as in [4]:

[TABLE]

When $q\to 1$, $\eta$ becomes a Gaussian colored noise, namely the Ornstein–Uhlenbeck process $\xi_{OU}(t)$, with correlations

[TABLE]

and probability distribution

[TABLE]

2.1 $\tau$ and $\langle\eta^{2}\rangle$ dependence on $q$

Equations Eqs. (5) and (7) tell us that for $q\neq 1$, $\langle\eta^{2}\rangle$ and $\tau_{q}$ do not attain their values ($D$ and $\tau$ respectively) in a normal distribution. Rather, they both diverge at $q=5/3$ (white squares in both panels of Fig. 2). It is however desirable to have a generator able to approximately keep constant the characteristics of these properties with respect of $q$, at least sufficiently far away from the divergence point. This can be very useful to study the effects of the statistics due to changes in $q$ keeping $\tau$ and variance constant.

This way an effective $\tau_{q}$ and $\langle\eta^{2}\rangle$ can be defined by dividing $\tau$ by Eq. (7) before integration and $\langle\eta^{2}\rangle$ by Eq. (5) after integration. The filled circles in both panels of Figure 2 show this dependence for both $\tau$ and $\langle\eta^{2}\rangle$, and how the system becomes independent of $q$ for the range $0<q<1.5$ approximately. For $q>1.5$, the proximity to the divergence point $q=5/3$ (shown with a dotted line) makes this approximation fail.

3 Generator <class> Description

The noise generator [10] is implemented in C++ as class, with dependencies on standard libraries only. It generates random numbers using functions in the built-in <random> class. The generator provides functions for Gaussian white noise, Gaussian colored noise (Orstein-Uhlenbeck), and two versions of non-Gaussian non-white noise. One where $\tau$ and $\langle\eta^{2}\rangle$ depend on $q$ (as in Eqs. 5 and 7) and a normalized version where this effect has been counterbalanced to the first order, sufficiently far away from $q=5/3$ (as shown as black dots in figure 2). This batch of functions would facilitate modeling a wide variety of scenarios and is suitable for many applications, some of which are detailed in the last section of this paper. By default, the tool uses the Mersenne-Twister generator [11] which provides a very long ($2^{19937}-1$) pseudo-random number cycle. Hence it is advised to seed the generator only once to avoid spurious correlations.

3.1 Functions

The class implements four public member functions as shown in Fig. 3.

The first function is a wrapper for the normal distribution, implemented in the <random> standard library. It is presented as a function of this class for convenience.

The second function is an implementation of the Orstein-Uhlenbeck noise. It accepts three parameters. The previous value of the noise (since it is a Markov process), the autocorrelation time $\tau$ of the noise and the integration time $H$ (necessary for setting the adequate timescale of the noise).

The third function implements the q-noise distribution. It accepts the same variables as the orsUhl function in addition to $q$ (the noise statistics), and sqrt_H as an optional variable. If $H$ is constant, explicitly setting sqrt_H = $H^{1/2}$ will avoid its calculation every time the function is called. A snippet of the function is shown in Listing 1 below.

Finally, the fourth function is a wrapper for the third function. Here $\tau$ is given by Equation (7) and the resulting noise is divided using Eq. (5) in order to counterbalance the dependence of both $\tau$ and the variance of the noise on $q$. See section 2.1 for an analysis of this effect and a discussion about its range of validity.

A generic unit test results are shown in Fig. 4. The test compares the generated average noise, of an ensemble of 10 qNoise runs for each set of $q$, $\tau$ and $N$, with the expected distribution of noise. As expected, only when $N$ is relatively small does the generated noise deviate from its theoretical distribution, particularly for high $\tau$. That is, it takes longer (higher $N$) for a highly correlated noise (high $\tau$) to explore the support and approximate the PDF.

3.2 Seeding

The presented class enables seeding the random number generator in two functions:

Timer seeding is provided as the default setting for random number generator, and it is done automatically. For lightweight single-threaded runs, manual seeding is not required. However, in multi-threading settings manual seeding (call seedManual) for each thread is recommended.

4 Properties of the generated noise

Bounded domain ($q<1$)

Bounded-domain noise is widespread in nature, and has multiple applications for modeling and control [9]111In practice, physical noise has bounded domain because arbitrarily large fluctuations are strongly suppressed. Nonetheless, Gaussian noise has many desirable theoretic properties which allow for analytical results.. The infra-Gaussian noise considered here can be addressed as a small deviation from Gaussianity, allowing a perturbative approach (Fig 5). In Sec. 5.1, an example of a infra-Gaussian noise is shown, in a resonant trap. Another use is as a source of noise whose distribution is quasi-normal but identically zero outside a boundary.

The noise generator algorithm does also ensure that noise domain is bounded, checking for out-of-bound values. This necessary test (especially for highly correlated noise) is implemented and documented accordingly in the provided source code.

Gaussian case ($q=1$)

The Gaussian case behaves exactly as an Orstein-Uhlenbeck noise, concurring perfectly with it for the whole range of $\tau$ (Fig 6). As shown in the introduction, the limit $q\rightarrow 1$ recovers the Gaussian noise, and all limits converge to it, see Eqs. (5)-(9). This limit allows to explore regions arbitrarily near the normal distribution. It can be used to model small deviations from it due to some underlying physical phenomenon. As the value of $q$ can be changed continuously and dynamically, this scheme also allows to model departures from the normal distribution due to long time-scale fluctuations, by slowly varying $1-\epsilon<q<1+\epsilon$ as a more realistic model for a small noisy system.

It is not generally recommended to compute the purely Gaussian case from the general case and set $q=1$, particularly due to potential computational complexity. An extensive batch of tests has been run in order to compare it to the Orstein-Uhlenbeck noise. All of the results were successfully recovered. As presented above, orsUhl, a function for generating Orstein-Uhlenbeck noise for a variable $\tau$ is included in this package and its results are equivalent to using the non-normalized qNoise function for $q=1$ at a fraction of the computation time.

Supra-Gaussian noise ($1<q<5/3$)

The supra-Gaussian (also called fat-tail) noise presented here is of the class of finite variance. This is usually an overlooked, modestly studied, class of noise. The Supra-Gaussian noise, generally considered in literature, tends to be Lévy-like, where the variance is infinite 222The q-noise presents infinite variance for $q>5/3$ but the description of this behavior is outside the scope of this article..

The noise presented here (Fig 7) is of a finite variance. The long excursions are however much longer and much more frequent than in the Gaussian case. This case is the most commonly used in the applications of non-Gaussian noise presented below, as it allows to model many realistic systems outside of equilibrium .

5 Applications

Some applications of the generated noise are shown below, both as an illustration and to show the potential of the generator both for research as well as a source of high quality noise for simulations or experiments (e.g. in electronics, optics, photosensitive chemical reactions etc…).

Stochastic resonance

This is a phenomenon occurring in some nonlinear systems, whereby enhancing the response to a weak external signal may require increasing the noise intensity. An often resorted-to measure is the signal-to-noise ratio at the input frequency $\omega$ (denoted by $R$).

The main numerical and theoretical results are [12, 13]: (1) for fixed $\tau$, the maximum $R$ increases with decreasing $q$; (2) for given $q$, the optimal noise intensity (the one maximizing $R$) decreases with $q$ and its value is approximately independent of $\tau$; (3) for fixed noise intensity, the optimal value of $q$ is independent of $\tau$ and in general turns out to be $q_{\mathrm{op}}\neq 1$. A simple stochastic resonance experiment with a non-Gaussian white noise [3] confirmed most of these predictions.

Brownian motors

A class of non-equilibrium systems with both potential technological applications and biological interest are the so called “ratchets”, in which the breakdown of spatial and/or temporal symmetry induces directional transport. Their transport properties can be studied by means of the Langevin equation

[TABLE]

with $m$ the particle’s mass, $\gamma$ the friction constant, $V(x)$ the (sawtooth-like) ratchet potential, $F$ a constant “load” force, and $\xi(t)$ the thermal noise, satisfying $\langle\xi(t)\xi(t^{\prime})\rangle=2\gamma T\delta(t-t^{\prime})$.

The system is kept out of thermal equilibrium by the time-correlated forcing $\eta(t)$ (with zero mean), allowing to rectify the motion. The $q$-dependence of the usual measures of performance has been studied: the mean current $J\equiv\langle dx/dt\rangle$ and the efficiency $\varepsilon$ (the ratio of the work per unit time done against $F$, to the mean power injected by $\eta$).

In the overdamped regime $(m=0,\gamma=1)$, $J$ is found to grow monotonically with $q$ whereas $\varepsilon$ is maximized for some $1<q<5/3$. For $m\neq 0$, ratchets exhibit mass-separation capabilities which are enhanced by non-Gaussian noise [14, 15]. In [16], effects of biological and technological relevance have been found in a model for the transport properties of motor proteins when departing from Gaussian behavior: $J$ is maximized not only by an optimal noise intensity but also by an optimal $q\neq 1$.

5.1 Resonant gated trapping

Stochastic resonance, which is essentially a threshold phenomenon, plays also a relevant role in ionic transport through cell membranes. In [17], a “toy model” considering the simultaneous action of a deterministic and a stochastic external field on the trapping rate of a gated imperfect trap, was studied by assuming Tsallis’ noise with $q<1$: the bounded character of the PDF contributed positively to the rate of overcoming the threshold, and such rate remained at about the same order within a larger range of values than if $\eta$ had been a white noise.

Noise-induced transition

A genetic model exhibiting a re-entrance from a disordered state to an ordered one, and again to a disordered state as $\tau$ varies from [math] to $\infty$ showed moreover a strong shift in the transition line, as $q$ departed from $q=1$. The transition was anticipated for $q>1$, while it was retarded for $q<1$ [18].

Noise-induced phase transition

In fact fat-tail noise distributions ($q>1$) counteract the effect of self-correlation (namely, they advance the ordering boundary as $D$ is increased at constant coupling), and compact-support ones ($q<1$) enhance it (they retard the ordering boundary). Particular interest rises the effect of ($q<1$) multiplicative noise on the susceptibility: it shifts from being larger on the ordering boundary to being larger on the disordering boundary [19, 20].

An example of this phenomenon can be found in climate change. Many climatic “Tipping points” are are, in fact, noise-induced phase transitions whose forcing them (including astronomic, natural and antropogenic noise) are not necessarily Gaussian. An “Early Warning”[Lenton2011] system of tipping points should include simulation considering the non-gaussianity of the stochastic forcing.

Broad-spectrum energy harvesting

In piezoelectric energy harvesting from noise, a system obeying a square-well potential can strongly profit from the large correlated excursion occurring for $q>1$ [21].

6 Conclusions

A lightweight software is presented that generates a class of non-Gaussian, colored noise. This noise can be handily generated during numerical experiments, or fed to experiments via an interface. The software, alongside documentation, is provided on the online repository Github including examples and unit test results, with an open source license. Instances of noise-induced phenomena arising when the system is submitted to (colored and non-Gaussian) noise sources with Tsallis’ q-statistics, as applications, have been briefly explored. The above discussed results show that non-Gaussian noise can significantly change the system’s response in many noise-induced phenomena, as compared with the Gaussian case. Moreover, in all the cases presented here, the system’s response was either enhanced or altered in a relevant way for values of $q$ departing from Gaussian behavior. In other words, the optimum response occurred for $q\neq 1$. Clearly, the study of the change in the response of other related noise-induced phenomena when subjected to such kind of non-Gaussian noise will be of great interest.

Acknowledgments

The authors express their recognition to H.S. Wio and R. Toral, as sources of inspiration and fruitful conversations.