{'fraction_non_alphanumeric': 0.03698715292167427, 'fraction_numerical': 0.00905511811023622, 'mean_word_length': 4.8131775475788965, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Complex-valued Gaussian Process Regression for Time Series Analysis

The construction of synthetic complex-valued signals from real-valued observations is an important step in many time series analysis techniques. The most widely used approach is based on the Hilbert transform, which maps the real-valued signal into its quadrature component. In this paper, we define a probabilistic generalization of this approach. We model the observable real-valued signal as the real part of a latent complex-valued Gaussian process. In order to obtain the appropriate statistical relationship between its real and imaginary parts, we define two new classes of complex-valued covariance functions. Through an analysis of simulated chirplets and stochastic oscillations, we show that the resulting Gaussian process complex-valued signal provides a better estimate of the instantaneous amplitude and frequency than the established approaches. Furthermore, the complex-valued Gaussian process regression allows to incorporate prior information about the structure in signal and noise and thereby to tailor the analysis to the features of the signal. As a example, we analyze the non-stationary dynamics of brain oscillations in the alpha band, as measured using magneto-encephalography.


Introduction
In this paper, we introduce the use of complex-valued Gaussian process (GP) regression for the construction of complex-valued signals. The principal aim of this construction is to quantify the instantaneous amplitude, phase and frequency of real-valued oscillatory signals. GP regression is a powerful tool for the analysis of time series that relies on a prior model of their temporal autocorrelations, as specified by a covariance function (Rasmussen, 2006;Brahim-Belhouari and Bermak, 2004;Reece and Roberts, 2010). We define two new families of complex-valued covariance functions that induce the appropriate crosscorrelation between the real and the imaginary part of the complex-valued signal. This statistical relation is then exploited for constructing the (unobservable) imaginary part from the (observable) real part of the signal. The resulting complex-valued GP regression (CGPR) is a very flexible tool that can be adapted to a wide range of signals by choosing a particular covariance function that is tailored to the specific properties of the signal.
In the following, we show that the complex-valued signal obtained by CGPR can greatly outperform existing methods in estimating the instantaneous amplitude and frequency of an oscillatory signal. See (Boashash, 1992b,c) for a review of the most commonly used methods for instantaneous frequency estimation). The improvement is greatest for highly non-stationary signals such as chirplets. We illustrate the method by analyzing real magnetoencephalographic (MEG) recordings from a human participant. We focus on the analysis of alpha oscillations, a brain rythm that has been associated to attention and perception (Kelly et al., 2006;Jensen and Mazaheri, 2010;van Ede et al., 2012van Ede et al., , 2011.

Methods
In this section, we introduce the basic concepts of CGPR and its relevance for time series analysis. Our main goal is to improve the estimation of the instantaneous frequency and amplitude of a real-valued signal. The most commonly used method for instantaneous frequency and amplitude estimation involves the complex-valued analytic representation of the measured real-valued signal (Boashash, 1992b). Therefore, we start the exposition by reviewing the construction of the analytic representation. Then, after a short review of the GP regression framework, we introduce the CGPR and the construction of the associated complex-valued signal. The properties of the resulting complex-valued signal depend on the choice of the covariance function. We introduce two new classes of covariance functions that reproduce the properties of the analytic representation in expectation. We denote these classes of functions as quadrature and quasi-quadrature covariance functions.

Theoretical background
In this subsection, we will briefly review some theoretical results that are relevant for the understanding the rest of the paper. In particular, we will introduce the quadrature filter as a tool for extracting the envelope and phase of a real-valued band limited signal.

Instantaneous amplitude and phase
Given a complex-valued signal z(t), its instantaneous amplitude and phase are defined through its polar representation: (1) Where the complex modulus |z(t)| = (Rz(t)) 2 + (Iz(t)) 2 is called instantaneous amplitude and arg[z(t)] = −i log z(t) |z(t)| is called instantaneous phase. Furthermore, the instantaneous (angular) frequency ω(t) is defined as the first derivative of the instantaneous phase (2) In order to extend these definitions to real-valued signals, it is natural to replace the complex exponential in Eq. 1 with a real-valued sinusoid: where A(t) is a positive-valued amplitude function and P(t) is a phase function. Unfortunately there are infinitely many ways of representing an arbitrary real-valued signal in this form. In other words, the problem of decomposing a real-valued signal into the product of an amplitude function and an oscillation with varying frequency is not well posed. Nevertheless, when the real-valued signal is band-limited, it is often intuitive and insightful to consider it as an oscillation subject to amplitude and frequency fluctuations.

Analysis of a simple family of band-limited signals
We will now see how to construct the instantaneous amplitude and phase for a simple family of real-valued signals. In doing so, we will follow a general strategy that will also be useful in the case of arbitrary band-limited signals. The idea is to define a transform that converts the original signal into a complex-valued signal whose instantaneous amplitude and phase can be extracted using Eq. 3. Consider the following deterministic signal where A(t) is an arbitrary positive-valued function that has most of its energy at frequencies lower than f 0 . It is intuitive to consider s(t) as an oscillation with instantaneous amplitude A(t) and instantaneous phase 2πf 0 t. Hence, we aim to construct a complex-valued signal z(t) such that |z(t)| = A(t) and arg[z(t)] = 2πf 0 t. This can be done by applying the formula cos(a) = 1 2 e ia + e −ia from which follows that s(t) decomposes in the sum of a counterclockwise part s + (t) and a clockwise part s − (t): where From this expression it is clear that the complex-valued signal z(t) constructed from s(t) should be its counterclockwise part. In fact: and arg[s + (t)] = 2πf 0 t .

Quadrature filter and analytic representation
Most real-world signals do not have a simple analytic expression such as Eq. 4. Nevertheless, it is possible to devise general purpose transforms that map arbitrary band limited realvalued signals to complex-valued signals that give an intuitively plausible quantification of instantaneous amplitude and phase. The most commonly used transform of this type is the quadrature filter (Boashash, 1992a). The action of the quadrature filter A on a real-valued signal s(t) is defined by the following formula where is called the Hilbert transform of s(t). Here and in the following, expressions such as Hs(t) denote the application of the linear operator H on the signal s(t). The filter kernel of the Hilbert transform is the function 1 πt and, by taking its formal Fourier transform, we obtain the frequency response function of the Hilbert transform: Hence, the Hilbert transform acts on the signal by phase-shifting the Fourier coefficients for positive and negative frequencies by, respectively, −π/2 and π/2. In the rest of the paper, we will refer to this relation between the Fourier coefficients of the real and the imaginary part as quadrature relation. Using Eq. 11, it is easy to show that the frequency response function of the quadrature filter is where h(ξ) is the Heaviside step function. This means that the transformed signal As(t) is always analytic, meaning that it has all its Fourier coefficients relative to the negative frequencies equal to zero. Therefore, the quadrature filter maps a signal into its unique analytic representation, defined as the complex-valued signal that has its real part identical to the original signal while having all the Fourier coefficients relative to the negative frequencies equal to zero (Boashash, 1992a). This mapping does not entail any information loss since, for real-valued signals, the Fourier coefficients of the negative frequencies are the complex conjugate of their positive frequency counterparts.

Bedrosian's theorem
Bedrosian's theorem (Wang, 2009) shows that the quadrature filter can perfectly separate the envelope function from the oscillatory part when their spectra are disjoint. More formally, consider a signal of the form given by Eq. 3. If there exists a positive number a such that the spectrum of A(t) is identically equal to zero for |ξ| > a while the spectrum of cos P(t) is identically equal to zero for |ξ| < a, then we have that This theorem justifies the use of the analytic representation for the quantification of the instantaneous amplitude and phase of a real-valued signal.
In the following, we will give an intuitive proof of Bedrosian's theorem for the special case of a signal of the form given by Eq. 4. From the previous analysis we know that, in order to correctly recover A(t) and P(t), the complex-valued signal constructed from s(t) should be its counterclockwise part s + (t). Therefore, we need to show that, when the conditions of Bedrosian's theorem are met, As(t) is equal to s + (t). This can be done by taking the Fourier transform of s(t). Using the convolution theorem, together with the fact that the Fourier transform of the complex exponential e iyt is the delta function δ(ξ − y), it is easy to see that this Fourier transform is . A) Spectrum of the envelope (blue) and oscillation (green). B) Signal in the frequency domain, the counterclockwise (clockwise) part has energy only in the positive (negative) side of the spectrum. wheres are the Fourier transforms of the counterclockwise part and of the clockwise part respectively (see Eq. 5). Recall that, in the frequency domain, the analytic representationÃs(ξ) is obtained by setting to zero all the Fourier coefficients relative to the negative frequencies and multiplying by two. Hence, it is easy to see thatÃs(ξ) is identically equal tos + (ξ) =Ã(ξ − 2πf 0 ) if and only ifÃ(ξ) vanishes when |ξ| > 2πf 0 . In fact, in this case, the counterclockwise and the clockwise parts of the signal have energy only at the positive and negative side of the spectrum respectively and, consequently, they can be perfectly separated by the quadrature filter. Figure 1 shows the spectral separation of the counterclockwise and clockwise parts for an example signal (in the frequency domain). Panel A shows the spectrum of the envelope (in blue) and of the oscillation (in green), note that both envelope and oscillation have energy only in a finite band of the spectrum. Panel B shows that in this example the counterclockwise (clockwise) part has energy only in the positive (negative) side of the spectrum.

Limitations of the quadrature filter
The Bedrosian's theorem is important for understanding the limitations of the quadrature filter. In particular, the condition of disjoint spectra is too restrictive for most real life signals and can lead to rather counterintuitive estimates of the instantaneous amplitude and phase. For example, consider the signal s(t) = e − |t| σ cos (2πf 0 t) where σ is a positive constant. This signal (with σ = 0.15 and f 0 = 3) is plotted in Fig.2A and its spectrum is plotted in Fig.2B. Importantly, the spectra of e − |t| σ and cos (2πf 0 t) are not disjoint because the Fourier transform of e − |t| σ is different from zero for all frequencies. This implies that the counterclockwise part has energy at negative frequencies and it is not equal to the analytic representation of the signal (see Fig.2C). Consequently, the analytic representation does not recover A(t) and P(t) exactly. However, the spectral overlap is very small if σ 1/f 0 and, under this condition, the estimates |As(t)| and arg[As(t)] are very close to respectively A(t) and P(t). Conversely when σ ≈ 1/f 0 , the estimates can differ dramatically from A(t) and P(t). Fig. 2D shows the amplitude constructed using the analytic representation (σ = 0.15, f 0 = 3). Note that this function is not unimodal and has oscillating tails.
More dramatic examples of this kind of behavior can be found in the Results section. In the following subsections, we will use CGPR to derive a class of probabilistic filters that improve on this sub-optimal behavior of the analytic representation.

Real-valued Gaussian process regression
We now review the basic concepts of real-valued GP regression, laying the groundwork for the definition of the GP-based complex signal. In fact, the complex-valued extension easily follows from the real-valued case and will allow us to obtain a probabilistic construction of the complex-valued signal.
A GP generalizes the multivariate Gaussian distribution to infinitely many degrees of freedom. More formally, a stochastic process ψ(t) follows a GP distribution when all its marginal distributions p(ψ(t 1 ), ..., ψ(t n )) are multivariate Gaussians (Rasmussen, 2006). Since these Gaussians can be parametrized in terms of their mean and covariance matrices, a GP is fully specified by its mean and covariance functions, respectively and In order to properly specify a GP distribution, k ψ (t 1 , t 2 ) has to be symmetric and nonnegative definite, meaning that all the covariance matrices of the marginal distributions are symmetric and non-negative definite.
In the remaining, we assume the mean function m ψ (t) to be zero for all values of t. Furthermore, we restrict our attention to stationary covariance functions, i.e. functions that solely depend on the time lag k(t 1 , t 2 ) = k(t 2 − t 1 ) = k(τ ). In this case, we can define the spectral density (or power spectrum) ofk(ξ) as the Fourier transform of k(τ ) Time ( The usefulness of GPs in machine learning follows from the fact that GP distributions can be used as priors for a non-parametric Bayesian regression. In particular, given a set of n sample points T = {t 1 , ....t n }, we can model a sampled signal s j as follows: where the vector , whose j-th entry is j , is a multivariate noise term that we assume to be Gaussian with covariance function Q. Note that the sample point do not need to be equally spaced. The aim of this regression problem is to estimate the latent function ψ(t) without resorting to a specific parametric form (e.g. linear, polynomial, sinusoidal, wavelet). Using GP regression, we can estimate the value of ψ(t) at any arbitrary set of m target points T × = {t × 1 , ..., t × m }, even if they are not part of our sample points. To this end, it is convenient to define the vector ψ with entries ψ j = ψ(t × j ). Since both the prior distribution of ψ and the likelihood p(s|ψ) are multivariate Gaussians, the posterior distribution p(ψ|s) is also a multivariate Gaussian. In particular, its expected value m ψ|s is given by: where the sample covariance matrix K ψ is defined by the entries mussen, 2006). From this formula we can see that the matrix K × ψ projects the information from the sample points to the target points by leveraging our prior model of the temporal correlations determined by the function k ψ (t 1 , t 2 ).
A similar formula applies if we aim to reconstruct a completely unobservable new process ρ(t) from the observation of ψ(t) once we assume their prior cross-covariance function c ρ,ψ (t 1 , t 2 ) = ρ(t 1 )ψ(t 2 ) . This point is crucial for our current goals as it will allow us to reconstruct the unobservable imaginary part from a measured real-valued signal. If we restrict our attention to the sampled time points, the posterior expectation of ρ given s is where the j, k-th entry of the cross-covariance matrix C ρ,ψ is c ρ,ψ (t j , t k ). We conclude this subsection by discussing the computational complexity of GP regression. The computational bottleneck of Eq. 20 and Eq. 21 is the inversion of the covariance matrix. The invesion of a positive-definite symmetric matrix is usually performed using the Cholesky decomposition. This decomposition of the covariance matrix has a cubic complexity with respect to the number of time points (Krishnamoorthy and Menon, 2011).

Complex-valued Gaussian Process regression
We can now generalize GP regression to the case in which the latent process ζ(t) = α(t) + iβ(t) is complex-valued. In order to construct a complex-valued GP we assume all the marginals of ζ(t) to be circularly-symmetric complex normal distributions (Picinbono, 1996). The resulting stochastic process is fully specified by the Hermitian, non-negative definite function k ζ (t 1 , t 2 ) = ζ(t 1 )ζ(t 2 ) , whereζ(t 2 ) is the complex conjugate of ζ(t 1 ). In analogy with the real-valued case, we refer to k ζ (t 1 , t 2 ) as the covariance function of the process.
In the remaining of the subsection we will derive the posterior expectation of the complex-valued signal ζ(t) for an arbitrary Hermitian covariance function k ζ (t 1 , t 2 ). We must assume that only the real part of ζ(t) generates s whereas the imaginary part is not observable Here, we assume that the measurements are corrupted by a noise term j . This is an important difference with the usual construction of the analytic representation, which cannot account for measurement noise. In fact, as already mentioned in the previous section, the noise has to be removed prior to the application of the quadrature filter. This two-step procedure ignores the features of the signal when removing the noise. Conversely, the GP approach can leverage the prior model of signal and noise, as expressed by their covariance functions.
The posterior expectation of ζ(t j ) can be obtained using Eq. 21. In fact, since the complex distribution is circularly symmetric, the prior cross-covariance matrix between the real and the imaginary part at the sample points is given by IK ζ , which denotes the entrywise imaginary part of the matrix K ζ (Picinbono, 1996). Analogously, the prior covariance matrix of both the real and the imaginary part is RK ζ , i.e. the entrywise real part of K ζ . Thus, from Eq. 21, it is easy to see that the posterior expectation of the complex-valued signal ζ(t) at the sample points is

Quadrature covariance functions
Our construction of a GP-based complex signal depends on the choice of the prior covariance function. The aim of this subsection is to introduce a family of covariance functions that are suitable for a probabilistic generalization of the analytic representation. One of the defining features of the analytic representation is that, for any given frequency, the Fourier coefficients of the imaginary part are obtained by a ± π 2 phase shift of the coefficients of the real part (see Eq. 11). Our aim is to reinterpret this relation in a probabilistic sense. In particular, we say that the process whereα(ξ) andβ(ξ) are the Fourier transforms of, respectively, the real and the imaginary part, and β (ξ) α(ξ) denotes the conditional expectation ofβ(ξ) givenα(ξ). In this expression, the positive proportionality constant c(ξ) ranges from zero to one and represent the strength of the statistical association between the Fourier coefficients of the real and the imaginary part. In order to construct a GP that is in quadrature relation in expectation, we define the quadrature covariance function induced by an arbitrary (real-valued) stationary covariance function k(τ ) as follows: This formula defines a valid covariance function, i.e. a function that is Hermitian and nonnegative definite (see appendix I). We will now show that the Fourier coefficients of any GP defined by a covariance function of this form are indeed quadrature related in expectation.
To this end, it is convenient to organize the real and the imaginary part of the process ζ(t) into the real-valued vector process Since ζ(t) is governed by a circularly-symmetric GP with quadrature covariance function, the (matrix valued) cross-covariance function of Ψ(t) is given by the following formula (see appendix II for a derivation) We can obtain the cross-spectral densityΘ(ξ) of Ψ(t) by taking the Fourier transform of Θ(t) using Eq. 11Θ The conditional expectation β (ξ) α(ξ) can be obtained from the cross-spectal density as follows (see (Picinbono, 1996)): where the coherency γ(ξ) is given by Therefore, the real and the imaginary part of the signal are quadrature related in expectation if and only if the phase of the coherency arg[γ(ξ)] is equal to −sgn(ξ) π 2 . For a cross-spectral density of the form given by Eq. 28, the coherency is Using this result and the fact thatΘ 11 (ξ) is equal toΘ 22 (ξ), it follows from Eq. 29 that n the Fourier coefficients of the imaginary part is Hence, the Fourier coefficients of the real and the imaginary part are quadrature related in expectation with c(ξ) = |γ(ξ)| = 1. So far, we have shown that the CGPR with a quadrature covariance function allows to construct a probabilistic version of the analytic representation. However, this result does not fully realize our aim since the correlation between the real and imaginary part is perfect and, consequently, the samples from the stochastic process are analytic functions with probability one. In order to understand the origin of this perfect correlation, we will now study the eigenvalues and eigenvectors of the cross-spectral density. We denote the two eigenvectors ofΘ(ξ) as v (1) (ξ) and v (2) (ξ) with eigenvalues λ (1) (ξ) and λ (2) (ξ) respectively. Note that, since the matrix is Hermitian and non-negative definite, the eigenvalues are real and non-negative. For quadrature covariance functions, the eigenvectors are v (1) = (i, 1) and v (2) = (−i, 1) (for all values of ξ), and their corresponding eigenvalues are, respectively, λ (1) (ξ) = h(ξ)k(ξ) and λ (2) (ξ) = h(−ξ)k(ξ), where h(ξ) is the Heaviside step function. Since the eigenvalues are never simultaneously non-zero, except for the irrelevant single point ξ = 0, all the bivariate complex distributions ofα(ξ) andβ(ξ) are degenerate and the quadrature relation is deterministic. Consequently, the resulting complex-valued signal does not alleviate the problems of the analytic representation in dealing with signals with overlapping amplitude and phase spectra.
We will conclude this subsection by introducing two quadrature covariance functions that will be used in the remaining of the paper. One of them, the periodic covariance function, is most suitable for the analysis of noise-free signals. The other, the harmonic covariance function, should be preferred in the analysis of noise corrupted signals. For a real-valued signal, the periodic covariance function is given by the following formula where ω 0 is the center frequency and ρ regulates the smoothness of the waveform. The associated quadrature covariance function can be obtained using Eq. 25. Fig. 3A shows the real and the imaginary part of the periodic quadrature covariance function (ω 0 = 2π × 3Hz, ρ = 0.5). The imaginary part has been obtained from a numerical approximation to the Hilbert transform of the real part). Fig. 3B shows the spectral density of this covariance function. Since the function is periodic, the density concentrates on a sequence of discrete frequencies. Importantly, the spectral density is zero for negative frequencies, implying that this covariance function is analytic. The periodic covariance function can capture complex waveforms but is sensitive to noise in the data since it does not suppress high order harmonics (see the spectrum in Fig. 3B). This is problematic as the instantaneous phase can be meaningfully quantified only when the signal is band limited. Consequently, in the analysis of noisy signals, we will use the harmonic quadrature covariance function, defined as follow This covariance function is much less flexible than 33 but it has the advantage of filtering out the high order harmonics and, consequently, being less prone to be affected by the noise in the signal.

Quasi-quadrature covariance functions
We will now show how to construct a non-degenerate covariance function that generates a process whose real and imaginary parts are quadrature related in expectation. More formally, we want to construct a complex-valued Gaussian process that satisfies Eq. 24 with c(ξ) < 1. To this end, we have to modify the quadrature covariance function such that both eigenvalues differ from zero for all frequencies whereas the phase of the coherency remains equal to −sgn(ξ) π 2 (see Eq. 29)). (ω 0 = 2π × 3Hz, ρ = 0.5). The imaginary part was obtained by numerical Hilbert transform. B) Spectral density of the periodic quadrature covariance function. C) Real part (blue line) and imaginary part (green line) of the periodic quasiquadrature covariance function (ω 0 = 2π×3Hz, ρ = 0.5) with squared exponential envelope (δ = 0.2s). D) Spectral density of the periodic quasi-quadrature covariance function with squared exponential envelope.
We define quasi-quadrature covariance functions as the product between a quadrature covariance function and a positive-valued envelope function f (τ ). More formally, we define the quasi quadrature operator A f as follows: For reasons that will be clear later, the envelope function has to be even and positive definite. By positive definite, we mean that its Fourier transform is positive-valued, so that the result is a valid covariance function. Furthermore, we assume f (τ ) to be smooth and rapidly decreasing, i.e. it tends to zero faster than any inverse power when τ tends to ±∞. This assumption guaranties that the Fourier transform of f (τ ) is itself rapidly decreasing.
We will now show that the eigenvalues of quasi quadrature covariance functions are simultaneously different from zero. To this end, we compute the cross-spectral densitỹ Θ f (ξ) of the new process by applying the Fourier convolution theorem on its cross-covariance This is an integral of matrices that share the same eigenvectors. Consequently, the eigenvectors ofΘ f (ξ) are still v (1) and v (2) for all frequencies but the eigenvalues are now given by Therefore, the first and the second eigenvalue functions are given by the convolution of, respectively, the analytic (k(ξ)h(+ξ)) and the anti-analytic (k(ξ)h(−ξ)) part of the spectrumk(ξ) with the Fourier transform of f (t). Note that, since f (t) is positive definite, f (ξ) is strictly positive and both eigenvalues are never exactly zero. Consequently, the cross-spectral density is never degenerate, and the Fourier coefficients of the real and the imaginary part are not perfectly correlated. We now need to show that quasi-quadrature covariance functions define complex-valued GPs whose real part and imaginary part are in quadrature relation in expectation. From Eq. 24 and Eq. 29, this implies that the argument of the coherency has to be equal to −i sgn(ξ). By direct computation, it is easy to see that the cross-spectrum (the numerator of the coherency) can be expressed in terms of eigenvalues as follows (see appendix II): where the complex numbers a and b solely depend on the eigenvectors v (1) and v (2) . From this expression, together with the fact that the denominator of the coherency is a positive real number, it follows that arg[γ(ξ)] does not change if both the eigenvalues are multiplied by a positive constant while the eigenvectors are left unchanged. More generally, arg[γ(ξ)] is left unchanged when both the positive-valued eigenvalue functions λ (1) (ξ) and λ (2) (ξ) are convolved with the same positive-valued function. Therefore, equations 36 and 37 imply that arg[γ(ξ)] is the same for quadrature and quasi-quadrature covariance functions and, consequently, the quadrature relation holds in expectation. Furthermore, as the Fourier transform of a smooth rapidly decreasing function such as f (τ ) is itself rapidly decreasing, the first (second) eigenvalue will always dominate when ξ is positive (negative) and has a large magnitude. Therefore, the Fourier coefficients of a stochastic process with quasi-quadrature covariance function will be in almost exact quadrature relation for high frequencies.
In the remaining, we will use a squared exponential envelope function: where δ is its characteristic time scale. Fig. 3C shows the quasi-quadrature covariance function obtained by multiplying a quadrature periodic covariance function with a squared exponential envelope (δ = 0.2s). The envelope breaks the exact periodicity of the function, thereby allowing to model signals that are not exactly periodic. Fig. 3D shows the spectrum of this covariance function. As effect of the envelope function, the negative frequencies are not exactly zero.

Results
We begin this section by comparing the performance of CGPR and quadrature filter in estimating instantaneous amplitude and phase of different kinds of stochastic oscillations. These simulations have been performed on noise free data as the quadrature filter is easily destabilized by broadband noise. Finally, in the third subsection, we report on the analysis of real magneto-encephalogram (MEG). In this analysis, we use the CGPR for extracting instantaneous amplitude and phase of human alpha oscillations.

Comparison between the analytic representation and the GP complex signal
In this subsection we compare the performance of a CGPR with the more conventional analytic representation for the purpose of quantifying instantaneous amplitude and frequency of an oscillatory signal. We focus on two particular classes of non-stationary signals: deterministic chirplets and stochastic oscillations. We make use of a periodic quadrature covariance with a squared exponential envelope. The same covariance function was used in the two simulations.

Analysis of deterministic chirplets
A chirplet is a signal whose instantaneous frequency increases as a function of time and whose amplitude is concentrated in a small time window (Yang et al., 2013). The analysis of chirplets is relevant in many scientific fields. For example, gravitational waves have a chirplet waveform (Adams et al., 2013). In this simulation study, we aim to asses the performance of the methods in recovering instantaneous amplitude and frequency of deterministic chirplets. We generate the chirplets by multiplying a cosine wave (with either linear, quadratic or exponential frequency increase) with an amplitude envelope.
where k is an integer from 1 to 3, b is a parameters that regulates the width of the envelope and h(t) is the Heaviside step function. The instantaneous phase Φ(t) is given by (ω 0 +at)t (linear frequency growth), (ω 0 +at 2 )t (quadratic frequency growth) or (ω 0 +ae t/2 )t (exponential frequency growth) plus a random initial phase. In total, we have nine possible combinations of envelope and phase functions, each having the width b, the initial frequency ω 0 and the frequency factor a as parameters. We run a simulation study by randomly selecting 500 of these combinations, with the parameters sampled from the uniform distributions b ∈ [0.1, 0.3], 2πω 0 ∈ [0.1, 2], k ∈ {1, 2, 3} and 2πa ∈ [0.1, 0.4]. For each of these chirplets, we extracted instantaneous amplitude and frequency using both the analytic representation and the GP-based complex signal. The latter was obtained using a periodic quasi-quadrature covariance with squared exponential envelope and parameters ω 0 = 0.5Hz, ρ = 0.1 and δ = 0.3s. In this analysis, the signal did not contain any noise, and therefore we set the noise covariance matrix in Eq. 23 to zero. Fig. 4 shows the results for a representative chirplet. Panel A shows the real and imaginary parts of the GP-based complex signal together with the amplitude computed as its complex modulus. Note that, since the GP regression was noise free, the real part is identical to the original signal. Panel B and C compare, respectively, the amplitude and frequency obtained through the complex GP with the "ground truth" given by and P(t) = Φ(t) (see Eq. 40).
Panels D, E and F show the same results for the analytic representation. From this simulation, it is clear that the analysis with the analytic representation fails to estimate the frequency in the low-amplitude/high-frequency segment and this failure is accompanied by a ringing of the estimated amplitude envelope. In addition, the analytic method exhibits serious boundary effects even though the signal is very close to zero at the boundaries. These effects arise from the fact that the quadrature filter is performed in the frequency domain and, consequently, requires a prior discrete Fourier transform of the signal. Conversely, the GP correctly tracks both amplitude and frequency and does not show any boundary effect. Fig. 5 shows the outcome of the comparison between the two methods. The errors of both methods where quantified as mean absolute deviation from the ground truth. Panels A and B show the comparison between GP-based complex signal and analytic representation for, respectively, amplitude and frequency. Every blue circle corresponds to an individual chirplet and the red dots indicate the mean error of the two methods. The GP method clearly outperforms the analytic representation for both amplitude and frequency estimation. The mean error for the GP-based complex signal is 0.02 and 0.31 for, respectively, amplitude and frequency, as compered to 0.04 and 1.24 for the analytic representation. Panels C and D show the performance of the methods as a function of the initial frequency. The difference between the GP-based complex signal and the analytic representation is maximal for very low frequencies and decreases with frequency.

Analysis of stochastic oscillatory processes
Now, we compare the methods with respect to the ability to reconstruct amplitude and frequency of a stochastic oscillatory signal. We generates stochastic oscillations by multiplying a smooth random amplitude function with a sinusoidal wave with a smooth random frequency function. Specifically, the oscillatory signal s(t) was generated as the following Figure 4: Analysis of a chirplet. A) Real part (blue), imaginary part (green) and amplitude (red) of the GP-based complex signal, since the GP regression was noise free, the real part is identical to the original signal. B,C) Comparison between the instantaneous amplitude and frequency obtained through the complex GP (green) and the ground truth (blue). D) Real part (blue), imaginary part (green) and amplitude (red) of the analytic representation. E,F) Comparison between the instantaneous amplitude and frequency obtained through the analytic representation (green) and the ground truth (blue). non-linear function of two GPs: In this formula, the amplitude process A(t) is given by a positive-valued non-linear transform of a GP: in which x(t) is a GP with a mean identically equal to √ 3 (so that the mean of A(t) is 1) and a covariance that is a squared exponential with time-scale δ x : 2δ 2 x .
The instantaneous phase function P(t) is a GP with mean ω 0 t (the phase of a stationary oscillation with frequency ω 0 ) and squared exponential covariance function with time-scale δ P : Finally, the initial phase φ 0 is a uniform random variable ranging from −π to π. Furthermore, the samples of the oscillatory process were multiplied by a Hann taper in order to reduce the boundary artifacts of the analytic representation. We ran a simulation study by generating 500 stochastic oscillatory signals. Each signal was generated by randomly initializing the parameters δ x , ω 0 and δ P , sampling x(t) and P(t) from the respective multivariate Gaussian distributions and combining them using Eq. 41 and 42. The ranges of the parameters were δ x ∈ [0.1, 0.4], 2πω 0 ∈ [1, 3] and δ P ∈ [0.1, 0.4]. Fig. 6 shows the results of the simulation study for amplitude (panel A) and frequency (panel B) estimation. The complex GP signal outperforms the analytic representation in most cases, although the difference is less pronounced than for the simulation study usingchirplets. T he mean error of the GP is 0.03 and 0.30 for, respectively, amplitude and frequency, which is to be comared with mean errors of 0.04 and 0.41 for the analytic representation.

Experimental results on brain data
In this subsection we apply the CGPR to real MEG data, and use it to study the nonstationary dynamics of brain waves. We focus on alpha oscillations (8)(9)(10)(11)(12)(13). In this analysis we use the harmonic quasi-quadrature covariance function (Eq. 34) since MEG recordings have a low SNR and alpha oscillations have an anharmonic waveform (Roberts et al., 2015).

Data acquisition and preprocessing
We collected resting state brain activity from an experimental participant that was instructed to fixate a cross at the center of a black screen. The study was conducted in 913 accordance with the Declaration of Helsinki and approved by the local ethics committee 914 (CMO Regio Arnhem-Nijmegen). Informed written consent was obtained from all 915 participants. Brain activity was recorded using a 275 axial gradiometer MEG setup (VSM/CTF Systems, Port Coquitlam, British Columbia, Canada). The acquisition sampling rate was 1200 Hz. Preprocessing was performed using the open source Matlab toolbox FieldTrip (Oostenveld et al., 2010). The continuous data was cut in segments of two seconds. Trial segments containing movement, muscle, and SQUID artifacts were discarded following a semi-automatic artifact detection routine. A fourth-order Butterworth band-stop filter (49-51 Hz) was used to remove 50 Hz line noise. Finally, residual eye blinks, heartbeat, line noise and muscular artifacts were isolated by independent component analysis (Comon, 1994) on the concatenated segments and subsequently removed from the data.

Analysis of alpha oscillations
We restricted our attention to the analysis of the MEG sensor with the greatest alpha power. This sensor was located in the occipital part of the helmet. We used the CGPR for constructing the real and the imaginary imaginary part of the oscillation and computing the instantaneous amplitude and frequency. MEG measurements are corrupted by strong low frequency noise, often described as pink noise (Chen et al., 1998). We modeled this component using the following noise covariance matrix: We used this covariance matrix in Eq. 23. The parameter σ(= 0.5) is the standard deviation of the white noise while ν(= 1) and b(= 1) are respectively standard deviation and relaxation parameter (inverse of time constant) of the low frequency broadband noise, here modeled using an exponential covariance function. The parameters of the noise and the oscillation were chosen in a way to fall within the range of real human MEG signals (see (Ambrogioni et al., 2016), where we used a least-squares method to determine the optimal parameter values for the oscillatory and the noise covariance functions.). Fig. 7A shows an example MEG signal. The alpha oscillation is clearly visible in the initial and the final segment of the trial, although its amplitude visibly decreases in the middle. These fluctuations in the alpha amplitude are well known in the neuroscience community (Klimesch et al., 2011). Fig. 7B shows the real and the imaginary part of the GP-based complex signal obtained from this signal, together with its instantaneous amplitude. We see that the combination of harmonic quadrature covariance function and noise model given by Eq. 43 allows to effectively filter out the white and low frequency broadband noise and to obtain a smooth estimate of the amplitude envelope. Fig. 7C shows the estimated instantaneous frequency. The frequency tends to fluctuate around 10 Hz except when the amplitude is very low. In the latter situations, the frequency shifts by several Hertz.
We further investigated the relation between instantaneous amplitude and frequency. This analysis could potentially constrain the possible network mechanisms that could generate the alpha oscillation. Fig 7D shows the density plot of instantaneous amplitude and frequency across all trials and time points. The mean frequency does not seem to systematically change as a function of the amplitude. Conversely, the frequency fluctuations clearly decreases when the instantaneous amplitude increases. The latter effect is probably due to the relation between instantaneous amplitude and signal-to-noise ratio (SNR).

Discussion
In this paper, we introduced a new method for estimating instantaneous amplitude, phase and frequency of narrow-band real-valued signals using a complex-valued version of GP regression. The main innovation of the paper is the definition of two new classes of complexvalued covariance functions (i.e., quadrature and quasi-quadrature covariance functions) that induce either an exact or an approximate quadrature relation between the real and the imaginary part of random signals. Using these covariance functions, we can construct the unobservable imaginary part from the observable real part of the signal. A quadrature covariance function is obtained by applying the quadrature filter on a stationary real-valued covariance function and it defines GPs whose samples are analytic functions with probability one. We argued that, when the spectra of amplitude envelope and oscillation are highly overlapping, the analytic complex valued signals result in rather counter-intuitive estimates of instantaneous amplitude, phase and frequency. To obtain more plausible estimates for this type of signals, we introduced quasi-quadrature covariance functions. A quasi-quadrature covariance function is obtained by multiplying a quadrature covariance function with a realvalued envelope. We showed that these covariance functions specify GPs whose samples are not analytic while maintaining the quadrature relations in expectation between the real and the imaginary part. By analyzing simulated chirplets and wavelets, we showed that the CGPR equipped with a quasi-quadrature covariance function outperforms the quadrature filter in recovering instantaneous amplitude and frequency. In particular, we found that the improvement is higher when the frequency of the signal was smaller (see Fig. 5 C and D). This result was expected as we showed that the deviation from the exact quadrature relation between real and imaginary part is largest for very low frequencies as for these frequencies the spectra of the amplitude envelope and the oscillation overlap most.

Practical considerations
We now discuss some practical considerations concerning the application of CGPR analysis on real data. An important feature of CGPR (and GP regression in general) is its ability to exploit prior knowledge about signal and noise processes through the choice of appropriate signal and noise covariance functions. For example, in our analysis of MEG data we modeled the signal as a complex-valued harmonic oscillation corrupted by two sources of real-valued noise: a white noise and an Ornstein-Uhlenbeck process. Importantly, the parameters of the different covariance functions (both for the noise and the oscillatory processes) can be inferred from the signal itself. For example, the peak frequency parameter of the harmonic or the periodic covariance function can be optimized in order to increase the performance of the method. The optimization can be done by maximizing the marginal likelihood of the model (Rasmussen, 2006) or by a non-linear least-squares approach (Ambrogioni et al., 2016). Since we assumed the imaginary part to not be directly observable, these optimization procedures are completely equivalent to the ones for the associated real-valued GP regression.
One of the shortcomings of the CGPR is its cubic time complexity. In fact, the computation of the posterior expectation (Eq. 23) requires the inversion of the prior covariance matrix between all sample points. The cubic complexity is problematic for the analysis of long signals. Fortunately, there are several methods that can effectively deal with this problem. A possible solution is to use an envelope function with a compact support, i.e. one that is different from zero only for a short time lag. This kind of covariance functions allows to work with more computationally advantageous sparse covariance matrices since many of the entries will be exacly zero (Wendland, 2004). An alternative approach is to convert the CGPR into a complex-valued infinite dimensional Kalman filter that can be solved in linear time (Sarkka et al., 2013). This conversion is exact for a large class of stationary covariance functions and approximate but arbitrarily accurate for the others. Another, possibly simpler, approach is to perform the CGPR in the frequency domain, this gives linear time complexity but introduces periodic boundary conditions as it requires a prior Fourier transform of the data (Paciorek, 2007).

Extensions of the CGPR method
The CGPR analysis can be adapted to non-Gaussian data by adopting different likelihood models. For example, the analysis can be applied to binary data by introducing a latent complex-valued Gaussian process whose real part determines the probability of a Bernoulli variable (Rasmussen, 2006). This approach can be used to estimate the instantaneous phase and frequency of a binary signal. Adopting a different likelihood, such as the studentt instead of the Gaussian, can also make the analysis more resistant to outliers (Jylänki et al., 2011).

Conclusion
In sum, we have described a new method for estimating instantaneous amplitude, phase and frquence of narrow-band real-valued signals using a complex-valued version of GP regression (CGPR). We showed that this CGPR equipped with quasi-quadrature covariance functions provides a much better estimate of the instantaneous amplitude and frequency than the established approach using the quadrature filter. CGPR is a versatile tool because it allows to incorporate prior information about the structure in signal and noise and thereby to tailor the analysis to the features of the signal.
Appendix I: Hermitianity and non-negative definiteness of quadrature and quasi-covariance functions In this appendix we show that quadrature and quasi-quadrature covariance functions are both positive definite and Hermitian. A kernel function K(τ ) is said to be Hermitian when: This implies that the real part of K(τ ) is a even function while its imaginary part is an odd function. The quadrature covariance function Ak(τ ) is obtained by applying the quadrature filter A on a stationary covariance function k(τ ). On one hand, the real part is Ak(τ ) is equal to k(τ ) and is therefore symmetric (otherwise k(τ ) would not be a valid covariance function). On the other hand, the imaginary part Hk(τ ) is odd since where we defined the new integration variable r = −s and leveraged the fact that k(τ ) is even. The Hermitianity implies that the Fourier transform is real-valued. This follow from In fact, the first Fourier transform on the right hand side is real-valued since k(τ ) is even while the second Fourier transform is purely imaginary since k(τ ) is odd. An Hermitian kernel function K(τ ) is said to be non-negative definite when its Fourier transform is non-negative (almost) everywhere. The Fourier transform Ak(τ ) is 2h(ξ)k(ξ) that is obviously non-negative if k(τ ) is a valid covariance function.
Quasi-quadrature covariance functions are obtained by multiplying a quadrature covariance function with a real-valued even positive definite envelope function. Clearly the pointwise product of an Hermitian function with a real-valued even function is itself Hermitian. Furthermore, a pointwise product of two non-negative definite functions is always non-negative definite because the convolution of two non-negative functions is always non-negative. Therefore, quasi-quadrature covariance functions are indeed valid covariance functions.