Colored noise in oscillators. Phase-amplitude analysis and a method to avoid the Ito-Stratonovich dilemma

TL;DR

This paper analyzes how colored noise affects oscillator phase fluctuations, introduces a transformation method to handle the Ito-Stratonovich dilemma, and improves phase noise modeling accuracy.

Contribution

It presents a novel transformation technique to convert colored noise into white noise for better phase noise analysis in oscillators.

Findings

Phase noise is a drift-diffusion process influenced by noise variance and correlation time.

The method reduces modeling errors compared to previous phase reduced models.

A new approach avoids the Ito-Stratonovich dilemma in stochastic oscillator analysis.

Abstract

We investigate the effect of time-correlated noise on the phase fluctuations of nonlinear oscillators. The analysis is based on a methodology that transforms a system subject to colored noise, modeled as an Ornstein-Uhlenbeck process, into an equivalent system subject to white Gaussian noise. A description in terms of phase and amplitude deviation is given for the transformed system. Using stochastic averaging technique, the equations are reduced to a phase model that can be analyzed to characterize phase noise. We find that phase noise is a drift-diffusion process, with a noise-induced frequency shift related to the variance and to the correlation time of colored noise. The proposed approach improves the accuracy of previous phase reduced models.

Equations

$${\dot{\mathbf{x}}}_t=\mathbf{a}({\mathbf{x}}_t)+\mathbf{B}({\mathbf{x}}_t){\eta }_t$$

$$\tau {\dot{\eta }}_t=-{\eta }_t+D{\dot{W}}_t$$

$$⟨{\eta }_t⟩={\eta }_0\exp \left(-\frac{t}{\tau }\right)$$

$$⟨{\eta }_t{\eta }_s⟩=\frac{D^2}{2\tau }\exp \left(-\frac{|t-s|}{\tau }\right)$$

$$d{\mathbf{x}}_t$$

$$\tau d{\eta }_t$$

$$d{\mathbf{x}}_t=\mathbf{a}({\mathbf{x}}_t)dt+D\mathbf{B}({\mathbf{x}}_t)d{W}_t$$

$$d{\mathbf{x}}_t$$

$$d{y}_t$$

$$\frac{\partial u}{\partial t}=\left({\Lambda }_0+\frac{1}{\epsilon }{\Lambda }_1+\frac{1}{{\epsilon }^2}{\Lambda }_2\right)u$$

$${\Lambda }_{0}u$$

$${\Lambda }_{1}u$$

$${\Lambda }_{2}u$$

$${\epsilon }^{-2}:{\Lambda }_2{u}_0$$

$${\epsilon }^{-1}:{\Lambda }_2{u}_1$$

$${\epsilon }^0:{\Lambda }_2{u}_2$$

$$⋮$$

$$\psi (y)=p_{\text{st}}(y)=\sqrt{\frac{1}{\pi D^2}}\exp \left(-\frac{y^2}{D^2}\right)$$

$$(b_n,\psi )=⟨b_n{⟩}_y={\int }_{\text{mathds}R}{b}_n{p}_{\text{st}}(y)dy=0$$

$$u_1(\mathbf{x},y,t)=y\frac{\partial u_0(\mathbf{x},t)}{\partial \mathbf{x}}\mathbf{B}(\mathbf{x})=-{\Lambda }_2^{-1}{\Lambda }_1{u}_0$$

$${⟨\frac{\partial u_0}{\partial t}⟩}_y-{⟨{\Lambda }_0{u}_0⟩}_y+{⟨{\Lambda }_1{\Lambda }_2^{-1}{\Lambda }_1{u}_0⟩}_y=0$$

$${\int }_{\text{mathds}R}\frac{\partial u_0(\mathbf{x},t)}{\partial t}{p}_{\text{st}}(y)dy=\frac{\partial u_0(\mathbf{x},t)}{\partial t}$$

$${\int }_{\text{mathds}R}{\Lambda }_0{u}_0{p}_{\text{st}}(y)dy=\frac{\partial u_0(\mathbf{x},t)}{\partial \mathbf{x}}\mathbf{a}(\mathbf{x})$$

$$\displaystyle\int_{\mathds{R}}\Lambda_{1}\Lambda_{2}^{-1}\Lambda_{1}u_{0}\,p_{\text{st}}(y)\,dy=-\dfrac{D^{2}}{2}\bigg{[}\dfrac{\partial u_{0}(\boldsymbol{x},t)}{\partial\boldsymbol{x}}\,\dfrac{\partial\boldsymbol{B}(\boldsymbol{x})}{\partial\boldsymbol{x}}\,\boldsymbol{B}(\boldsymbol{x})$$

$$\displaystyle\hskip 99.58464pt+\boldsymbol{B}^{T}(\boldsymbol{x})\,\dfrac{\partial^{2}u_{0}(\boldsymbol{x},t)}{\partial\boldsymbol{x}^{2}}\,\boldsymbol{B}(\boldsymbol{x})\bigg{]}$$

$$\frac{\partial u_0(\mathbf{x},t)}{\partial t}=$$

$$\displaystyle+\boldsymbol{B}^{T}(\boldsymbol{x})\,\dfrac{\partial^{2}u_{0}(\boldsymbol{x},t)}{\partial\boldsymbol{x}^{2}}\,\boldsymbol{B}(\boldsymbol{x})\bigg{]}$$

$$d{\mathbf{x}}_t=\mathbf{a}({\mathbf{x}}_t)dt+D\mathbf{B}({\mathbf{x}}_t)\circ d{W}_t$$

$$d{\mathbf{x}}_t=\left[\mathbf{a}({\mathbf{x}}_t)+\frac{D^2}{2}\frac{\partial \mathbf{B}({\mathbf{x}}_t)}{\partial \mathbf{x}}\mathbf{B}({\mathbf{x}}_t)\right]dt+D\mathbf{B}({\mathbf{x}}_t)d{W}_t$$

$$\frac{d\mathbf{x}}{dt}=\mathbf{a}(\mathbf{x})$$

$${\mathbf{u}}_1(t)=\frac{\mathbf{a}({\mathbf{x}}_s(t))}{|\mathbf{a}({\mathbf{x}}_s(t))|}$$

$$\frac{d\tilde{\mathbf{x}}(t)}{dt}={\mathbf{J}}_{\mathbf{a}}(t)\tilde{\mathbf{x}}(t)$$

$${\mathbf{x}}_t={\mathbf{x}}_s({\theta }_t)+\mathbf{Y}({\theta }_t){\mathbf{R}}_t$$

$$d{\theta }_t=$$

$$+B_{\theta }({\theta }_t,{\mathbf{R}}_t)d{W}_t$$

$$d{\mathbf{R}}_t=$$

$$+{\hat{\mathbf{a}}}_{\mathbf{R}}({\theta }_t,{\mathbf{R}}_t)+{\mathbf{b}}_{\mathbf{R}}({\theta }_t,{\mathbf{R}}_t)]dt+$$

Full text

I Introduction

Oscillators and phase locked loops (PLLs) are fundamental components of electronic and optical systems. For instance, in digital systems they are used to establish a reference time to synchronize operations. In communication systems, they are used for frequency coding and decoding, and for channel selection.

Noise sources, both intrinsic and external, are a major nuisance plaguing oscillator and PLL performance. They can be classified as white (frequency independent) fluctuations, such as thermal noise in electrical circuits with resistive elements or shot noise in semiconductor devices, and time-correlated (colored) noise sources. Among the latter, particularly relevant in oscillators based on bipolar transistors and MOSFET devices used as radio frequency sources, we find Lorentzian low-frequency noise and flicker noise. Flicker or $1/f$ fluctuations can be in several cases traced back to the superposition of low-frequency Lorentzian noise sources, that therefore are of paramount importance in assessing oscillator random variations.

The performance and reliability of oscillators depend crucially on noise sources, which deteriorate the oscillator response and are responsible for phase noise and time jitter. Phase noise and time jitter are strictly related concepts, defining oscillator short term frequency instabilities. In particular, phase noise is a frequency domain measure of the oscillator spectral purity, while time jitter describes the time domain accuracy of the oscillator waveforms. Phase noise spreads the bandwidth of the fundamental frequency of oscillators and may produce interference with neighboring channels, thus degrading the whole system performances. Therefore, characterizing phase noise in oscillators is a major problem for practical applications.

Since the seminal work [1], linear time invariant models (LTI) have been applied to high-$Q$ resonant and quartz-crystal oscillators. While of great practical importance, such a technique is often too simplistic and fails to capture essential features such as spectral dispersion. Inclusion of linear time variant effects (LTV) can yield more accurate results [2, 3].

The most rigorous treatment of phase noise in nonlinear oscillator perhaps dates back to [4], where the author decomposes the oscillator response into phase and magnitude components, and successfully derived a differential equation for the phase deviation. In [5], the oscillator response is decomposed into orthogonal components, and equations for purely phase and amplitude deviations are derived. Unfortunately, as shown in [6], using an orthogonal decomposition to separate phase and amplitude deviations leads to inaccurate results. The method proposed in [6], based on using Floquet vectors to project the response of the noisy oscillator onto a shifted version of the unperturbed response, exploits a linear periodically time varying approximation of the oscillator behavior leading, ultimately, to a nonlinear phase equation. The idea that Floquet vectors constitute the ideal basis to decouple phase and amplitude dynamics was confirmed in [7, 8]. Phase domain models based on the ideas introduced in [6], have been extensively used to derive an analytical stochastic characterization of oscillator phase noise due to both white and colored noise sources [9, 10, 11, 12, 13, 14, 15, 16].

In [4, 6, 9], simplified scalar stochastic differential equations for the phase variable were derived, neglecting all contributions of noise and amplitude fluctuations, i.e small deviations from the limit cycle, beyond the first order. This assumption is well justified for most electronic oscillators, that are typically subject to small noise and exhibit strongly stable limit cycle, thus leading to very accurate phase models. LTV techniques yield reliable and precise predictions for the phase diffusion process, but they fail to capture the frequency shift phenomenon [11, 17, 18]. Such a frequency shift is usually negligibly small for strongly stable oscillators, such as those customarily used in electronic systems, but may play a relevant role in autonomous systems from other fields, e.g. system biology and neuroscience [19].

This paper proposes a novel methodology to analyze phase noise in nonlinear oscillators subject to time-correlated noise, modeled as an Ornstein-Uhlenebeck process. Making use of the method proposed in [20, 21] the system with colored noise is first transformed into an equivalent system subject to white Gaussian noise. The advantage is twofold: first the transformation allows to use the whole machinery of stochastic differential equations on a reduced dimensional system. Second, and more important, the transformation avoids the Itô-Stratonovich dilemma [22]. Phase and amplitude deviation equations are derived for the transformed system, and then reduced to a phase model, that describes the oscillator dynamics in terms of the phase variable only. This description is the ideal tool for phase noise analysis, since it gives an approximate yet very accurate description of the phase dynamics. We show that phase noise in nonlinear oscillators is a drift-diffusion process, that is, noise sources not only induce a spread in the oscillator spectrum, but also a shift in the oscillation mean frequency. The frequency shift is related to the variance of the colored noise and to the noise correlation time. Some examples are presented to assess the validity of the approach.

II Modeling

Consider the nonlinear system subject to random noise source

[TABLE]

where $\boldsymbol{x}_{t}:\mathds{R}\mapsto\mathds{R}^{n}$ denotes the state of the system, $\boldsymbol{a}:\mathds{R}^{n}\mapsto\mathds{R}^{n}$ is a smooth vector field that defines the system internal dynamics, $\boldsymbol{B}:\mathds{R}^{n}\mapsto\mathds{R}^{n}$ is a smooth vector valued function and $\eta_{t}:\mathds{R}\mapsto\mathds{R}$ is a scalar function describing random fluctuations, both internal and external.

Random fluctuations are often modeled as zero mean, Gaussian distributed white noise. A zero mean Gaussian white noise $\eta_{t}=\dot{W}_{t}$ is characterized by $\langle\dot{W}_{t}\rangle=0$ and $\langle\dot{W}_{t}\dot{W}_{s}\rangle=\delta(t-s)$, where $\delta$ is the Dirac delta function, while $t$ and $s$ denote two time instants. White Gaussian noise is a reasonable approximation in the case where the typical time scales of the underlying deterministic dynamics are much larger than the noise correlation time (quasi-white approximation). Unfortunately, white noise is rarely an accurate model to represent all of the noise sources that induce fluctuations in real electronic devices and systems. As a consequence, a Dirac delta correlation in time is more justified by mathematical convenience than being physically plausible. Real processes have typically finite correlation times, and often $1/f$ power spectra, e.g. flicker noise.

A more realistic description of noise in electronic systems is given by an exponentially correlated process. Better known as colored noise, it can be modeled as an Ornstein-Uhlenbeck process [23]. In this case the noise source is modeled according to

[TABLE]

where $\tau$ is a parameter proportional to the finite noise correlation time, and $D$ is the diffusion constant. The Ornstein-Uhlenbeck process with deterministic initial state $\eta_{0}$ is characterized by the expectation value

[TABLE]

and by the correlation function

[TABLE]

Thus we consider equations (1) and (2), that we rewrite adopting the standard notation of stochastic differential equations (SDEs)

[TABLE]

where $W_{t}$ is a Wiener process, i.e. the integral of a white noise.

The SDEs (5)-(6) describe a diffusion process with unmodulated (additive) noise. Therefore the equations can be interpreted using any of the two main interpretation schemes, namely Itô or Stratonovich, obtaining the same solution. In the full space $(\boldsymbol{x},\eta)$ the system is Markovian, that is, future states are completely determined by the current state and the stochastic process shows no memory. Methods for analysis of Markovian systems, based on the Fokker-Planck or Kolmogorov equations, are well developed. However, the practical solution of the equations obtained using these approaches may become unbearable because of the large number of state variables involved.

A possible solution strategy amounts to study the system dynamics in a reduced dimension space. However, if we consider the $(\boldsymbol{x})$ space only the system is non Markovian due to the presence of multiplicative noise, as the increments of the state variables depend on the past history of the noise process. Problems arise even in the simpler case of a quasi-white approximation: for $\tau\rightarrow 0$, (6) shows that the external fluctuations reduce to a white noise $\eta_{t}dt=DdW_{t}$. Substituting this approximation into (5) yields

[TABLE]

The white noise is now modulated by the state dependent function $\boldsymbol{B}(\boldsymbol{x})$ (multiplicative noise), therefore the question arises whether (7) should be interpreted according to Itô or Stratonovich.

Applying the procedure presented in [20, 21], in the next section we shall derive a reduced description in the $(\boldsymbol{x})$ space of problem (5)-(6) where the SDE system is transformed into an SDE (for the state vector $\boldsymbol{x}$ only) subject to white Gaussian noise. The main results are the following

• •

The reduced system holds for small, but not necessarily vanishing correlation time $\tau$.

• •

The reduced system resolves the Itô-Stratonovich dilemma, that is, we shall derive equivalent SDEs for the two interpretations. By equivalent, we mean two different SDEs, interpreted following different rules, having the same solution. Because the solution is unique, it is just a matter of personal preference to choose one interpretation rather than the other.

III White noise approximation

Dividing both sides of (6) by $\tau$, substituting $\tau=\varepsilon^{2}$ and introducing $\eta_{t}=y_{t}/\varepsilon$, equations (5)-(6) become

[TABLE]

Usually, in electronic systems the correlation time $\tau$ is small compared to the characteristic time constants for the deterministic part of the dynamics. Therefore, we can assume the correlation time $\tau$ small enough such that $\varepsilon\ll 1$. Under this assumption, equations (8) and (9) show a time scale separation, since the Ornstein-Uhlenbeck process $y_{t}$ is one order, in the parameter $\varepsilon$, faster than the state variables $\boldsymbol{x}_{t}$. Notice that a straightforward application of stochastic averaging [24] would lead to inconsistent results. In fact, since asymptotically $\langle y_{t}\rangle=0$, the averaged equation would simply coincide with the deterministic (noiseless) system.

Now, let $u(x,t)=\mathsf{E}[f(\boldsymbol{x}_{t},t)]$ denote the expected value for a generic, smooth enough function $f(\boldsymbol{x}_{t},t)$. The Kolmogorov backward equation corresponding to (8)-(9) takes the form [22]

[TABLE]

where

[TABLE]

and $\partial u/\partial\boldsymbol{x}$ is a row vector denoting the gradient of the scalar function $u$ with respect to the vector $\boldsymbol{x}$ (thus $\partial u/\partial\boldsymbol{x}\;\boldsymbol{a}(\boldsymbol{x})$ is the scalar product of the two vectors). We look for a solution of (10) in the form of a power series expansion $u=u_{0}+\varepsilon u_{1}+\varepsilon^{2}u_{2}+\ldots$

Introducing this ansatz into (10) and equating the coefficients of the same powers of $\varepsilon$ yields the hierarchy of equations

[TABLE]

The first equation in the hierarchy implies that $u_{0}$ does not depend on $y$, so that $u_{0}=u_{0}(\boldsymbol{x},t)$.

The other equations are of the type $\Lambda_{2}u_{n}=b_{n}$. According to Fredholm alternative theorem, these equations are solvable provided that a function $\psi$ exists such that: (1) $\Lambda_{2}^{*}\,\psi=0$, where $\Lambda_{2}^{*}$ is the conjugate operator of $\Lambda_{2}$, and (2) each $b_{n}$ satisfies $(b_{n},\psi)=0$, where $(\cdot,\cdot)$ denotes the inner product in the $L^{2}$ Banach space. Taking into account (13), $\Lambda_{2}^{*}\,\psi=0$ implies that $\psi$ is the stationary distribution of the Ornstein-Uhlenbeck process (9), thus [25]

[TABLE]

As a consequence, condition $(b_{n},\psi)=0$ amounts to require that each term $b_{n}$ averages to zero with respect to $y$

[TABLE]

Taking into account that $\langle y\rangle_{y}=0$, it is straightforward to verify that $\langle\Lambda_{1}u_{0}\rangle_{y}=0$. Thus equation (15) is solvable, and direct substitution shows that

[TABLE]

Similarly, equation (16) is solvable if

[TABLE]

The three averages on the left hand side can be expressed as

[TABLE]

where $\partial\boldsymbol{B}(\boldsymbol{x})/\partial\boldsymbol{x}$ is the Jacobian matrix of the vector function $\boldsymbol{B}$ with respect to the $\boldsymbol{x}$ variables with elements $(\partial\boldsymbol{B}(\boldsymbol{x})/\partial\boldsymbol{x})_{ij}=\partial B_{i}(\boldsymbol{x})/\partial x_{j}$, while $\partial^{2}u_{0}(\boldsymbol{x},t)/\partial\boldsymbol{x}^{2}$ is the Hessian matrix of the scalar function $u_{0}$ with respect to $\boldsymbol{x}$ with elements $(\partial^{2}u_{0}(\boldsymbol{x},t)/\partial\boldsymbol{x}^{2})_{ij}=\partial^{2}u_{0}(\boldsymbol{x},t)/(\partial x_{i}\partial x_{j})$. Furthermore, to solve the last integral we used the fact that, from (17), $\langle y^{2}\rangle_{y}=D^{2}/2$.

Substituting equations (21)-(23) into (20) yields the Kolmogorov backward equation

[TABLE]

The corresponding Stratonovich SDE is (we use the symbol $\circ$ to denote Stratonovich stochastic integral)

[TABLE]

while the equivalent Itô SDE is

[TABLE]

By equivalent we mean that the SDEs (25) and (26) are interpreted according to different rules but they have the same solution. The Stratonovich SDE (25) can be transformed into the Itô SDE (26) (and vice versa) by addition (respectively, subtraction) of the Wong-Zakai drift correction term [22]. The solution of (25) and (26) is also a weak solution for the original system with colored noise (5)-(6). Weak means that the two solutions for a specific realization of the the noise process $dW_{t}$ are different in details, but they converge in probability, i.e. they have the same probability density function and therefore the same statistical properties.

Whether to use the Stratonovich SDE (25) or the equivalent Itô SDE (26) is at this point a matter of personal taste. As a rule of thumb, Stratonovich interpretation may be better suited for algebraic manipulations, since traditional calculus rules apply. By contrast, Itô interpretation requires a whole new set of calculus rules, known as Itô calculus [22], but it is more suitable for numerical simulations and for the calculation of expected quantities.

IV Phase-amplitude equations for nonlinear oscillators with colored noise

Let us now consider the case where equations (5)-(6) describe a nonlinear oscillator subject to colored noise. In the absence of random fluctuations, equation (5) reduces to the autonomous ordinary differential equation (ODE)

[TABLE]

We assume that (27) admits of an asymptotically stable $T$-periodic solution $\boldsymbol{x}_{s}(t)$, represented by a limit cycle in its state space.

We shall derive an equivalent description of system (25) (or (26)) in terms of phase and amplitude deviation variables, analogous to the one derived in [26, 18, 27]. The phase function used in our description coincides locally, in the neighborhood of the limit cycle, with the asymptotic phase defined in [6, 13, 28]. As a second step, we shall derive a phase reduced model, that describes the oscillator dynamics in terms of the phase variable alone.

For our purpose it is more convenient to work with the Itô SDE (26). The reason to prefer the Itô over the Stratonovich interpretation is that Itô integrals are adapted processes, i.e. state variables and the noise increment are independent. By contrast, in the Stratonovich interpretation state variables and noise increments are correlated, a property known as “anticipating nature” or “look in the future property” of the Stratonovich integral. The far reaching consequence is that when one tries to describe the dynamics by using only a subset of state variables, an additional piece of information is lost, represented by the correlation between eliminated variables and noise increments [17].

To make the paper self-contained, we introduce some notation. We consider a set of time dependent vectors $\{\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)\}$, forming a basis for $\mathds{R}^{n}$, for all $t$. These basis vectors can be conveniently constructed as follows: the vector $\boldsymbol{u}_{1}(t)$ is chosen as the unit vector tangent to the limit cycle at any $t$

[TABLE]

The remaining $n-1$ vectors $\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)$ can be chosen as the Floquet vectors (apart from the limit cycle tangent $\boldsymbol{u}_{1}(t)$) of the linearized variational equation [26, 18, 27]

[TABLE]

where $\boldsymbol{J}_{\boldsymbol{a}}(t)=\partial\boldsymbol{a}/\partial\boldsymbol{x}$ is the Jacobian matrix of the vector function $\boldsymbol{a}(\boldsymbol{x}(t))$ evaluated on the limit cycle $\boldsymbol{x}_{s}(t)$. Thus the vectors $\{\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)\}$ are independent, although in general they are not orthogonal. We construct the matrix $\boldsymbol{U}(t)=[\boldsymbol{u}_{1}(t),\ldots,\boldsymbol{u}_{n}(t)]$, and we define the reciprocal vectors $\boldsymbol{v}_{1}^{T}(t),\ldots,\boldsymbol{v}_{n}^{T}(t)$ to be the rows of the inverse matrix $\boldsymbol{V}(t)=\boldsymbol{U}^{-1}(t)$. Thus $\{\boldsymbol{v}_{1}(t),\ldots,\boldsymbol{v}_{n}(t)\}$ also span $\mathds{R}^{n}$ and the bi-orthogonality condition $\boldsymbol{v}_{i}^{T}\boldsymbol{u}_{j}=\boldsymbol{u}_{i}^{T}\boldsymbol{v}_{j}=\delta_{ij}$ for all $t$, holds. Finally we introduce matrices $\boldsymbol{Y}(t)=[\boldsymbol{u}_{2}(t),\ldots,\boldsymbol{u}_{n}(t)]$, $\boldsymbol{Z}(t)=[\boldsymbol{v}_{2}(t),\ldots,\boldsymbol{v}_{n}(t)]$, and the magnitude (in the $L^{2}$ norm) of the vector field evaluated on the limit cycle, $r(t)=|\boldsymbol{a}(\boldsymbol{x}_{s}(t))|$.

Following [26, 18, 27], we decompose the solution of (26) into two components

[TABLE]

The first component $\boldsymbol{x}_{s}(\theta_{t})$ represents the projection of the stochastic process $\boldsymbol{x}_{t}$ onto the limit cycle, evaluated at an unknown time instant $\theta_{t}$. The second component $\boldsymbol{Y}(\theta_{t})\boldsymbol{R}_{t}$ represents the distance between the solution and the limit cycle, measured along the directions spanned by the vectors $\boldsymbol{v}_{2},\ldots,\boldsymbol{v}_{n}$ at the random time $\theta_{t}$. Because $\boldsymbol{x}_{t}$ is a stochastic process, both $\theta:\mathds{R}\mapsto\mathds{R}$ and $\boldsymbol{R}:\mathds{R}\mapsto\mathds{R}^{n-1}$ are stochastic processes as well. Itô equations defining the time evolution of these stochastic processes can be found following the procedure given in [18, Theorem 3.1] and [27, Theorem 1], obtaining

[TABLE]

where (the ′ sign denotes the derivative with respect to $\theta$)

[TABLE]

[TABLE]

and

[TABLE]

The SDEs (31)-(32) describe phase noise in nonlinear oscillators with colored noise as a drift-diffusion process. The responsibility of random fluctuations to phase diffusion does not come as a surprise, since, contrary to the amplitude, phase deviation does not have a self-limiting mechanism. Phase deviations are not damped, and may eventually grow unbounded as time passes. However, noise is also responsible for phase drift, that is, it produces a shift in the position of the peaks of the oscillator frequency spectrum. Because of the nonlinear response of the oscillator, random forces applied at a certain angle are amplified, while other are reduced. This results in a net, non null contribution to the expected angular frequency. Noise induced frequency shift is also observed in nonlinear oscillators subject to white Gaussian noise [17, 18, 27], but in presence of colored noise there is an additional shift contribution, represented by the term $b_{\theta}$, that can be ascribed to the finite correlation time of the noise source.

V Phase equation

The phase and amplitude deviation SDEs (31)-(32) have the same solution as the Stratonovich SDE (25) or the Itô SDE (26), that in turn are characterized by the same statistical properties of the solution of the SDEs (5)-(6). Since (31)-(32) are exact, they are not easier to solve than the white noise approximated SDE (26). However, they can be used to derive a phase reduced model [26, 29, 30] that in turn can form the basis to find useful, albeit approximate, results. The main advantage of a phase reduced model is that methods for Markovian systems, e.g. Fokker-Planck and Kolmogorov equations, can be (comparatively) easily applied, and the obtained equations can be more easily solved, being one dimensional (single variable).

To derive a simplified phase equation, we exploit a stochastic averaging technique. First we observe that if the Floquet basis is used as vectors $\boldsymbol{u}_{i},\,\boldsymbol{v}_{j}$, then

[TABLE]

where $\mathcal{O}(\boldsymbol{R}^{2})$ denotes terms quadratic in the amplitude deviation components, see [27, Theorem 3] for a proof. Then in the neighborhood of the limit cycle the “deterministic” frequency shift $a_{\theta}$ becomes negligibly small, and all points rotate with a uniform angular frequency. In many of the practical applications noise perturbations are small if compared to deterministic effects, that is, $\hat{a}_{\theta}(\theta,\boldsymbol{R})$, $b_{\theta}(\theta,\boldsymbol{R})$ and $B_{\theta}(\theta,\boldsymbol{R})$ can be considered as perturbation terms. They either include some explicit small parameter, or the condition $D\ll 1$ holds. As a consequence, the stationary distribution for the angle is expected to remain close to the uniform distribution $p_{\text{st}}(\theta)=1/(2\pi)$, that describes the phase diffusion process in a nonlinear oscillator with uniform angular frequency [18, 27].

Similar considerations can be made for the amplitude deviation SDE (32). It can be shown that (see [27, Theorem 3])

[TABLE]

where $\boldsymbol{D}=\textrm{diag}[\nu_{2},\ldots,\nu_{n}]$ is a diagonal matrix whose entries are the Floquet characteristic exponents, with the exception of the structural one $\nu_{1}=0$. The amplitude deviation dynamics is the balance of two competing forces: random fluctuations drive the system out of the limit cycle, while the asymptotic stability of the limit cycle implies that the system is continuously pushed toward the periodic orbit. Electronic systems are usually strongly stable, meaning that $\text{Re}\{\nu_{i}\}\ll 0$ for all $i=2,\ldots,n$, and as a consequence amplitude fluctuations remain confined to a small neighborhood of the limit cycle. Thus we can linearize the amplitude deviation SDE around the noiseless solution $\boldsymbol{R}=0$, and after averaging with respect to the phase stationary distribution $p_{\text{st}}(\theta)=1/(2\pi)$ we obtain

[TABLE]

where (as usual $\partial\hat{\boldsymbol{a}}_{\boldsymbol{R}}/\partial\boldsymbol{R}$, $\partial\boldsymbol{b}_{\boldsymbol{R}}/\partial\boldsymbol{R}$ and $\partial\boldsymbol{B}_{\boldsymbol{R}}/\partial\boldsymbol{R}$ are the Jacobian matrices with respect to $\boldsymbol{R}$)

[TABLE]

and

[TABLE]

In general, the solution of the linear SDE (45) is not a Gaussian process, but the vector of the expected values $\boldsymbol{\mu}(t)=\langle\boldsymbol{R}_{t}\rangle$ and the matrix of second moments $\boldsymbol{P}(t)=\langle\boldsymbol{R}_{t}\,\boldsymbol{R}_{t}^{T}\rangle$ can be found solving the linear ODE [31]

[TABLE]

Finally, the first and second moment are used to obtain a phase reduced equation. Expanding the terms of the phase SDE (31) in Taylor series around $\boldsymbol{R}=0$, and averaging with respect to the amplitude deviation, yields

[TABLE]

where the functions $a_{\theta},\,\hat{a}_{\theta},\,b_{\theta},\,B_{\theta}$ and their derivatives are evaluated at $(\theta,0)$, and $R_{i}$ denotes the $i$-th component of vector $\boldsymbol{R}$.

We compare the phase equation (52) with the analogous equations obtained for a nonlinear oscillator subject to white Gaussian noise [18, 27]. Apart from 1 that represents the oscillator’s normalized angular frequency, the terms in the first two rows describe a frequency shift, whereas the terms in the last row describe a diffusion. $\hat{a}_{\theta}$ and its derivatives resolve the correlation between the phase and noise increments. These terms were already discussed for a nonlinear oscillator subject to white Gaussian noise [18, 27]. By contrast, $b_{\theta}$ and its derivatives describe the different action that colored noise exerts on the phase with respect to white noise only, due to the non null noise correlation time.

It is worth noticing that in the weak noise limit, if higher order contributions of amplitude fluctuations and the correlation resolving term are neglected (implying $\mu_{i}=0,\,P_{ij}=0$ and $\hat{a}_{\theta}(\theta,\boldsymbol{R})=0$, respectively), then the simplified phase equation is obtained

[TABLE]

Eq. (53) is the equivalent of the phase equations derived in [4, 6] for the case of colored noise, where $b_{\theta}(\theta,0)$ is a zero order approximation of the frequency shift effect produced by the finite noise correlation time.

In order to compare our results against previous literature on phase noise in oscillators subject to colored noise sources, we consider here the approach developed in [9] where higher order contributions of amplitude fluctuations are neglected. For a strongly stable limit cycle, fluctuations are expected to keep the trajectory in a small neighborhood of the limit cycle, so that amplitude noise plays no influence on the phase dynamics. The noisy solution is then approximated as a time shifted version of the noiseless limit cycle $\boldsymbol{x}_{t}=\boldsymbol{x}_{S}(\theta_{t})$, and a phase equation is readily derived (see [9], eq. (3)111The division by $r(\theta)$, absent in [9], is a normalization required to guarantee that the oscillator’s free running frequency is equal to one.). For our system (5), (6) the phase model in [9] reads

[TABLE]

Notice that this phase equation is further approximated in [9] to derive the noise spectra.

Fully neglecting the amplitude fluctuations may be a reasonable approximation for strongly stable oscillators, however a more detailed analysis shows that in some cases the amplitude noise impacts on the cycle frequency inducing a non-negligible shift [11, 17, 18], especially for some autonomous systems exploited in computational biology and neuroscience [19].

VI Examples

VI-A Stuart-Landau oscillator with colored noise

As a first example we consider a Stuart-Landau oscillator with colored noise. The reason to choose such a simple system is twofold. First, most of the analysis can be made analytically, making the example useful to illustrate the theory and techniques described in the previous sections. Second, because many of the equations admit of an exact solution, the example permits to assess the accuracy of exploited approximations.

The state equations are the following

[TABLE]

where $\alpha$ and $\beta$ are real parameters that define the oscillator free running frequency.

Applying the methodology described in section III we obtain the following SDE with modulated white Gaussian noise

[TABLE]

We can now take advantage of the particularly simple structure of the Stuart-Landau system. Because the SDE for the amplitude is independent on the phase, the Fokker-Planck equation for the amplitude is single variable

[TABLE]

The stationary distribution can be found analytically

[TABLE]

where $\mathcal{N}$ is a constant determined through a normalization condition $\int_{0}^{+\infty}p_{\text{st}}(\rho)d\rho=1$. It is worth noticing that the same result cannot be obtained if the original problem is considered, because in the SDE (56) the amplitude equation and the Ornstein-Uhlenbeck process are coupled.

The theoretical prediction (59) is compared to the amplitude stationary distribution obtained through numerical integration of (56) in figure 1. Milstein numerical integration scheme has been used in the simulation, and the probability to find the amplitude in the interval $\rho+d\rho$ has been evaluated as the fraction of time spent in that interval divided by the total simulation length. As expected, the accuracy of the white noise approximation increases as the noise correlation time $\tau$ decreases.

In absence of noise the Stuart-Landau oscillator admits of an asymptotically stable limit cycle $\boldsymbol{x}_{s}(t)=[(\alpha-\beta)t,1]^{T}$. The Floquet vectors are $\boldsymbol{u}_{1}(t)=[1,0]^{T}$, $\boldsymbol{u}_{2}(t)=[\beta,1]^{T}$, while the co-vectors are $\boldsymbol{v}_{1}(t)=[1,-\beta]^{T}$, $\boldsymbol{v}_{2}(t)=[0,1]^{T}$. It is straightforward to derive the phase and amplitude deviation equations

[TABLE]

The linearized SDE for the amplitude is

[TABLE]

and the stationary distribution obtained solving the associated Fokker-Planck equation is

[TABLE]

where $\hat{\mathcal{N}}$ is the normalization constant.

Figure 2 shows the comparison between the amplitude deviation stationary distribution for the full system (59), and its counterpart for the linearized system (63), for different values of $D$. As expected, the distribution of the linearized system approximates well the full distribution for small values of $D$. The stationary distribution for the phase, obtained using numerical integration, is shown in figure 3 for two different values of $D$. It confirms that the uniform distribution is a good approximation even for fairly large values of $D$.

Solving equations (50), (51), we find the first two moments

[TABLE]

Finally, taking stochastic expectation on both sides of the first of (60) and neglecting $\mathcal{O}(R^{3})$ terms we obtain the expected angular frequency

[TABLE]

A comparison with [9] can be made making explicit (54)

[TABLE]

and taking the stochastic expectation on both sides of (67). Using (3) we find the expected angular frequency

[TABLE]

Therefore, according to the model in [9], noise has no influence at all on the asymptotic expected angular frequency.

Figure 4 shows the expected normalized angular frequency for the Stuart-Landau oscillator with colored noise (56), the equivalent system with white noise (57), and the theoretical predictions (69) and (66), as functions of $D$. Under the hypothesis that the system is ergodic, the normalized expected frequencies for systems (56) and (57) have been obtained through numerical integration, using the time average

[TABLE]

for $t_{2}\gg t_{1}$.

An analysis of the phase equation (60) provides further information. Averaging over the amplitude, substituting (64) and (65), and neglecting $\mathcal{O}(R^{3})$ terms, proves that the angle variable is well approximated by a Brownian motion with drift222This process is sometime referred to simply as Brownian motion, whereas the case $\mu=0$ is called Standard Brownian motion, in accordance with the corresponding PDFs. $\theta(t)=\mu t+\sigma W_{t}$, where

[TABLE]

For the sake of simplicity we assume a perfectly localized initial condition $\theta(0)=0$, and that $\theta\in(-\infty,+\infty)$, with boundary conditions $p(\pm\infty,t)=0$. Then the phase has a normal distribution $p(\theta,t)\sim N(\mu t,\sigma^{2}t)$, and the auto-correlation is

[TABLE]

Figure 5 shows the PDF $p(\theta-t,t)$ for the phase deviation (i.e. the difference between $\theta(t)$ and the phase in absence of noise) at three different time instants.

The power spectral density (PSD) for the orbital noise component, calculated using Welch’s method together with an average over 200 realizations of the corresponding time-domain process, is shown in figure 6. The blue line is the PSD for the orbital noise $x_{s}=\rho(t)\cos\phi(t)$ in presence of the colored noise source as defined in (56). The noise correlation time is $\tau=0.5$ and the noise intensity is $D=0.4$. The red line is the PSD for the orbital noise component $x_{s}=(1+\mu)\cos[(\alpha-\beta)\theta(t)+\beta\mu]$, where $\theta$ is the solution of the reduced phase equation (52), again for $D=0.4$. The two PSDs show excellent agreement around the carrier frequency, while the system with colored noise shows a significantly reduced power content at high frequency (not shown in this figure). The black solid line is the PSD for the system without noise, here shown to put in evidence the frequency shift induced by noise. PSDs were obtained considering $10^{4}$ oscillations, divided into $2^{28}$ points, corresponding to a sampling rate of 5369 samples per second for the discrete Fourier transform calculation. Finally, the black dotted line marks the theoretical prediction for the expected angular frequency $\left\langle{d\theta}/{dt}\right\rangle$ as given by (66).

VI-B van der Pol oscillator with colored noise

As a second example we consider the nonlinear oscillator (van der Pol) shown in figure 7. The nonlinear resistor $N_{R}$ is assumed to be noiseless, with characteristic $i_{G}=g(v)={v^{3}}/{3}-v$. The random source on the right models environment and internal noise.

Using Kirchhoff current law it is straightforward to derive the state equations

[TABLE]

where $i_{N}$ is the integral of the stochastic process $s_{n}$. With the change of variables

[TABLE]

and assuming that the random source $s_{n}$ is a colored noise modulated by the current through the capacitor

[TABLE]

we obtain the state equations

[TABLE]

where $\alpha=\sqrt{{L}/{C}}$.

Transformation to the equivalent system with white noise yields

[TABLE]

Figures 8 and 9 show the PDF $p(x_{1},x_{2},t)$, at the same time instant, for the van der Pol oscillator with colored noise (77) and noise correlation time $\tau=0.5$ and $\tau=0.1$, respectively. Figure 10 show the PDF for the equivalent system with white Gaussian noise (78). The PDF has been computed from numerical simulations, the same initial condition has been used in all cases333Parameter values are conveniently chosen to simplify calculations and to highlight higher order contributions of noise..

Because the limit cycle and the Floquet vectors cannot be found analytically for the van der Pol oscillator, we resort to semi analytical techniques and numerical methods for the analysis. In particular, we have developed a methodology based on the following steps:

• •

First, the limit cycle of the noiseless system is determined in a semi analytical form using the Harmonic Balance technique [32, 33].

• •

The semi analytical expression of the limit cycle is used to determine the Floquet vectors and co-vectors. Because the example under investigation is a second order system, Floquet vectors an co-vectors are computed using the formulas given in [34]. For higher order systems, they can be found exploiting efficient numerical techniques [33, 35, 36, 37].

• •

The limit cycle, and the Floquet vectors and co-vectors are used to determine the functions in equations (33)-(42).

• •

The functions $\boldsymbol{M}$, $\boldsymbol{m}$, $\boldsymbol{N}$ and $\boldsymbol{n}$ given by (46)-(48) are calculated, using the uniform distribution $p_{\text{st}}(\theta)=1/2\pi$ for averaging. Equations (50) and (51) are solved to find $\langle R_{i}\rangle$ and $\langle R_{i}R_{j}\rangle$ for all $i,j$.

• •

The reduced phase equation is written and can be analyzed to determine the expected normalized angular frequency, frequency shift, diffusion constant and so on.

Figure 11 shows the expected normalized angular frequency for the van der Pol system with colored noise (77), as a function of $D$. The expected frequency exhibits little dependence on the noise correlation time, in fact the blue curves (the solid curve corresponds to $\tau=0.5$, while the dashed one to $\tau=0.1$) are very close. The red line represents the expected normalized angular frequency for the equivalent system with white Gaussian noise (78). The curve has been determined through numerical integration of the full phase-amplitude deviation SDEs. Finally the black curves represents the theoretical predictions. The dotted line is the expected angular frequency determined through numerical integration of (54), (55) [9]. The solid line is the theoretical prediction given by our phase reduced model, obtained using the methods described above.

As a further confirmation, we have computed the power spectral density for the numerical solution of the SDEs (77) and (78). Figure 12 shows the power spectral density calculated combining Welch’s method with an averaging over 200 realizations to the numerical solution of the van der Pol system with colored and white noise, respectively. The discrete Fourier transform was calculated on a signal with time length $10^{4}T_{0}$ ($T_{0}$ being the period of the noiseless oscillator determined exploiting the harmonic balance technique) and a sampling rate of $4207$ samples per second. The blue line is the PSD in the absence of noise, while the red line is PSD for the system with noise ($D=0.4$). Noise clearly shifts the position of the power peak. The black dashed line identifies the expected frequency found using our theoretical model.

VII Conclusions

In this paper we have investigated the effect of colored noise on phase noise in nonlinear oscillators. We used a general method that transforms nonlinear systems subject to colored noise, modeled as an Ornstein-Uhlenbeck process, into equivalent systems subject to white Gaussian noise. The original system has unmodulated (additive) noise, and as such it is free of the Itô-Stratonovich dilemma. The transformation leads to two equivalent stochastic differential equations for the Itô and the Stratonovich interpretations. Since the two equations are equivalent, they have the same solution and therefore which one should be used is just a matter of personal preference. An alternative description in terms of stochastic differential equations for the phase and the amplitude deviation is derived for the transformed system. The phase variable used coincide locally, in the neighborhood of the unperturbed limit cycle, with the asymptotic phase defined using the concept of isochrons. Using stochastic averaging, a reduced phase model is derived, where the system dynamics is described only in terms of the phase variable.

The reduced phase model describes phase noise as a drift-diffusion process. The shift in the expected frequency is related to the variance of the colored noise and noise correlation time. A numerical procedure is presented for the solution of the phase equation, that provides more accurate results than other previously proposed models.