Real-time estimation of phase and amplitude with application to neural data

Computation of the instantaneous phase and amplitude via the Hilbert Transform is a powerful tool of data analysis. This approach finds many applications in various science and engineering branches but is not proper for causal estimation because it requires knowledge of the signal’s past and future. However, several problems require real-time estimation of phase and amplitude; an illustrative example is phase-locked or amplitude-dependent stimulation in neuroscience. In this paper, we discuss and compare three causal algorithms that do not rely on the Hilbert Transform but exploit well-known physical phenomena, the synchronization and the resonance. After testing the algorithms on a synthetic data set, we illustrate their performance computing phase and amplitude for the accelerometer tremor measurements and a Parkinsonian patient’s beta-band brain activity.


Introduction
Estimation of the instantaneous phase and amplitude is essential for electrical and mechanical engineering, synchronization studies of oscillatory systems of different nature, time series analysis of physiological data, and, in particular, for neuroscience [1][2][3][4] . A special application is developing efficient sensing algorithms for adaptive deep brain stimulation, a recent advancement of a widely used treatment option for Parkinson's disease and other neurological disorders 5 . One of the directions in this development is to adjust stimulation parameters according to a peripheral or neurophysiological signal's phase and amplitude computed on the fly [6][7][8][9][10][11][12] . The popular approach to phase and amplitude estimation is to exploit the analytic signal approach based on the Hilbert Transform (HT) or, equivalently, the wavelet transform with a complex wavelet 1,[13][14][15] . However, this widely-used tool is non-causal and therefore not appropriate for real-time analysis, whereas causal estimation of phase and amplitude is often crucial for closed-loop control of complex systems. Despite several attempts to adapt the HT for a causal measurement 16,17 , the reliable and fast estimation of phase and amplitude of real-world signals remains challenging.
In this study, we do not rely on the Hilbert Transform. Instead, we follow another development line and extend our previous approach 18,19 for real-time estimation based on the oscillation theory and nonlinear dynamics. The main idea is as follows. Suppose we have an oscillator, e.g., an electronic circuit, which amplitude and phase we can monitor. Next, suppose we feed our measurement to this oscillator. We choose the oscillator so that there is a one-to-one correspondence between the oscillator's phase and amplitude (what we can monitor) and the measured signal (what we want to determine). If this correspondence is achieved, we recompute the monitored state of the oscillator into desired quantities. Thus, the oscillator acts as a measuring device. Indeed, we will not exploit a physical oscillator but use a simple computer program that simulates the oscillator's dynamics. In our approach, we rely on two well-known physical effects: linear resonance and synchronization. We present and compare three techniques that provide a fast estimation of phases and amplitudes, using only the past of the time series. We test the algorithms on model data and apply them to neural time series. Namely, demonstrating our approach's efficiency, we causally estimate the phase and amplitude of the accelerometer tremor measurements and Parkinsonian patients' beta-band brain activity and compare the results with non-causal HT-based analysis.

Hilbert Transform vs causal estimation
A common practice in extracting the amplitude and the phase of an input signal s(t) is based on application of the Hilbert transform. The HT analysis provides proper estimates for the instantaneous phase of s(t) for narrow-band one-component signals with slowly varying amplitude and frequency. We emphasize that, though formally one can compute the HT for an arbitrary signal s(t), not in all cases the extraction of the amplitude and the phase will lead to reasonable results. We also stress that HT is a non-causal operation: to extract features at time instant t , one has to know both the signals' past t < t and future t > t . For a further detailed discussion of the HT's properties and practical implementation, see Methods.
The main goal of this paper is not to extend HT approach to complex signals (see 20,21 for examples of such an extension), but to provide and explore causal alternatives to the HT approach. In what follows, we use the HT-based amplitude a H and the HT-based phase ϕ H as a "gold standard" for testing our algorithms. Strictly speaking, this is reasonable for narrow-band signals only. Definition and determination of the proper phase and amplitude for a complex signal represent a challenging problem that we do not address here. Instead, we use an operational approach: a proper amplitude should correctly represent an envelope of the signal; a proper phase should gain 2π at each oscillatory cycle. Below, we present three techniques for causal computations and compare the computed amplitudes and phases with non-causally obtained a H , ϕ H .

Phase locking approach
The first technique exploits the ideas from the synchronization theory. It is well-known that a force s(t) acting on a limit-cycle oscillator can entrain (lock) it. It means that the oscillator's frequency becomes equal to that of the force, and their phases differ by a constant. Thus, the phase of the locked limit-cycle oscillator will correspond to the phase of the signal. For our purposes, it is helpful to use the simplest oscillator model, the so-called phase oscillator. To ensure the phase-locking to the force, we have to adjust the oscillator's frequency to the signal's frequency. We assume that we do not know the latter a priori, but can only roughly estimate it. We propose a simple approach that starts with this estimate and automatically tunes the "device's" frequency to ensure the locking and thus provides the instantaneous phase ϕ L (t). The amplitude is not determined with this approach. One can treat the suggested scheme as a software implementation of a phase-locked loop 22 . Technically, the algorithm reduces to solving differential equation incorporating measured data given at discrete time points; for details of the technique and its implementation, see Methods.

Nonresonant linear filter
The second technique relies on the resonance effect. Our measuring "device" consists now of two linear damped oscillators. The oscillators' frequency is much larger than the frequency of the signal, i.e., the system is far from resonance. We choose the damping parameters to ensure that (i) phase of the first linear oscillator equals that of the input and that (ii) amplitude of the second one and the input relate by a known constant multiplicator. The technique yields both phase ϕ N (t) and amplitude a N (t), where the index N stands for "non-resonant".

Resonant linear filter
Our third approach adopts the technique used for model studies in our previous publications 18,19 . The corresponding "device" consists of a linear oscillator in resonance with the measured signal and of an integrating unit. It also provides both the phase and the amplitude that we denote as ϕ R and a R , respectively. The method exploits the known relation between the resonant oscillator's phase and amplitude and those of the input. Additionally, the resonant oscillator acts as a bandpass filter for experimental data.
Technically, both latest techniques require the numerical solution of linear differential equations with the input signal s(t). We present a detailed description of the techniques and the developed numerical schemes in the Methods section. Like the phase-locking algorithm, both techniques include an automated frequency-tuning algorithm to adjust the systems to the a priori unknown signal's frequency.

Artificial data
First, we test our approach on artificial data. All algorithms work well with simple narrow-band signals like slowly modulated sine waves. So, we do not show these results and proceed with a more complicated case. The signal is where ψ(t) = t + 5 sin(Ω 2 t) , its waveform has three harmonics. The signal is amplitude-and phase-modulated, see Fig. 1a; it is sampled with ∆ = 0.01 to yield s(∆k) = s k . The frequencies are Ω 1 = √ 2/30, Ω 2 = √ 5/60. Another test signal iss k = s k + ξ k , where ξ k are Gaussian random numbers with zero mean and standard deviation 0.05, see Fig. 2a. For testing all three algorithms, we set the initial frequency of the device to be 10% higher than the actual value; in this way, we imitate the imprecision on initial frequency estimation. The other parameters are given in the Methods.   Fig. 1a). When the signal's amplitude is very small, the noise dominates and phase determination becomes complicated. The HT approach fails here (blue line in (b,c)), while both locked and resonant oscillator "devices" provide reasonable results. The performance of the non-resonant technique is poor: when the amplitude nearly vanishes, the phase estimated by this technique (not shown) is not better than the Hilbert phase.
performance is comparable or even better than that of two other techniques. In the presence of noise the performance is poor because noisy perturbations excite high-frequency oscillation with the frequency ω. The amplitude estimation is slightly worse than that via HT, see Fig. 2a. Thus, for processing without a filter, this technique is not optimal.

Resonant linear oscillator (filter).
This technique also demonstrates efficient phase estimation, see Fig. 1d for the noise-free data and Fig. 2c for the noisy case. Actually, the results for the latter are hardly distinguishable from those in Fig. 2b, though the frequency adaptation of the resonant oscillator is faster. The technique also provides the instantaneous amplitude, see Fig. 2a.
Notice that the Hilbert amplitude a H for the multi-component signals is not the perfect envelope. The real-time amplitude a R is not perfect either, but is much more smooth than a H .

Estimation without an additional filter: tremor data
Tremor data was acquired from a patient with essential tremor using a 3D accelerometer (TMSi, Oldenzaal, The Netherlands) attached to the right index finger. The signal was sampled at 2048 Hz sampling rate using a Porti amplifier (TMSi). The patient stretched their arms out in front of their chest holding a bottle in order to provoke postural tremor. Only data from one axis was analyzed. The tremor time series (Fig. 3a) is relatively simple: it resembles the test signal from the previous section. The tremor data is amplitude modulated, and its power spectrum also contains several harmonics (Fig. 3b). However, the real-world signals are naturally more complicated than the simple test data series Eq. (1). Analyzing them, one frequently faces additional difficulties. For the tremor signal under consideration these features are: a drift of the baseline; presence of epochs when the amplitude vanishes; and outliers. To cope with the baseline fluctuation, we exploit the following causal detrending algorithm. For each new measurement point, we remove the mean value computed over several previous cycles. For the sake of computational speed, we update this value several times per period. We denote the detrended signal bys k (Fig. 3c). Panels (e,f) demonstrate the Hilbert phase ϕ H (blue) and causal phases ϕ L (red) and ϕ R (magenta). For epochs with vanishing amplitude, causal phases do not exhibit noisy jumps as the Hilbert phase. The causal phases are stable with respect to outliers (f). Figure 3c compares the Hilbert amplitude a H with the causally computed a R . Like in the case of test data, the causally computed amplitude a R is a better estimation of the envelope than a H is. Moreover, computation of a R is stable with respect to outliers (Fig. 3d). Next, in Fig. 3e,f we show the Hilbert phase ϕ H along with two causally-obtained phases ϕ L and ϕ R . We see that for large-amplitude oscillation, all the phases practically coincide. When the amplitude almost vanishes, the "device" of the locking-based method falls out of synchrony and makes one cycle less. On the contrary, the resonant-oscillator phase ϕ R and the Hilbert phase reveal the same number of cycles. We emphasize, that both causal techniques are not sensitive to the outliers in the original data, while the Hilbert approach is (see Fig. 3f). Finally, we mention that the non-resonant oscillator technique works poorly for the unfiltered tremor data. Only if we use the bandpass filter 4.5 ± 2 Hz, then this technique works perfectly.

Wide-band signal: beta-band brain activity
Elevated beta-band activity has been established as a marker for rigidity and bradykinesia in patients with Parkinson's Disease 23 . As such, it has been employed in several studies of adaptive deep brain stimulation aiming for automatic adjustment of 4/11 stimulation parameters in response to beta-band amplitude 7,9 . Here we recorded LFP data from a patient with Parkinson's Disease on medication, 2 days after surgery for deep brain stimulation with electrode cables being externalized. LFPs were acquired from the left subthalamic nucleus via a 4-contact electrode (Model 3389, Medtronic, Minneapolis, USA) using a D360 amplifier (Digitimer Ltd., Welwyn Garden City, UK) and were digitized at 1 kHz sampling rate with a 1401 analog-to-digital converter (CED Ltd., Cambridge, UK). A bipolar referencing scheme was adopted, with the signal shown here originating from contacts one to two. During the recording session, the patient was comfortably seated and was asked to rest quietly.
The analysis of the beta-band brain activity requires a bandpass filter preprocessing. We use a simple FIR filter with bandwidth 17 ± 4 Hz (281-point filter was generated by the Matlab fir1 function). A causal filtration introduces a delay, but this is a necessary price to be paid. Our tests with this signal demonstrate that the non-resonant-oscillator technique outperforms the other two that fail because of the substantial variation of the signal's amplitude. Noteworthy, due to the bandpass, the approach works without the baseline correction and the frequency adaptation. Hence, the algorithm implementation requires only a few lines of code. The results are illustrated by Fig. 4.

Discussion
To summarize, in this paper, we considered the problem of causal instantaneous amplitude and frequency estimation. Contrary to attempts to adjust the inherently non-causal analytic signal approach, we used oscillation theory ideas. We presented and tested three algorithms. The first one exploits the phase-locking property of nonlinear limit-cycle oscillators, while the other two rely on linear resonance. Technically, we encountered solving differential equations incorporating experimental data; we suggested efficient numerical schemes to solve them. Below, we discuss the advantages and disadvantages of all three techniques. Since these techniques aim at real-time applications, we pay special attention to their computational efficiency. The locking-based technique is similar to the phase-locked loop approach and provides the phase only. It works well for relatively narrow-band signals like the tremor data and does not require an additional bandpass filter. It can easily incorporate the frequency adaptation that makes the algorithm able to cope with the signal's frequency's slow variation. The technique does not work well with the signals like beta-band activity, where both the frequency and the amplitude vary relatively fast; this results in frequent loss of synchrony. The algorithm is computationally efficient: the only demanding operations for a new phase value are several computations of the sine function.
As an advantage of the resonant oscillator technique, we mention that it provides both the phase and the amplitude. It is stable to high-frequency noise because the resonant oscillator acts as a weak bandpass filter. However, it requires baseline correction because the integrator unit enhances the low-frequency perturbations. The implementation described in the Methods section is very efficient, provided the oscillator's frequency does not change. Adaptation of the frequency to that of the signal requires recomputation of algorithm's coefficients what slightly reduces the efficiency.
If extraction of the rhythm of our interest requires a bandpass filter, the non-resonant oscillator approach is the best choice. Using the beta-band brain activity as an illustrative example, we demonstrated that this causal algorithm provides phase and amplitude very close to non-causally computed Hilbert Transform-based values. An essential advantage of the algorithm is its efficiency: computation of the new phase and amplitude values requires only several arithmetic operations and two function calculations. This property makes the algorithm especially appropriate for the real-time processing of high-frequency signals.
We emphasize that all three techniques are much faster than approaches based on the HT. For example, the technique introduced in 17 requires the computation of the direct and inverse Fourier transform in a running window plus additional filtration to compensate for the boundary effects and non-causal nature of the HT; altogether, this computation is about thirty times more demanding than an oscillator-based approach.
Taken together, the techniques suggested here may provide useful means for real-time detection of phase and amplitude in the context of adaptive deep brain stimulation. For example, estimation of instantaneous amplitude has been implemented by interpolating beta-filtered and squared LFP data 24 . However, this approach is not able to recover the signal's phase simultaneously. Additionally tracking phase may augment therapeutic benefit, though, by delivering stimulation pulses at specific points within the oscillatory cycle. This has been shown to effectively disrupt pathological oscillations [10][11][12]17 . So far, to our knowledge, adaptive stimulation algorithms integrating both the phase and amplitude of a signal have not been established yet.
Finally, we mention that our techniques do not cause any additional delay. The processing delay is only due to a causal filter if the latter is required, e.g., for brain activity data. To minimize the delay, one has to exploit specially designed filters 25 .

Methods
In all cases the input is an oscillatory signal s(t) is sampled with interval ∆. Thus, available is a sequence s(t k ) = s(k∆) = s k .
The non-causal Hilbert-Transform based approach. In this approach one constructs the corresponding complex-valued analytic signal Z(t) = s(t) + iŝ(t), whereŝ(t) = π −1 P.V. ∞ −∞ s(τ) t−τ dτ is the HT of s(t). Obviously, the HT is a non-causal operation; HT of a finite-length time series yields spurious values for the boundaries. The absolute value and argument of Z(t) provide the instantaneous Hilbert amplitude a H (t) and phase ϕ H (t), respectively. For a narrow-band one-component signals like s(t) = A(t) cos[ω(t)t], where A(t), ω(t) are slow functions of time, the analytic signal approach provides a H ≈ A and ϕ H ≈ ω(t)dt. For a discussion of the HT's practical implementation and technical hints we refer, e.g., to 1,4,26 . In a practical implementation either a discrete evaluation of the integral is performed, or a discrete Fourier transform is used.

Causal estimation of phase and amplitude
In these methods, we obtain the value of the instantaneous phase ϕ(t k ) = ϕ k and amplitude a(t k ) = a k by using only the current and the previous values of the signal, i.e., s k , s k−1 , . . ..

Measuring "device": phase-locked oscillator
The synchronization theory says that an oscillatory force s(t) acting on a limit-cycle oscillator can entrain it, if the frequency of the force is close to the natural frequency of the limit-cycle oscillator. It means that the oscillator's frequency becomes equal to that of the force, and their phases fulfill the locking condition ϕ − θ ≈ const, where θ and ϕ are the oscillator's and the signal's phases. For our purposes, it is appropriate to use the so-called phase oscillator. Its forced dynamics is described bẏ where ε is a parameter that determines the coupling strength. Consider first a harmonic force, s = a cos(νt). If ω = ν, then in the locked state ϕ ≈ θ . However, ν is not known a priori, but we assume that we can roughly estimate it. Let this initial guess be ν 0 . Thus, we set initially ω = ν 0 and start our computation with this value. To adapt the phase oscillator to the signal, we estimate the frequency of the forced oscillator ν e on the fly. (The index e stands for "estimated"). To this end, for the time instant t, we take previously computed unwrapped phases θ for the time interval [t − T e ,t], where T e ∼ 2π/ν 0 (approximately 6/11 one or two cycles). Assuming that within this time interval θ (t ) = θ (t − T e ) + ν e (t − t + T e ), we compute the frequency ν e via the linear fit. Next, we update the oscillator's autonomous frequency as where K is a constant update factor. Adapting the frequency in this way, we ensure a transition from the unlocked to locked dynamics. Performing frequency estimation several times per cycle, we successfully estimate the phase of signals with slowly drifting frequency. We expect that this algorithm also works if the force's amplitude a slightly varies with time (so that the oscillator still remains locked). One can treat the suggested scheme as a software implementation of a phase-locked loop 22 , cf. also 27 . Practically, we have to solve the differential Eq. (2) numerically, whereas its right-hand side is known only in discrete time points t k . The easiest way is to exploit Euler's technique with the integration step ∆ to advance from the known value θ k at t k to the new phase θ k+1 at t k+1 = t k + ∆, given a new measurement s k+1 . This technique may work properly if ∆ is sufficiently small; otherwise, the numerical solution becomes unstable. (We recall that ∆ is the fixed sampling interval that cannot be made arbitrarily small). Therefore, to advance the solution of Eq. (2) from ϕ k to ϕ k+1 given new measurement s k+1 , we use the parabolic approximation of s(t) on the interval [t k ,t k+1 ] and then exploit the standard Runge-Kutta algorithm. (Notice that the Runge-Kutta step can be much smaller than ∆.) The coefficients of the parabolic fit are computed from s k−1 , s k , s k+1 . Starting with some initial condition, e.g., θ 0 = 0, we achieve, after a short transient (see examples in Fig. 1), the synchronous state where ϕ k ≈ θ k . Notice that the phase oscillator can be complimented by a low-pass filter: then the scheme becomes the simplest traditional phase-locked loop 4,22 ; for τ → 0 it reduces to Eq. (2). This extension may be useful for the data with strong high-frequency noise but it does not show any advantages for the tremor data we use.

Measuring "device": non-resonant linear oscillator
Our second measuring "device" is linear damped oscillator: where ω and α are the frequency and the damping parameter of the oscillator, respectively. Suppose first that s(t) is harmonic, s(t) = a cos(νt) = a cos(ϕ(t)). The well-known stationary solution of linear Eq. (4) is (ω 2 −ν 2 ) 2 +(αν) 2 and β = arctan −αν ω 2 −ν 2 . Thus, the forced system (4) oscillates with the force's frequency ν and the amplitude b. The dependencies of the amplitude ratio b/a and of the phase shift β on the forcing frequency ν reflect the well-known resonance effect. Knowing the oscillator's state x(t) = b cos(νt + β ),ẋ(t) = −bν sin(νt + β ) we find the amplitude and phase of the external force: Thus, if our "device" yields x(t) andẋ(t), then we easily compute a(t) and ϕ(t), provided the frequency ν of the force is known. However, it is not known but can only be roughly estimated. Moreover, generally, it varies with time. Before we discuss how to cope with this fact, we mention how we practically obtain x(t) andẋ(t). Here, we make use of the linearity of Eq. (4) and develop and exploit an efficient numerical scheme, presented below. With this scheme, starting with some initial conditions, e.g., x 0 =ẋ 0 = 0, we obtain, after a short transient, the forced solution of Eq. (4) and compute a(t k ) and ϕ(t k ) from Eqs. (5). Now, we discuss how to choose the parameters of the measuring oscillator. Recall that to compute a, ϕ, we need the value of the signal frequency ν. First, let us consider the amplitude measurement. Inspecting Fig. 5, we see that if we take a large value of α (strongly damped oscillator) and ν ω then the ratio b/a is practically independent on ν. Thus, to compute the second square root in Eq. (5) we need only a very rough estimate of ν. For the phase estimation, the damping parameter shall be different. Indeed, as follows from Fig. 5b, now we have to choose α to be small, then the phase shift β ≈ 0 in a wide range of ν and the phase of the external force equals the phase of the oscillator. A reasonable choice is to take the oscillator's frequency ω about five times larger than ν. Thus, we neglect β in the expression for ϕ and compute (ω 2 − ν 2 ) 2 + (αν) 2 in the expression for a only once, using an initial guess ν 0 for frequency. However, the terms x(t) 2 + [ẋ(t)/ν] 2 and arctan −ẋ(t) νx(t) essentially depend on ν. Imprecise estimation of ν and/or its variation with time results in spurious oscillations of the estimated amplitude and phase. To remove these oscillations and to make the results less dependent on the initial frequency estimation, we improve this estimation by computing the signal frequency in real-time as discussed in the precious section. We start with an initial guess  Figure 5. Resonance curves for amplitude (left) and phase shift (right) for the linear oscillator with frequency ω = 5. Blue and red curves correspond to the weakly (α = 0.2) and strongly (α = 6) damped oscillator, used for the phase and amplitude measurement, respectively (these parameters were used to process the artificial data). The domains where b/a ≈ const (for the red curve) and β ≈ 0 (for the blue curve) are marked by red and blue arrows, respectively. These are intervals of signal frequency ν where the algorithm's performance is good. We see that ω ≈ 5ν is a reasonable choice.
ν 0 and begin the computation of a(t), ϕ(t). After a transient time of the order of several cycles, we begin with the frequency estimation to obtain ν e . (The estimation is performed in exactly same way as for the phase-locked oscillator.) We use ν e instead of ν 0 and, in this way, significantly improve the real-time computation of the amplitude and frequency. Performing frequency estimation several times per cycle, we track the slowly drifting signal frequency. (Notice that for the bandpass filtered signal used in the last example, the algorithms works well without any frequency correction. Here we set nu equal to the center of the bandpass.) In summary, for the amplitude measurement, we need a strongly damped oscillator, while for the phase estimation, we need an oscillator with small damping. If we need both measurements simultaneously, then the solution is to exploit two "devices" concurrently. Though Eqs. (5) are derived for the harmonic force, we expect that they approximately hold for signals with a slow variation of the amplitude and frequency. The presented numerical tests confirm this expectation.

Measuring "device": resonant linear oscillator
Here we adopt the technique used for model studies in our previous publications 18,19 . The device consists of a linear oscillator in resonance with the measured signal and an integrating unit: µż + z =ẋ .
The role of the harmonic oscillator Eq. (6) is twofold. First, it acts as a bandpass filter and the damping factor α determines the width δ f ≈ α/2π of the bandpass. Second, the harmonic oscillator yields signalẋ which phase is close to that of the input x(t), provided the frequency ω is close to the mean frequency of s(t). (This condition also ensures that the band-pass is centered at the frequency of the signal.) Next, for µ 1 Eq. (7) acts as the integrating unit. Its output is shifted by π/2 with respect toẋ. It is useful to rescaleẋ, z so that their amplitudes are close to that of s(t). For this, we compute u = αẋ and w = αω 0 µz and obtain the instantaneous phase and amplitude of s(t) as ϕ(t) = arctan(w/u), a(t) = √ u 2 + w 2 .

Solving the linear oscillator equation for discrete input signal
Here we present the numerical scheme for solving Eq. (4), adopted to the situation where the input signal is available at a finite sampling rate. Using the substitution y = xe γt , where γ = α/2, one obtains y + η 2 y = s(t)e γt = S(t) , with η 2 = ω 2 − γ 2 .
Another standard substitution where A * denotes complex conjugate of A, yieldsȦ = −iη −1 S(t)e −iηt . Integrating we obtain Next, locally interpolating the measured signal s(t) by a parabola going through the points s k−1 , s k , s k+1 we compute the integral I = ∆ 0 S(t k + τ)e −iητ dτ. As a result, we obtain the following practical scheme for integration of Eq. (4). First, we pre-compute 8/11 the coefficients C 1,2,3 : where Next, for given x k ,ẋ k and the new measurement s k+1 we compute x k+1 ,ẋ k+1 in three steps: We emphasize that all coefficients can be pre-computed, so that the integration step requires only a few summations and multiplications.
To summarize the presentation of numerical techniques, we emphasize that if the sampling rate is very high, the differential equations in all presented algorithms can be fast and easily solved by the midpoint or predictor-corrector technique 28 . However, the stability of these simple methods is not guaranteed.

Choosing parameters
Locking-based oscillator. The main parameter is the forcing coefficient ε (Eq. 2). The larger ε the faster the system synchronizes. On the other hand, the forcing term ε sin θ · s(t) shall be small enough to ensure the monotonicity of the phase θ . So, if s(t) = a cos ϕ then (Eq. 2) readsθ = ω + εa 2 sin(ϕ − θ ) − εa 2 sin(ϕ + θ ); in the locked state ϕ ≈ θ the first sine term vanishes andθ > 0 if εa < 2ω. In the artificial data example we checked that he values 0.4 ≤ ε ≤ 0.8 provide good result. In this example, the amplitude goes up to ≈ 2 and the frequency is about 1, so that the above condition is satisfied. The parameter of the frequency adaptation shall be K ≤ 1; the values between 0.5 and 1 work well. The data in Fig. 5 correspond to ε = 0.8, K = 1. For the tremor data ω ≈ 2π · 4.5, ε = 30, and K = 0.5. In both cases we updated the frequency ω 20 times per period using the preceding phases θ in the interval that is about one cycle long.

Resonant oscillator.
Here we always take α = 0.3ω what corresponds to bandpass width about 30% of the central frequency ω. The requirement for the parameter µ is µ 1; in all computations we choose µ = 500.