Colored noise in oscillators. Phase-amplitude analysis and a method to avoid the Ito-Stratonovich dilemma
11 Colored noise in oscillators.Phase-amplitude analysis and a method to avoid theIt ˆo-Stratonovich dilemma
Michele Bonnin, Fabio L. Traversa, Fabrizio Bonani,
Senior Member, IEEE
Abstract —We investigate the effect of time-correlated noise onthe phase fluctuations of nonlinear oscillators. The analysis isbased on a methodology that transforms a system subject tocolored noise, modeled as an OrnsteinUhlenbeck process, into anequivalent system subject to white Gaussian noise. A descriptionin terms of phase and amplitude deviation is given for thetransformed system. Using stochastic averaging technique, theequations are reduced to a phase model that can be analyzedto characterize phase noise. We find that phase noise is a drift-diffusion process, with a noise-induced frequency shift relatedto the variance and to the correlation time of colored noise.The proposed approach improves the accuracy of previous phasereduced models.
Index Terms —Oscillator noise, phase noise, colored noise,stochastic differential equations (SDEs), Fokker-Planck equation,stochastic averaging, phase models.
I. I
NTRODUCTION
Oscillators and phase locked loops (PLLs) are fundamentalcomponents of electronic and optical systems. For instance,in digital systems they are used to establish a reference timeto synchronize operations. In communication systems, theyare used for frequency coding and decoding, and for channelselection.Noise sources, both intrinsic and external, are a majornuisance plaguing oscillator and PLL performance. They canbe classified as white (frequency independent) fluctuations,such as thermal noise in electrical circuits with resistiveelements or shot noise in semiconductor devices, and time-correlated (colored) noise sources. Among the latter, partic-ularly relevant in oscillators based on bipolar transistors andMOSFET devices used as radio frequency sources, we findLorentzian low-frequency noise and flicker noise. Flicker or /f fluctuations can be in several cases traced back to thesuperposition of low-frequency Lorentzian noise sources, thattherefore are of paramount importance in assessing oscillatorrandom variations.The performance and reliability of oscillators depend cru-cially on noise sources, which deteriorate the oscillator re-sponse and are responsible for phase noise and time jitter.Phase noise and time jitter are strictly related concepts, defin-ing oscillator short term frequency instabilities. In particular,phase noise is a frequency domain measure of the oscillator M. Bonnin and F. Bonani are with the Dipartimento di Elettronica eTelecomunicazioni, Politecnico di Torino, Corso Duca degli Abruzzi 24,10129 Torino, Italy. F.L. Traversa is with MemComputing Inc., San Diego,CA, 92093-0319, USA. spectral purity, while time jitter describes the time domainaccuracy of the oscillator waveforms. Phase noise spreadsthe bandwidth of the fundamental frequency of oscillatorsand may produce interference with neighboring channels,thus degrading the whole system performances. Therefore,characterizing phase noise in oscillators is a major problemfor 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 atechnique is often too simplistic and fails to capture essentialfeatures such as spectral dispersion. Inclusion of linear timevariant effects (LTV) can yield more accurate results [2], [3].The most rigorous treatment of phase noise in nonlinearoscillator perhaps dates back to [4], where the author de-composes the oscillator response into phase and magnitudecomponents, and successfully derived a differential equationfor the phase deviation. In [5], the oscillator response isdecomposed into orthogonal components, and equations forpurely phase and amplitude deviations are derived. Unfortu-nately, as shown in [6], using an orthogonal decomposition toseparate phase and amplitude deviations leads to inaccurateresults. The method proposed in [6], based on using Floquetvectors to project the response of the noisy oscillator ontoa shifted version of the unperturbed response, exploits alinear periodically time varying approximation of the oscillatorbehavior leading, ultimately, to a nonlinear phase equation.The idea that Floquet vectors constitute the ideal basis todecouple 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 stochasticcharacterization of oscillator phase noise due to both whiteand colored noise sources [9]–[16].In [4], [6], [9], simplified scalar stochastic differentialequations for the phase variable were derived, neglecting allcontributions of noise and amplitude fluctuations, i.e smalldeviations from the limit cycle, beyond the first order. Thisassumption is well justified for most electronic oscillators,that are typically subject to small noise and exhibit stronglystable limit cycle, thus leading to very accurate phase models.LTV techniques yield reliable and precise predictions for thephase diffusion process, but they fail to capture the frequencyshift phenomenon [11], [17], [18]. Such a frequency shift isusually negligibly small for strongly stable oscillators, suchas those customarily used in electronic systems, but may playa relevant role in autonomous systems from other fields, e.g. a r X i v : . [ n li n . AO ] M a y system biology and neuroscience [19].This paper proposes a novel methodology to analyze phasenoise in nonlinear oscillators subject to time-correlated noise,modeled as an Ornstein-Uhlenebeck process. Making use ofthe method proposed in [20], [21] the system with colorednoise is first transformed into an equivalent system subjectto white Gaussian noise. The advantage is twofold: first thetransformation allows to use the whole machinery of stochas-tic differential equations on a reduced dimensional system.Second, and more important, the transformation avoids theItˆo-Stratonovich dilemma [22]. Phase and amplitude deviationequations are derived for the transformed system, and then re-duced to a phase model, that describes the oscillator dynamicsin terms of the phase variable only. This description is the idealtool for phase noise analysis, since it gives an approximate yetvery accurate description of the phase dynamics. We show thatphase noise in nonlinear oscillators is a drift-diffusion process,that is, noise sources not only induce a spread in the oscillatorspectrum, but also a shift in the oscillation mean frequency.The frequency shift is related to the variance of the colorednoise and to the noise correlation time. Some examples arepresented to assess the validity of the approach.II. M ODELING
Consider the nonlinear system subject to random noisesource ˙ x t = a ( x t ) + B ( x t ) η t (1)where x t : R (cid:55)→ R n denotes the state of the system, a : R n (cid:55)→ R n is a smooth vector field that defines the systeminternal dynamics, B : R n (cid:55)→ R n is a smooth vector valuedfunction and η t : R (cid:55)→ R is a scalar function describingrandom fluctuations, both internal and external.Random fluctuations are often modeled as zero mean,Gaussian distributed white noise. A zero mean Gaussianwhite noise η t = ˙ W t is characterized by (cid:104) ˙ W t (cid:105) = 0 and (cid:104) ˙ W t ˙ W s (cid:105) = δ ( t − s ) , where δ is the Dirac delta function,while t and s denote two time instants. White Gaussian noiseis a reasonable approximation in the case where the typicaltime scales of the underlying deterministic dynamics are muchlarger than the noise correlation time (quasi-white approxima-tion). Unfortunately, white noise is rarely an accurate modelto represent all of the noise sources that induce fluctuations inreal electronic devices and systems. As a consequence, a Diracdelta correlation in time is more justified by mathematicalconvenience than being physically plausible. Real processeshave typically finite correlation times, and often /f powerspectra, e.g. flicker noise .A more realistic description of noise in electronic systemsis given by an exponentially correlated process. Better knownas colored noise, it can be modeled as an Ornstein-Uhlenbeckprocess [23]. In this case the noise source is modeled accord-ing to τ ˙ η t = − η t + D ˙ W t (2)where τ is a parameter proportional to the finite noise cor-relation time, and D is the diffusion constant. The Ornstein- Uhlenbeck process with deterministic initial state η is char-acterized by the expectation value (cid:104) η t (cid:105) = η exp (cid:18) − tτ (cid:19) (3)and by the correlation function (cid:104) η t η s (cid:105) = D τ exp (cid:18) − | t − s | τ (cid:19) (4)Thus we consider equations (1) and (2), that we rewrite adopt-ing the standard notation of stochastic differential equations(SDEs) d x t = [ a ( x t ) + B ( x t ) η t ] dt (5) τ dη t = − η t dt + D dW t (6)where W t is a Wiener process, i.e. the integral of a whitenoise.The SDEs (5)-(6) describe a diffusion process with un-modulated (additive) noise. Therefore the equations can beinterpreted using any of the two main interpretation schemes,namely Itˆo or Stratonovich, obtaining the same solution. Inthe full space ( x , η ) the system is Markovian, that is, futurestates are completely determined by the current state and thestochastic process shows no memory. Methods for analysisof Markovian systems, based on the Fokker-Planck or Kol-mogorov equations, are well developed. However, the practicalsolution of the equations obtained using these approachesmay become unbearable because of the large number of statevariables involved.A possible solution strategy amounts to study the systemdynamics in a reduced dimension space. However, if weconsider the ( x ) space only the system is non Markoviandue to the presence of multiplicative noise, as the incrementsof the state variables depend on the past history of thenoise process. Problems arise even in the simpler case of aquasi-white approximation: for τ → , (6) shows that theexternal fluctuations reduce to a white noise η t dt = DdW t .Substituting this approximation into (5) yields d x t = a ( x t ) dt + D B ( x t ) dW t (7)The white noise is now modulated by the state dependentfunction B ( x ) (multiplicative noise), therefore the questionarises whether (7) should be interpreted according to Itˆo orStratonovich.Applying the procedure presented in [20], [21], in the nextsection we shall derive a reduced description in the ( x ) spaceof problem (5)-(6) where the SDE system is transformed intoan SDE (for the state vector x only) subject to white Gaussiannoise. The main results are the following • The reduced system holds for small, but not necessarilyvanishing correlation time τ . • The reduced system resolves the Itˆo-Stratonovichdilemma, that is, we shall derive equivalent SDEs for thetwo interpretations. By equivalent, we mean two differentSDEs, interpreted following different rules, having thesame solution. Because the solution is unique, it is just amatter of personal preference to choose one interpretationrather than the other.
III. W
HITE NOISE APPROXIMATION
Dividing both sides of (6) by τ , substituting τ = ε andintroducing η t = y t /ε , equations (5)-(6) become d x t = (cid:20) a ( x t ) + 1 ε B ( x t ) y t (cid:21) dt (8) dy t = − ε y t dt + Dε dW t (9)Usually, in electronic systems the correlation time τ is smallcompared to the characteristic time constants for the deter-ministic part of the dynamics. Therefore, we can assume thecorrelation time τ small enough such that ε (cid:28) . Under thisassumption, equations (8) and (9) show a time scale separation,since the Ornstein-Uhlenbeck process y t is one order, inthe parameter ε , faster than the state variables x t . Noticethat a straightforward application of stochastic averaging [24]would lead to inconsistent results. In fact, since asymptotically (cid:104) y t (cid:105) = 0 , the averaged equation would simply coincide withthe deterministic (noiseless) system.Now, let u ( x, t ) = E [ f ( x t , t )] denote the expected value fora generic, smooth enough function f ( x t , t ) . The Kolmogorovbackward equation corresponding to (8)-(9) takes the form[22] ∂u∂t = (cid:18) Λ + 1 ε Λ + 1 ε Λ (cid:19) u (10)where Λ u = ∂u∂ x a ( x ) (11) Λ u = y ∂u∂ x B ( x ) (12) Λ u = − y ∂u∂y + D ∂ u∂y (13)and ∂u/∂ x is a row vector denoting the gradient of the scalarfunction u with respect to the vector x (thus ∂u/∂ x a ( x ) isthe scalar product of the two vectors). We look for a solutionof (10) in the form of a power series expansion u = u + εu + ε u + . . . Introducing this ansatz into (10) and equating the coefficientsof the same powers of ε yields the hierarchy of equations ε − : Λ u = 0 (14) ε − : Λ u = − Λ u (15) ε : Λ u = ∂u ∂t − Λ u − Λ u (16)...The first equation in the hierarchy implies that u does notdepend on y , so that u = u ( x , t ) .The other equations are of the type Λ u n = b n . Accordingto Fredholm alternative theorem, these equations are solvableprovided that a function ψ exists such that: (1) Λ ∗ ψ = 0 ,where Λ ∗ is the conjugate operator of Λ , and (2) each b n satisfies ( b n , ψ ) = 0 , where ( · , · ) denotes the inner productin the L Banach space. Taking into account (13), Λ ∗ ψ = 0 implies that ψ is the stationary distribution of the Ornstein-Uhlenbeck process (9), thus [25] ψ ( y ) = p st ( y ) = (cid:114) πD exp (cid:18) − y D (cid:19) (17)As a consequence, condition ( b n , ψ ) = 0 amounts to requirethat each term b n averages to zero with respect to y ( b n , ψ ) = (cid:104) b n (cid:105) y = (cid:90) R b n p st ( y ) dy = 0 (18)Taking into account that (cid:104) y (cid:105) y = 0 , it is straightforward toverify that (cid:104) Λ u (cid:105) y = 0 . Thus equation (15) is solvable, anddirect substitution shows that u ( x , y, t ) = y ∂u ( x , t ) ∂ x B ( x ) = − Λ − Λ u (19)Similarly, equation (16) is solvable if (cid:28) ∂u ∂t (cid:29) y − (cid:104) Λ u (cid:105) y + (cid:10) Λ Λ − Λ u (cid:11) y = 0 (20)The three averages on the left hand side can be expressed as (cid:90) R ∂u ( x , t ) ∂t p st ( y ) dy = ∂u ( x , t ) ∂t (21) (cid:90) R Λ u p st ( y ) dy = ∂u ( x , t ) ∂ x a ( x ) (22) (cid:90) R Λ Λ − Λ u p st ( y ) dy = − D (cid:20) ∂u ( x , t ) ∂ x ∂ B ( x ) ∂ x B ( x )+ B T ( x ) ∂ u ( x , t ) ∂ x B ( x ) (cid:21) (23)where ∂ B ( x ) /∂ x is the Jacobian matrix of the vectorfunction B with respect to the x variables with elements ( ∂ B ( x ) /∂ x ) ij = ∂B i ( x ) /∂x j , while ∂ u ( x , t ) /∂ x is theHessian matrix of the scalar function u with respect to x with elements ( ∂ u ( x , t ) /∂ x ) ij = ∂ u ( x , t ) / ( ∂x i ∂x j ) .Furthermore, to solve the last integral we used the fact that,from (17), (cid:104) y (cid:105) y = D / .Substituting equations (21)-(23) into (20) yields the Kol-mogorov backward equation ∂u ( x , t ) ∂t = ∂u ( x , t ) ∂ x a ( x ) + D (cid:20) ∂u ( x , t ) ∂ x ∂ B ( x ) ∂ x B ( x )+ B T ( x ) ∂ u ( x , t ) ∂ x B ( x ) (cid:21) (24)The corresponding Stratonovich SDE is (we use the symbol ◦ to denote Stratonovich stochastic integral) d x t = a ( x t ) dt + D B ( x t ) ◦ dW t (25)while the equivalent Itˆo SDE is d x t = (cid:20) a ( x t ) + D ∂ B ( x t ) ∂ x B ( x t ) (cid:21) dt + D B ( x t ) dW t (26)By equivalent we mean that the SDEs (25) and (26) areinterpreted according to different rules but they have the samesolution. The Stratonovich SDE (25) can be transformed intothe Itˆo SDE (26) (and vice versa) by addition (respectively, subtraction) of the Wong-Zakai drift correction term [22]. Thesolution of (25) and (26) is also a weak solution for the originalsystem with colored noise (5)-(6). Weak means that the twosolutions 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 thereforethe same statistical properties.Whether to use the Stratonovich SDE (25) or the equivalentItˆo SDE (26) is at this point a matter of personal taste. Asa rule of thumb, Stratonovich interpretation may be bettersuited for algebraic manipulations, since traditional calculusrules apply. By contrast, Itˆo interpretation requires a wholenew set of calculus rules, known as Itˆo calculus [22], but it ismore suitable for numerical simulations and for the calculationof expected quantities.IV. P HASE - AMPLITUDE EQUATIONS FOR NONLINEAROSCILLATORS WITH COLORED NOISE
Let us now consider the case where equations (5)-(6)describe a nonlinear oscillator subject to colored noise. In theabsence of random fluctuations, equation (5) reduces to theautonomous ordinary differential equation (ODE) d x dt = a ( x ) (27)We assume that (27) admits of an asymptotically stable T -periodic solution x s ( t ) , represented by a limit cycle in itsstate 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 [18], [26], [27]. The phasefunction used in our description coincides locally, in theneighborhood of the limit cycle, with the asymptotic phasedefined in [6], [13], [28]. As a second step, we shall derivea phase reduced model, that describes the oscillator dynamicsin terms of the phase variable alone.For our purpose it is more convenient to work with the ItˆoSDE (26). The reason to prefer the Itˆo over the Stratonovich in-terpretation is that Itˆo integrals are adapted processes, i.e. statevariables and the noise increment are independent. By contrast,in the Stratonovich interpretation state variables and noiseincrements are correlated, a property known as “anticipatingnature” or “look in the future property” of the Stratonovichintegral. The far reaching consequence is that when one tries todescribe the dynamics by using only a subset of state variables,an additional piece of information is lost, represented by thecorrelation between eliminated variables and noise increments[17].To make the paper self-contained, we introduce somenotation. We consider a set of time dependent vectors { u ( t ) , . . . , u n ( t ) } , forming a basis for R n , for all t . Thesebasis vectors can be conveniently constructed as follows: thevector u ( t ) is chosen as the unit vector tangent to the limitcycle at any t u ( t ) = a ( x s ( t )) | a ( x s ( t )) | (28) The remaining n − vectors u ( t ) , . . . , u n ( t ) can be chosen asthe Floquet vectors (apart from the limit cycle tangent u ( t ) )of the linearized variational equation [18], [26], [27] d (cid:101) x ( t ) dt = J a ( t ) (cid:101) x ( t ) (29)where J a ( t ) = ∂ a /∂ x is the Jacobian matrix of the vectorfunction a ( x ( t )) evaluated on the limit cycle x s ( t ) . Thusthe vectors { u ( t ) , . . . , u n ( t ) } are independent, although ingeneral they are not orthogonal. We construct the matrix U ( t ) = [ u ( t ) , . . . , u n ( t )] , and we define the reciprocalvectors v T ( t ) , . . . , v Tn ( t ) to be the rows of the inverse matrix V ( t ) = U − ( t ) . Thus { v ( t ) , . . . , v n ( t ) } also span R n and the bi-orthogonality condition v Ti u j = u Ti v j = δ ij for all t , holds. Finally we introduce matrices Y ( t ) =[ u ( t ) , . . . , u n ( t )] , Z ( t ) = [ v ( t ) , . . . , v n ( t )] , and the mag-nitude (in the L norm) of the vector field evaluated on thelimit cycle, r ( t ) = | a ( x s ( t )) | .Following [18], [26], [27], we decompose the solution of(26) into two components x t = x s ( θ t ) + Y ( θ t ) R t (30)The first component x s ( θ t ) represents the projection of thestochastic process x t onto the limit cycle, evaluated at anunknown time instant θ t . The second component Y ( θ t ) R t represents the distance between the solution and the limitcycle, measured along the directions spanned by the vectors v , . . . , v n at the random time θ t . Because x t is a stochasticprocess, both θ : R (cid:55)→ R and R : R (cid:55)→ R n − arestochastic processes as well. Itˆo equations defining the timeevolution of these stochastic processes can be found followingthe procedure given in [18, Theorem 3.1] and [27, Theorem1], obtaining dθ t = [1 + a θ ( θ t , R t ) + ˆ a θ ( θ t , R t ) + b θ ( θ t , R t )] dt + B θ ( θ t , R t ) dW t (31) d R t = [ L ( θ t ) R t + a R ( θ t , R t )++ˆ a R ( θ t , R t ) + b R ( θ t , R t )] dt ++ B R ( θ t , R t ) dW t (32)where (the (cid:48) sign denotes the derivative with respect to θ ) a θ ( θ, R ) = κ v T (cid:2) a ( x s + Y R ) − a ( x s ) − Y (cid:48) R (cid:3) (33) ˆ a θ ( θ, R ) = − κ v T (cid:20) Y (cid:48) B R ( θ, R ) B θ ( θ, R )+ 12 B θ ( θ, R )( x (cid:48)(cid:48) s + Y (cid:48)(cid:48) R ) (cid:21) (34) b θ ( θ, R ) = D κ v T ∂ B ( x s + Y R ) ∂ x B ( x s + Y R ) (35) B θ ( θ, R ) = D κ v T B ( x s + Y R ) (36) L ( θ ) = − Z T Y (cid:48) (37) a R ( θ, R ) = Z T (cid:2) a ( x s + Y R ) − Y (cid:48) R a θ ( θ, R ) (cid:3) (38) ˆ a R ( θ, R ) = − Z T (cid:20) Y (cid:48) R ˆ a θ ( θ, R ) + Y (cid:48) B R ( θ, R ) B θ ( θ, R )+ 12 B θ ( θ, R )( x (cid:48)(cid:48) s + Y (cid:48)(cid:48) R ) (cid:21) (39) b R ( θ, R ) = − Z T Y (cid:48) R b θ ( θ, R )+ D Z T ∂ B ( x s + Y R ) ∂ x B ( x s + Y R ) (40) B R ( θ, R ) = − Z T Y (cid:48) R B θ ( θ, R ) + D Z T B ( x s + Y R ) (41)and κ = (cid:0) r + v T Y (cid:48) R (cid:1) − (42)The SDEs (31)-(32) describe phase noise in nonlinearoscillators with colored noise as a drift-diffusion process.The responsibility of random fluctuations to phase diffusiondoes 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 growunbounded as time passes. However, noise is also responsiblefor phase drift, that is, it produces a shift in the position ofthe peaks of the oscillator frequency spectrum. Because of thenonlinear response of the oscillator, random forces applied ata certain angle are amplified, while other are reduced. Thisresults in a net, non null contribution to the expected angularfrequency. Noise induced frequency shift is also observed innonlinear oscillators subject to white Gaussian noise [17], [18],[27], but in presence of colored noise there is an additionalshift contribution, represented by the term b θ , that can beascribed to the finite correlation time of the noise source.V. P HASE EQUATION
The phase and amplitude deviation SDEs (31)-(32) havethe same solution as the Stratonovich SDE (25) or the ItˆoSDE (26), that in turn are characterized by the same statisticalproperties of the solution of the SDEs (5)-(6). Since (31)-(32) are exact, they are not easier to solve than the whitenoise approximated SDE (26). However, they can be usedto derive a phase reduced model [26], [29], [30] that inturn can form the basis to find useful, albeit approximate,results. The main advantage of a phase reduced model isthat methods for Markovian systems, e.g. Fokker-Planck andKolmogorov equations, can be (comparatively) easily applied,and the obtained equations can be more easily solved, beingone dimensional (single variable).To derive a simplified phase equation, we exploit a stochas-tic averaging technique. First we observe that if the Floquetbasis is used as vectors u i , v j , then a θ ( θ, R ) = O ( R ) (43)where O ( R ) denotes terms quadratic in the amplitude de-viation components, see [27, Theorem 3] for a proof. Thenin the neighborhood of the limit cycle the “deterministic” frequency shift a θ becomes negligibly small, and all pointsrotate with a uniform angular frequency. In many of thepractical applications noise perturbations are small if com-pared to deterministic effects, that is, ˆ a θ ( θ, R ) , b θ ( θ, R ) and B θ ( θ, R ) can be considered as perturbation terms. They eitherinclude some explicit small parameter, or the condition D (cid:28) holds. As a consequence, the stationary distribution for theangle is expected to remain close to the uniform distribution p st ( θ ) = 1 / (2 π ) , that describes the phase diffusion processin a nonlinear oscillator with uniform angular frequency [18],[27].Similar considerations can be made for the amplitude devi-ation SDE (32). It can be shown that (see [27, Theorem 3]) L ( θ ) R + a R ( θ, R ) = DR + O ( R ) (44)where D = diag [ ν , . . . , ν n ] is a diagonal matrix whoseentries are the Floquet characteristic exponents, with the ex-ception of the structural one ν = 0 . The amplitude deviationdynamics is the balance of two competing forces: randomfluctuations drive the system out of the limit cycle, while theasymptotic stability of the limit cycle implies that the system iscontinuously pushed toward the periodic orbit. Electronic sys-tems are usually strongly stable, meaning that Re { ν i } (cid:28) forall i = 2 , . . . , n , and as a consequence amplitude fluctuationsremain confined to a small neighborhood of the limit cycle.Thus we can linearize the amplitude deviation SDE around thenoiseless solution R = 0 , and after averaging with respect tothe phase stationary distribution p st ( θ ) = 1 / (2 π ) we obtain d R t = ( M R t + m ) dt + n dW t (45)where (as usual ∂ ˆ a R /∂ R , ∂ b R /∂ R and ∂ B R /∂ R are theJacobian matrices with respect to R ) M = D + (cid:28) ∂ ˆ a R ∂ R (cid:29) θ + (cid:28) ∂ b R ∂ R (cid:29) θ (46) m = (cid:104) ˆ a R (cid:105) θ + (cid:104) b R (cid:105) θ (47) n = (cid:104) B R (cid:105) θ (48)and (cid:104) f ( θ ) (cid:105) θ = (cid:90) π f ( θ ) p st ( θ ) dθ = 12 π (cid:90) π f ( θ ) dθ (49)In general, the solution of the linear SDE (45) is not a Gaussianprocess, but the vector of the expected values µ ( t ) = (cid:104) R t (cid:105) and the matrix of second moments P ( t ) = (cid:104) R t R Tt (cid:105) can befound solving the linear ODE [31] d µ dt = M µ + m (50) d P dt = M P + P M T + mµ T + µm T + nn T (51)Finally, the first and second moment are used to obtain aphase reduced equation. Expanding the terms of the phase SDE (31) in Taylor series around R = 0 , and averaging withrespect to the amplitude deviation, yields dθ = (cid:20) a θ + b θ + (cid:88) i (cid:18) ∂ ˆ a θ ∂R i + ∂b θ ∂R i (cid:19) µ i + 12 (cid:88) i,j (cid:18) ∂ a θ ∂R i ∂R j + ∂ ˆ a θ ∂R i ∂R j + ∂ b θ ∂R i ∂R j (cid:19) P ij (cid:21) dt + (cid:18) B θ + (cid:88) i ∂B θ ∂R i µ i + 12 (cid:88) i,j ∂ B θ ∂R i ∂R j P ij (cid:19) dW t (52)where the functions a θ , ˆ a θ , b θ , B θ and their derivatives areevaluated at ( θ, , and R i denotes the i -th component ofvector R .We compare the phase equation (52) with the analogousequations obtained for a nonlinear oscillator subject to whiteGaussian noise [18], [27]. Apart from 1 that represents theoscillator’s normalized angular frequency, the terms in the firsttwo rows describe a frequency shift, whereas the terms in thelast row describe a diffusion. ˆ a θ and its derivatives resolvethe correlation between the phase and noise increments. Theseterms were already discussed for a nonlinear oscillator subjectto white Gaussian noise [18], [27]. By contrast, b θ and itsderivatives describe the different action that colored noiseexerts on the phase with respect to white noise only, due tothe non null noise correlation time.It is worth noticing that in the weak noise limit, if higher or-der contributions of amplitude fluctuations and the correlationresolving term are neglected (implying µ i = 0 , P ij = 0 and ˆ a θ ( θ, R ) = 0 , respectively), then the simplified phase equationis obtained dθ = [1 + b θ ( θ, dt + B θ ( θ, dW t (53)Eq. (53) is the equivalent of the phase equations derived in[4], [6] for the case of colored noise, where b θ ( θ, is a zeroorder approximation of the frequency shift effect produced bythe finite noise correlation time.In order to compare our results against previous literatureon phase noise in oscillators subject to colored noise sources,we consider here the approach developed in [9] where higherorder contributions of amplitude fluctuations are neglected.For a strongly stable limit cycle, fluctuations are expectedto keep the trajectory in a small neighborhood of the limitcycle, so that amplitude noise plays no influence on the phasedynamics. The noisy solution is then approximated as a timeshifted version of the noiseless limit cycle x t = x S ( θ t ) , anda phase equation is readily derived (see [9], eq. (3) ). For oursystem (5), (6) the phase model in [9] reads dθ = (cid:20) v T ( θ ) B ( x S ( θ )) r ( θ ) η t (cid:21) dt (54) τ dη t = − η t dt + D dW t (55)Notice that this phase equation is further approximated in [9]to derive the noise spectra. The division by r ( θ ) , absent in [9], is a normalization required to guaranteethat the oscillator’s free running frequency is equal to one. Fully neglecting the amplitude fluctuations may be a reason-able approximation for strongly stable oscillators, however amore detailed analysis shows that in some cases the amplitudenoise impacts on the cycle frequency inducing a non-negligibleshift [11], [17], [18], especially for some autonomous systemsexploited in computational biology and neuroscience [19].VI. E
XAMPLES
A. Stuart-Landau oscillator with colored noise
As a first example we consider a Stuart-Landau oscillatorwith colored noise. The reason to choose such a simple systemis twofold. First, most of the analysis can be made analyti-cally, making the example useful to illustrate the theory andtechniques described in the previous sections. Second, becausemany of the equations admit of an exact solution, the examplepermits to assess the accuracy of exploited approximations.The state equations are the following dφ = (cid:0) α − βρ + ρ η t (cid:1) dtdρ = ( ρ − ρ + ρ η t ) dtτ dη t = − η t dt + D dW t (56)where α and β are real parameters that define the oscillatorfree running frequency.Applying the methodology described in section III weobtain the following SDE with modulated white Gaussiannoise dφ = (cid:20) α + (cid:18) D − β (cid:19) ρ (cid:21) dt + D ρ dW t dρ = (cid:2) ρ + (cid:0) D − (cid:1) ρ (cid:3) dt + D ρ dW t (57)We can now take advantage of the particularly simple structureof the Stuart-Landau system. Because the SDE for the ampli-tude is independent on the phase, the Fokker-Planck equationfor the amplitude is single variable ∂p∂t = − ∂∂ρ (cid:8)(cid:2) ρ + (cid:0) D − (cid:1) ρ (cid:3) p (cid:9) + D ∂ ∂ρ (cid:0) ρ p (cid:1) (58)The stationary distribution can be found analytically p st ( ρ ) = N ρ − ( D ) exp (cid:18) − D ρ (cid:19) (59)where N is a constant determined through a normalizationcondition (cid:82) + ∞ p st ( ρ ) dρ = 1 . It is worth noticing that thesame result cannot be obtained if the original problem isconsidered, because in the SDE (56) the amplitude equationand the Ornstein-Uhlenbeck process are coupled.The theoretical prediction (59) is compared to the amplitudestationary distribution obtained through numerical integrationof (56) in figure 1. Milstein numerical integration scheme hasbeen used in the simulation, and the probability to find theamplitude in the interval ρ + dρ has been evaluated as thefraction of time spent in that interval divided by the totalsimulation length. As expected, the accuracy of the whitenoise approximation increases as the noise correlation time τ decreases.In absence of noise the Stuart-Landau oscillator admits ofan asymptotically stable limit cycle x s ( t ) = [( α − β ) t, T . Reference (56)Analytical (59)
Reference (56)Analytical (59)
Figure 1. Stationary distribution for the amplitude of a Stuart-Landauoscillator, for different values of the correlation time τ . Blue bars are obtainedfrom a numerical solution of the oscillator subject to colored noise (56). Thered line is the theoretical stationary distribution (59) obtained with the whitenoise approximation. On the left the results for τ = 0 . , on the right for τ = 0 . . Other parameters are α = 4 , β = 2 , D = 0 . . -0.5 0 0.5 1012345 Figure 2. Stationary distribution for the amplitude deviation of a Stuart-Landau oscillator, for different values of D . Solid and dashed blue lines arethe stationary distributions (59) and (63), respectively, for D = 0 . . Solidand dashed red lines are the stationary distributions (59) and (63), respectively,for D = 0 . . The Floquet vectors are u ( t ) = [1 , T , u ( t ) = [ β, T ,while the co-vectors are v ( t ) = [1 , − β ] T , v ( t ) = [0 , T . Itis straightforward to derive the phase and amplitude deviationequations dθ = (cid:26) α − β (cid:20) − βR + (cid:18) D − β (cid:19) (1 + R ) − β (cid:0) D − (cid:1) (1 + R ) (cid:21)(cid:27) dt + Dα − β (1 + R )[1 − β (1 + R )] dW t (60) dR = (cid:20) R + (cid:0) D − (cid:1) (1 + R ) (cid:21) dt + D (1 + R ) dW t (61)The linearized SDE for the amplitude is dR = (cid:2)(cid:0) − D (cid:1) R + D (cid:3) dt + D dW t (62)and the stationary distribution obtained solving the associatedFokker-Planck equation is ˆ p st ( R ) = ˆ N exp (cid:18) D − D R + 2 R (cid:19) (63)where ˆ N is the normalization constant. Figure 3. Stationary distribution for the phase of a Stuart-Landau oscillator,for two different values of D . Left: D = 0 . . Right D = 0 . . Parametersare α = 4 , β = 2 . Figure 2 shows the comparison between the amplitudedeviation stationary distribution for the full system (59), andits counterpart for the linearized system (63), for differentvalues of D . As expected, the distribution of the linearizedsystem approximates well the full distribution for small valuesof D . The stationary distribution for the phase, obtained usingnumerical integration, is shown in figure 3 for two differentvalues of D . It confirms that the uniform distribution is a goodapproximation even for fairly large values of D .Solving equations (50), (51), we find the first two moments (cid:104) R (cid:105) = D − D (64) (cid:104) R (cid:105) = 2 D (cid:104) R (cid:105) + D − D (65)Finally, taking stochastic expectation on both sides of the firstof (60) and neglecting O ( R ) terms we obtain the expectedangular frequency (cid:28) dθdt (cid:29) =1 + 1 α − β (cid:26) D − β ) + D (1 − β ) (cid:104) R (cid:105) + (cid:20) D − β ) + 2 β (cid:21) (cid:104) R (cid:105) (cid:27) (66)A comparison with [9] can be made making explicit (54) dθ = (cid:18) − βα − β η t (cid:19) dt (67) τ dη t = − η t dt + D dW t (68)and taking the stochastic expectation on both sides of (67).Using (3) we find the expected angular frequency (cid:28) dθdt (cid:29) = 1 + η e − tτ (69)Therefore, according to the model in [9], noise has no influ-ence at all on the asymptotic expected angular frequency.Figure 4 shows the expected normalized angular frequencyfor the Stuart-Landau oscillator with colored noise (56), theequivalent system with white noise (57), and the theoreticalpredictions (69) and (66), as functions of D . Under thehypothesis that the system is ergodic, the normalized expectedfrequencies for systems (56) and (57) have been obtainedthrough numerical integration, using the time average (cid:28) dθdt (cid:29) = 1 α − β φ ( t ) − φ ( t ) t − t (70) Figure 4. Expected normalized angular frequency for a Stuart-Landauoscillator versus D . Blue lines: Numerical result for the system with colorednoise (56). Solid: τ = 0 . . Dashed: τ = 0 . . Red line: numerical resultfor the equivalent system with white Gaussian noise (57). Black dotted line:Asymptotic theoretical prediction (69). Black solid line: Theoretical prediction(66). Parameters are α = 4 , β = 2 . for t (cid:29) t .An analysis of the phase equation (60) provides furtherinformation. Averaging over the amplitude, substituting (64)and (65), and neglecting O ( R ) terms, proves that the anglevariable is well approximated by a Brownian motion withdrift θ ( t ) = µt + σW t , where µ =1 + 1 α − β (cid:26) D − β ) + D (1 − β ) (cid:104) R (cid:105) (71) + (cid:20) D − β ) + 2 β (cid:21) (cid:104) R (cid:105) (cid:27) (72) σ = Dα − β (cid:2) − β + (1 − β ) (cid:104) R (cid:105) − β (cid:104) R (cid:105) (cid:3) (73)For the sake of simplicity we assume a perfectly localizedinitial condition θ (0) = 0 , and that θ ∈ ( −∞ , + ∞ ) , withboundary conditions p ( ±∞ , t ) = 0 . Then the phase hasa normal distribution p ( θ, t ) ∼ N ( µt, σ t ) , and the auto-correlation is R ( θ ( t ) , θ ( s )) = (cid:115) min( t, s )max( t, s ) (74)Figure 5 shows the PDF p ( θ − t, t ) for the phase deviation(i.e. the difference between θ ( t ) and the phase in absence ofnoise) at three different time instants.The power spectral density (PSD) for the orbital noisecomponent, calculated using Welch’s method together withan average over 200 realizations of the corresponding time-domain process, is shown in figure 6. The blue line is thePSD for the orbital noise x s = ρ ( t ) cos φ ( t ) in presence of thecolored noise source as defined in (56). The noise correlationtime is τ = 0 . and the noise intensity is D = 0 . . Thered line is the PSD for the orbital noise component x s =(1 + µ ) cos[( α − β ) θ ( t ) + βµ ] , where θ is the solution of thereduced phase equation (52), again for D = 0 . . The two PSDs This process is sometime referred to simply as Brownian motion, whereasthe case µ = 0 is called Standard Brownian motion, in accordance with thecorresponding PDFs. -10 -5 0 500.511.52 Figure 5. PDF p ( θ − t, t ) at three different time instants. Solid lines: PDFwith (cid:104) R (cid:105) and (cid:104) R (cid:105) given by (64) and (65), respectively. Dashed lines: PDFobtained neglecting amplitude fluctuations, that is imposing (cid:104) R (cid:105) = (cid:104) R (cid:105) = 0 .Other parameters are α = 4 , β = 2 , D = 0 . . PS D ( d B / H z ) Figure 6. Power spectral density for the orbital noise component x s of theStuart–Landau oscillator computed using the Welch’s method. Blue line: PSDfor the full system with colored noise (56), noise intensity D = 0 . and noisecorrelation time τ = 0 . . Black solid line: PSD for the full system withoutnoise. Red line: PSD obtained from the reduced phase equation. Black dashedline: expected angular frequency (cid:68) dθdt (cid:69) given by (66). Parameters are α = 4 and β = 2 . show excellent agreement around the carrier frequency, whilethe system with colored noise shows a significantly reducedpower 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 inducedby noise. PSDs were obtained considering oscillations,divided into points, corresponding to a sampling rate of5369 samples per second for the discrete Fourier transformcalculation. Finally, the black dotted line marks the theoreticalprediction for the expected angular frequency (cid:104) dθ/dt (cid:105) as givenby (66). 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 / − v . The random source on the right models environmentand internal noise. Using Kirchhoff current law it is straight-forward to derive the state equations dv = u dt (cid:48) (75) du = (cid:20) − LC v − C g (cid:48) ( v ) u − C s n ( t (cid:48) ) (cid:21) dt (cid:48) (76) i G i L i C i N N R L C + − v Figure 7. Second order nonlinear oscillator. where i N is the integral of the stochastic process s n . With thechange of variables t = 1 √ LC t (cid:48) x = v x = √ LC u and assuming that the random source s n is a colored noisemodulated by the current through the capacitor s n ( t ) = (cid:114) CL u η t we obtain the state equations dx = x dtdx = (cid:2) − x + α (cid:0) − x (cid:1) x + x η t (cid:3) dtτ dη t = − η t dt + D dW t (77)where α = (cid:112) L/C .Transformation to the equivalent system with white noiseyields dx = x dtdx = (cid:20) − x + α (cid:0) − x (cid:1) x + D x (cid:21) dt + D x dW t (78)Figures 8 and 9 show the PDF p ( x , x , t ) , at the sametime instant, for the van der Pol oscillator with colored noise(77) and noise correlation time τ = 0 . and τ = 0 . ,respectively. Figure 10 show the PDF for the equivalent systemwith white Gaussian noise (78). The PDF has been computedfrom numerical simulations, the same initial condition hasbeen used in all cases .Because the limit cycle and the Floquet vectors cannot befound analytically for the van der Pol oscillator, we resortto semi analytical techniques and numerical methods for theanalysis. In particular, we have developed a methodologybased on the following steps: • First, the limit cycle of the noiseless system is determinedin a semi analytical form using the Harmonic Balancetechnique [32], [33]. • The semi analytical expression of the limit cycle is usedto determine the Floquet vectors and co-vectors. Becausethe example under investigation is a second order system,Floquet vectors an co-vectors are computed using theformulas given in [34]. For higher order systems, theycan be found exploiting efficient numerical techniques[33], [35]–[37]. • The limit cycle, and the Floquet vectors and co-vectorsare used to determine the functions in equations (33)-(42). Parameter values are conveniently chosen to simplify calculations and tohighlight higher order contributions of noise. Figure 8. Probability density function for the van der Pol oscillator withcolored noise. Parameters are α = 0 . , D = 0 . . Noise correlation time is τ = 0 . .Figure 9. Probability density function for the van der Pol oscillator withcolored noise. Parameters are α = 0 . , D = 0 . . Noise correlation time is τ = 0 . . • The functions M , m , N and n given by (46)-(48) arecalculated, using the uniform distribution p st ( θ ) = 1 / π for averaging. Equations (50) and (51) are solved to find (cid:104) R i (cid:105) and (cid:104) R i R j (cid:105) for all i, j . • The reduced phase equation is written and can be an-alyzed to determine the expected normalized angular
Figure 10. Probability density function for the van der Pol oscillator withwhite Gaussian noise. Parameters are α = 0 . , D = 0 . . Figure 11. Expected normalized frequency for a van der Pol oscillatorversus D . Blue lines: Numerical result for the system with colored noise(77). Solid: τ = 0 . . Dotted: τ = 0 . . Red line: numerical result for theequivalent system with white Gaussian noise (78). Black dotted line: expectedfrequency obtained through numerical integration of (54), (55). Black solidline: Theoretical prediction using the proposed method. Parameter α = 0 . . frequency, frequency shift, diffusion constant and so on.Figure 11 shows the expected normalized angular frequencyfor the van der Pol system with colored noise (77), as afunction of D . The expected frequency exhibits little depen-dence on the noise correlation time, in fact the blue curves(the solid curve corresponds to τ = 0 . , while the dashedone to τ = 0 . ) are very close. The red line representsthe expected normalized angular frequency for the equivalentsystem with white Gaussian noise (78). The curve has beendetermined through numerical integration of the full phase-amplitude deviation SDEs. Finally the black curves representsthe theoretical predictions. The dotted line is the expectedangular frequency determined through numerical integrationof (54), (55) [9]. The solid line is the theoretical predictiongiven by our phase reduced model, obtained using the methodsdescribed above.As a further confirmation, we have computed the powerspectral density for the numerical solution of the SDEs (77)and (78). Figure 12 shows the power spectral density cal-culated combining Welch’s method with an averaging over200 realizations to the numerical solution of the van derPol system with colored and white noise, respectively. Thediscrete Fourier transform was calculated on a signal with timelength T ( T being the period of the noiseless oscillatordetermined exploiting the harmonic balance technique) and asampling rate of samples per second. The blue line isthe PSD in the absence of noise, while the red line is PSDfor the system with noise ( D = 0 . ). Noise clearly shifts theposition of the power peak. The black dashed line identifiesthe expected frequency found using our theoretical model.VII. C ONCLUSIONS
In this paper we have investigated the effect of colorednoise on phase noise in nonlinear oscillators. We used ageneral method that transforms nonlinear systems subject tocolored noise, modeled as an Ornstein-Uhlenbeck process,into equivalent systems subject to white Gaussian noise. Theoriginal system has unmodulated (additive) noise, and as suchit is free of the Itˆo-Stratonovich dilemma. The transformation PS D ( d B / H z ) PS D ( d B / H z ) Figure 12. Power Spectral Density for numerical solution of a van der Poloscillator computed using Welch’s algorithm. Left: PSD for the system withcolored noise (77), τ = 0 . . Right: PSD for the equivalent system with whiteGaussian noise (78). The black dashed line identifies the expected frequencyobtained using our theoretical model. Parameter α = 0 . . leads to two equivalent stochastic differential equations forthe Itˆo and the Stratonovich interpretations. Since the twoequations are equivalent, they have the same solution andtherefore which one should be used is just a matter of personalpreference. An alternative description in terms of stochasticdifferential equations for the phase and the amplitude deviationis derived for the transformed system. The phase variable usedcoincide locally, in the neighborhood of the unperturbed limitcycle, with the asymptotic phase defined using the concept ofisochrons. Using stochastic averaging, a reduced phase modelis derived, where the system dynamics is described only interms of the phase variable.The reduced phase model describes phase noise as a drift-diffusion process. The shift in the expected frequency is relatedto the variance of the colored noise and noise correlation time.A numerical procedure is presented for the solution of thephase equation, that provides more accurate results than otherpreviously proposed models.R EFERENCES[1] D. B. Leeson, “A simple model of feedback oscillator noise spectrum,”
Proceedings of the IEEE , vol. 54, no. 2, pp. 329–330, 1966.[2] E. Hafner, “The effects of noise in oscillators,”
Proceedings of the IEEE ,vol. 54, no. 2, pp. 179–198, 1966.[3] K. Kurokawa, “Noise in synchronized oscillators,”
IEEE Transactionson microwave theory and techniques , vol. 16, no. 4, pp. 234–240, 1968.[4] F. X. Kaertner, “Analysis of white and f − α noise in oscillators,” International Journal of Circuit Theory and Applications , vol. 18, no. 5,pp. 485–519, 1990.[5] A. Hajimiri and T. H. Lee, “A general theory of phase noise in electricaloscillators,”
IEEE Journal of Solid-State Circuits , vol. 33, no. 2, pp.179–194, 1998.[6] A. Demir, A. Mehrotra, and J. Roychowdhury, “Phase noise in oscil-lators: A unifying theory and numerical methods for characterization,”
IEEE Transactions on Circuits and Systems I: Fundamental Theory andApplications , vol. 47, no. 5, pp. 655–674, May 2000.[7] G. J. Coram, “A simple 2-d oscillator to determine the correct de-composition of perturbations into amplitude and phase noise,”
IEEETransactions on Circuits and Systems I: Fundamental Theory andApplications , vol. 48, no. 7, pp. 896–898, 2001.[8] F. L. Traversa and F. Bonani, “Oscillator noise: a nonlinear perturbativetheory including orbital fluctuations and phase-orbital correlation,”
IEEETransactions on Circuits and Systems I: Regular Papers , vol. 58, no. 10,pp. 2485–2497, October 2011.[9] A. Demir, “Phase noise and timing jitter in oscillators with colored-noisesources,”
IEEE Transactions on Circuits and Systems I: FundamentalTheory and Applications , vol. 49, no. 12, pp. 1782–1791, 2002.[10] J. P. Gleeson, “Phase diffusion due to low-frequency colored noise,”
IEEE Transactions on Circuits and Systems II: Express Briefs , vol. 53,no. 3, pp. 183–186, 2006. [11] F. O’Doherty and J. P. Gleeson, “Phase diffusion coefficient for oscil-lators perturbed by colored noise,” IEEE Transactions on Circuits andSystems II: Express Briefs , vol. 54, no. 5, pp. 435–439, 2007.[12] P. Maffezzoni, “Frequency-shift induced by colored noise in nonlinearoscillators,”
IEEE Transactions on Circuits and Systems II: ExpressBriefs , vol. 54, no. 10, pp. 887–891, 2007.[13] T. Djurhuus, V. Krozer, J. Vidkjær, and T. K. Johansen, “Oscillatorphase noise: a geometrical approach,”
IEEE Transactions on Circuitsand Systems Part I: Regular Papers , vol. 56, no. 7, pp. 1373–1382,2009.[14] T. Nagashima, X. Wei, H.-A. Tanaka, and H. Sekiya, “Locking rangederivations for injection-locked class-e oscillator applying phase reduc-tion theory,”
IEEE Transactions on Circuits and Systems I: RegularPapers , vol. 61, no. 10, pp. 2904–2911, 2014.[15] P. Maffezzoni, Z. Zhang, and L. Daniel, “A study of deterministic jitterin crystal oscillators,”
IEEE Transactions on Circuits and Systems I:Regular Papers , vol. 61, no. 4, pp. 1044–1054, 2014.[16] P. Maffezzoni, B. Bahr, Z. Zhang, and L. Daniel, “Reducing phase noisein multi-phase oscillators,”
IEEE Transactions on Circuits and SystemsI: Regular Papers , vol. 63, no. 3, pp. 379–388, 2016.[17] M. Bonnin and F. Corinto, “Phase noise and noise induced frequencyshift in stochastic nonlinear oscillators,”
IEEE Transactions on Circuitsand Systems I: Regular Papers , vol. 60, no. 8, pp. 2104–2115, 2013.[18] M. Bonnin, F. L. Traversa, and F. Bonani, “Influence of amplitudefluctuations on the noise-induced frequency shift of noisy oscillators,”
IEEE Transactions on Circuits and Systems II: Express Briefs , vol. 63,no. 7, pp. 698–702, 2016.[19] L. Tait, K. Wedgwood, K. Tsaneva-Atanasova, J. T. Brown, andM. Goodfellow, “Control of clustered action potential firing in a mathe-matical model of entorhinal cortex stellate cells,”
Journal of TheoreticalBiology , vol. 449, pp. 23–34, 2018.[20] D. Givon, R. Kupferman, and A. Stuart, “Extracting macroscopicdynamics: model problems and algorithms,”
Nonlinearity , vol. 17, no. 6,p. R55, 2004.[21] G. Pavliotis and A. Stuart,
Multiscale methods: averaging and homog-enization . Springer Science & Business Media, 2008.[22] B. Øksendal,
Stochastic Differential Equations , 6th ed. Berlin:Springer–Verlag, 2003.[23] D. T. Gillespie, “Exact numerical simulation of the Ornstein-Uhlenbeckprocess and its integral,”
Physical Review E , vol. 54, no. 2, p. 2084,1996.[24] R. Z. Khasminskii, “On the principle of averaging the Ito’s stochasticdifferential equations,”
Kibernetika , vol. 4, no. 3, pp. 260–279, 1968.[25] C. W. Gardiner,
Handbook of Stochastic Methods . Berlin: Springer,1985.[26] F. L. Traversa, M. Bonnin, F. Corinto, and F. Bonani, “Noise inoscillators: a review of state space decomposition approaches,”
Journalof Computional Electronics , vol. 14, no. 1, pp. 51–61, 2015.[27] M. Bonnin, “Amplitude and phase dynamics of noisy oscillators,”
International Journal of Circuit Theory and Applications , vol. 45, pp.636–659, 2017.[28] A. Mauroy and I. Mezic, “Global computation of phase-amplitudereduction for limit-cycle dynamics,” arXiv preprint arXiv:1803.07379 ,2018.[29] H. Nakao, “Phase reduction approach to synchronisation of nonlinearoscillators,”
Contemporary Physics , vol. 57, no. 2, pp. 188–214, 2016.[30] L. Freitas, L. A. Torres, and L. A. Aguirre, “Phase definition to assesssynchronization quality of nonlinear oscillators,”
Physical Review E ,vol. 97, no. 5, p. 052202, 2018.[31] P. E. Kloeden and E. Platen,
Numerical solution of stochastic differentialequations springer-verlag . New York: Springer, 1992.[32] M. Bonnin, F. Corinto, and M. Gilli, “Periodic oscillations in weaklyconnected cellular nonlinear networks,”
IEEE Transactions on Circuitsand Systems I: Regular Papers , vol. 55, no. 6, pp. 1671–1684, 2008.[33] F. L. Traversa, F. Bonani, and S. Donati Guerrieri, “A frequency-domainapproach to the analysis of stability and bifurcations in nonlinear systemsdescribed by differential-algebraic equations,”
International Journal ofCircuit Theory and Applications , vol. 36, no. 4, pp. 421–439, 2008.[34] M. Bonnin, F. Corinto, and M. Gilli, “Phase space decomposition forphase noise and synchronization analysis of planar nonlinear oscillators,”
IEEE Transactions on Circuits and Systems II: Express Briefs , vol. 59,no. 10, pp. 638–642, 2012.[35] F. L. Traversa and F. Bonani, “Frequency domain evaluation of theadjoint Floquet eigenvectors for oscillator noise characterization,”
IETCirc Device Syst , vol. 5, no. 1, pp. 46–51, 2011. [36] ——, “Improved harmonic balance implementation of Floquet analysisfor nonlinear circuit simulation,”
AEU - International Journal of Elec-tronics and Communications , vol. 66, no. 5, pp. 357–363, may 2012.[37] ——, “Selective determination of Floquet quantities for the efficientassessment of limit cycle stability and of oscillator noise,”