Bayesian filtering framework for noise characterization of frequency combs

Amplitude and phase noise correlation matrices are of fundamental importance for studying noise properties of frequency combs. They include information about the origin of noise sources as well as the scaling and correlation of the noise across the comb lines. These matrices provide an insight that is essential for obtaining low-noise performance which is important for, e.g., applications in optical communication, low–noise microwave signal generation, and distance measurements. Estimation of amplitude and phase noise correlation matrices requires highly–accurate measurement technique which can distinguishes between noise sources coming from the frequency comb and the measurement system itself. Bayesian filtering provides a theoretically optimum approach for filtering of measurement noise and thereby, the most accurate measurement of phase and amplitude noise. In this paper, a novel Bayesian filtering based framework for joint estimation of amplitude and phase noise of multiple frequency comb lines is proposed, and demonstrated for phase noise characterization. Compared to the conventional approaches, that do not employ any measurement noise filtering, the proposed approach provides significantly more accurate measurements of correlation matrices, operates over a wide range of signal–to–noise–ratios and gives an insight into comb’s dynamics at short scales (<10−8 s). © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Measurement of amplitude and phase noise of optical frequency combs can be performed using a variety of analog and digital measurements techniques, see e.g. [1,2]. Typically, analog techniques require complex calibration procedures and phase-locked loops [3]. Most importantly, as the functionalities are implemented in the analog domain, it becomes very challenging to implement any sophisticated signal processing techniques, especially, at high-speeds (beyond 10 GHz). Digital measurement techniques on the other hand, employ a balanced receiver, an analogue-to-digital converter (ADC) for sampling the signal of interest and a computer for storing and processing the data offline. The measurement of amplitude and phase noise is thereby performed offline, on the sampled data, allowing for the implementation of advanced digital signal processing and machine learning techniques. What makes digital measurement techniques even more attractive is that recent advances in optical communication systems have enabled balanced receivers operating with up to 100 GHz of analog electrical bandwidth and ADCs with sampling rates in the range of 160 Gs/s [4]. This is opening up for new opportunities for ultra-broadband noise characterization of frequency combs.
Irrespective of the measurement technique, an effective method for filtering the measurement noise must be employed. This is because the thermal noise floor will set a limit on the magnitude and the frequency range in which the amplitude and phase noise can be measured. Most of the conventional techniques do not provide effective filtering of measurement noise and are therefore limited to high signal-to-noise ratio (SNR) of frequency comb lines, MHz measurement ranges and approximately −140 dB rad 2 /Hz phase noise levels [1,2] [5][6][7]. Most importantly, for correlation matrix computation joint amplitude and phase noise estimation of multiple frequency lines, which may easily have different SNRs, is needed.
In many practical cases, the requirement on high SNR per frequency comb line is difficult to satisfy. This is because the output power of the frequency comb is divided among multiple carriers, resulting in lower available power, and hence low SNR, per comb line. This makes accurate noise characterization using conventional approaches rather challenging. A common approach is to employ optical or electrical amplification prior to the phase noise measurement, to increase the signal power. However, this approach also increases the thermal noise level due to amplifier noise. It is therefore of great significance to have an amplitude and phase noise measurement technique that is highly-sensitive and works well across wide range of SNRs.
It has recently been demonstrated that Bayesian filtering provides an optimum filtering of the measurement noise allowing for a record sensitive optical phase noise measurement. In Ref. [8], phase noise measurements from ultra-low power signals (−73 dBm), as well sensitivities down to −200 dB rad 2 /Hz and measurement ranges up to 10 GHz offset frequencies have been demonstrated for a single frequency laser. In short, the results in [8] demonstrated that Bayesian filtering can be used to surpass the thermal noise limit of the detection system. Bayesian filtering thus holds a great potential for providing a unifying framework for noise characterization of lasers and optical frequency combs overcoming the limitations of measurement noise.
In this paper, a Bayesian filtering framework for joint estimation of amplitude and phase noise of multiple frequency comb lines is proposed. The theoretical foundation and solution to the Bayesian filtering equations are presented. The extracted amplitude and phase noise traces are used to compute the corresponding amplitude and phase noise covariance and correlation matrices. The subsequent eigenvalue decomposition is also presented and it provides an insight into how various noise sources (from e.g. the input continuous wave laser and radio frequency oscillators) contribute to the amplitude and phase noise across difference comb lines. Simultaneous detection of multiple lines is realized using a digital measurement technique in which the comb under test (CUT) is heterodyned with another local oscillator (LO) frequency comb, similar to a dual-comb interferometer [9,10]. In this manner, the noise performance of multiple lines is captured at once [7], allowing for computation of the amplitude and phase noise correlation and covariance matrices with single line resolution. The proposed framework in this paper, is much more general and complete compared to [8] in which only phase noise of a single frequency laser is estimated. Moreover, compared to the proof-of-principle results we presented in [11], this paper includes significant extensions in terms of the method itself, as well as numerical and experimental results.
We numerically and experimentally validate the performance of the framework with electrooptic frequency comb sources reported in [9], focusing only on phase noise characterization. These frequency combs have well-defined phase noise characteristics that are dictated by the performance of the radio-frequency oscillator driving the modulators and the optical phase noise of the input continuous-wave laser [12], thus providing a reliable source to assess the framework. The performance in terms of accuracy is investigated numerically by computing the mean square error (MSE) as a function of SNR. Experimentally, we demonstrate accurate joint phase noise tracking and estimation of 49 frequency comb lines, simultaneously. The estimated phase noise traces are later used for computation of the phase noise correlation matrix and compared to the state-of-the-art method that does not employ any measurement noise filtering. Finally, eigenvalue decomposition of the phase correlation matrix is performed to get an insight into magnitude of noise sources and their evolutions.
The rest of the paper is structured as follows. In section 2, the state-space model for joint estimation of amplitude and phase noise of multiple frequency comb lines is presented. The state-space model provides the necessary framework for the solution of the Bayesian filtering equations. In section 3, numerical simulations are presented. The benefits of the proposed framework, in terms of the mean square error, eigenvectors and eigenvalues estimation, are investigated and compared to the state-of-the-art method. In section 4, experimental results for the phase noise correlation matrix estimation as well as the corresponding eigenvalue decomposition are presented for an electro-optic comb.

State-space model
In Fig. 1, a dual-comb heterodyne set-up for frequency comb noise characterization, using a digital measurement technique, is shown. It employs a comb under test (CUT) that is mixed with a local oscillator (LO) comb in a balanced receiver. The resulting signal, y(t), is a down-converted comb, characterized by M spectral lines. The spacing ∆f is equal to the difference in repetition rates between the CUT and the LO comb. The down-converted signal y(t) is then sampled by an ADC and the product M∆f is kept within the bandwidth set by the minimum between the Nyquist condition for sampling and the 3 dB bandwidth of the balanced receiver BW rec 3dB , i.e. M∆f ≤ min(F s /2, BW rec 3dB ), where F s is the sampling frequency of the ADC. The sampled photocurrent y k , where k is a discrete-time index, is then stored and passed to the Bayesian filtering framework for off-line signal processing. By heterodyning the two combs, the m-th line phase, of the downconverted comb, φ m k is given by the phase difference: In this paper, we assume that the CUT and the LO comb are generated by two identical and independent laser and RF sources, hence they equally contribute to the downconverted comb amplitude and phase noise. Given the stored data, y 1:T := [y 1 , . . . , y T ], where T is the length of the measurement data in samples, the objective is to perform joint tracking of the downconverted frequency combs' phases φ k := [φ 1 k , . . . , φ M k ] and amplitudes a k := [a 1 k , . . . , a M k ] for k = 1 · · · T, where denotes the transpose operator. We consider amplitude and phase noise of the donwconverted comb as a discrete-time stochastic processes, indexed by k and m in time and line number, respectively. Implementation of the Bayesian filtering framework requires a state-space model (SSM) that consists of: 1) a measurement equation, which describes the relation between measurements y k and the time-varying phase and amplitude noise, i.e. [φ k , a k ] and a state equation 2) which is an approximate model describing the evolution of the optical phase and amplitude noise. In this paper, a random walk model is employed. For any given time index k, the proposed SSM is: where initial amplitudes and phases are set to zero, i.e. φ m 1 = 0 and a m 1 = 0. Here F s = 1/T s denotes the sampling frequency of the ADC and n k is the measurement noise. It is assumed that n k has a Gaussian distribution with zero mean and variance σ 2 n . It contains contributions from thermal and quantum noise sources associated with photodetection, amplifiers and ADC. The angular frequency difference between the comb lines is expressed as ∆ω m = 2π∆fm.
are multivariate Gaussian distributed random variables with zero mean and covariance matrices Q A and Q φ , respectively. The diagonal elements of Q A and Q φ describe our assumptions on the variance (strength) of amplitude and phase noise sources, respectively. The off-diagonal elements describe the phase and amplitude correlation between M frequency combs lines, respectively. The cross-covariance matrix between amplitude and phase noise of comb lines is denoted Q A,φ . All these matrices together, form the complete noise matrix Q, which is defined as: The matrix Q includes all the necessary information to get an insight into noise sources and correlations of the frequency combs under consideration.

Bayesian filtering
As we only have access to the measurements y k , the objective of Bayesian filtering is to estimate amplitude and phase noise which are explicitly hidden in y k . (A detailed treatment of Bayesian filtering is beyond the scope of this article and the interested reader is referred to [13]). We assume that the underlying model of a m k and φ m k is governed by (1) and (2), and the measurements depend on a m k and φ m k via (3). Strictly speaking, given y k , we would like to estimate amplitudes and phases that minimize the following MSE problem: denotes the statistical expectation and x true := [φ true k ; a true k ] are true amplitude and phase noise profiles. The solutions are then said to be optimal in the mean square error sense and are theoretically, the most accurate estimations of amplitude and phase noise in the presence on measurement noise. We denote these solutions as x According to the Bayesian filtering theory, the optimal estimate x opt k are obtained given by computing the expectation of the filtering distribution [13]: In (5) p(x k |y 1:k ) is the posterior probability distribution, or filtering distribution of x k given the measurement from 1 to k, i.e. y 1:k . Typically, all solutions to (5) are iterative in k. An exact solution to (5) can be obtained if the measurement Eq. (3) is linear and is expressed via Kalman filtering equations. Since the measurement Eq. (3) in our case is nonlinear in x k , only an approximate solution to (5) can be obtained. There are several approximate solutions, however, and in this paper we consider a solution based on an extended Kalman filter (EKF) due to its ability to handle large traces of data (we consider the size y k of up to 60 · 10 6 ). Taking into considerations the state-space model specified in (1)-(3), the EKF filtering equations for every time step k = 1 · · · T are expressed as following: Implementing the system of equations in (6) requires initial conditions for x 0 . Typically, is the gradient of the measurement function (3) with respect to the hidden states x k . h x (·) contains the gradient with respect to the phases h φ (·) and the amplitudes h A (·) expressed as: The EKF Eqs. (6), once executed for every time-step k, return the optimal estimations of the amplitude and phase noise x This is under assumption that we have the knowledge of the correlation matrix Q, variance of the measurement noise σ 2 n and the frequency offset ∆f . These are so called static parameters and their estimation needs to be performed in conjunction with the EKF equations as illustrated in Fig. 1 and explained next.

Static parameter estimation
The frequency spacing ∆f and the measurement noise variance σ 2 n can be extracted by inspection of the power spectral density of the photocurrent prior to the Bayesian filtering. The matrix Q, however, needs to be jointly estimated from the observed data y 1:T together with x 1:T . The Bayesian treatment for the optimal estimation of Q, in the MSE sense, is obtained by maximizing the posterior distribution of Q given the observations y 1:T : where p(y 1:T |Q) is the likelihood function. It describes how likely it is that a particular choice of matrix Q is responsible for generating the observations y 1:T . The prior probability describes our belief (initial values) on Q prior to any observation y k . Domain knowledge can typically be very useful in determining good priors which can significantly reduce the convergence time of the static parameter search algorithm.
To avoid numerical approximation errors due to limited precision, it is common practice to work in the log domain of the probabilities, transforming then the maximization of the likelihood in (8) into the minimization of the negative-log-likelihood LL − (Q) := − log p(y 1:T |Q). Then, LL − (Q) is the cost function we want to minimize. Typically, for a general state'space model LL − (Q) can be approximated as [13]: in which v k and s k are computed by the EKF Eqs. (6) for a known matrix Q. Equation (8) is then reformulated as: An efficient approach to minimize (9) is to employ the expectation maximization (EM) algorithm [13], see Fig. 1. The advantage of EM over other optimization techniques is that at each iteration the optimal parameter, in our case the covariance matrix Q opt that minimizes (9), is returned in a closed form, and the procedure is theoretically guaranteed to converge [14]. EM involves a two step procedure, starting from an initial guess of the matrix Q (0) : the expectation step (E) followed by the maximization step (M), which are then iterated N times. During the n-th EM iteration the E step requires computation of the smoothing state estimate x s k := [φ s k ; a s k ] . This smoothing estimate can be computed once the sequence has been filtered with the EKF using the Extended Kalman Smoother (EKS). The EKS runs backwardly in time and produces a smooth state estimate sequence x s 1:K with associated smoothed covariance sequence Σ s 1:K . Here K represent the number of data samples used for parameter estimation, and typically N T. The EKS is initialized by setting Then, for every k = K − 1 · · · 0 the following equations are computed where Q (n) is the estimation of the matrix Q during the n-th EM iteration. Once completed the filtering and smoothing procedure, the M step concludes by calculating the updated parameter Q (n+1) using where P, C and B are defined as EM iterates by repeating the E and M steps, invoking EKF and updating the values of Q at the end of each M step. EM can be halted when the improvement on the cost function becomes negligible, that is when a minimum of (9) is found. Typically, the EM does not need that many iterations to converge. Our experience is that the EM will converge before N = 1000 iterations, and Q reaches its optimal value, i.e. Q (N) ≈ Q opt . Once the EM has converged, Q (N) is used in the EKF equations to estimate amplitude, a opt 1:T , and phase noise, φ opt 1:T , respectively, see Fig. 1(c).

Covariance and correlation matrix computations
Noise characterization of the frequency comb can the be performed by computing phase noise and RIN spectrum using φ opt 1:T and a opt 1:T , respectively. Additionally, the sample amplitude and phase noise covariance matrix can then be calculated on a window of samples of length K, embracing the contiguous temporal indices from l to n: where u = a opt 1:T or u = φ opt 1:T . K = n − l + 1,ū l:n denotes the mean estimate of a k or φ k over the current window, i.e.φ l:n = K −1 n k=l φ k andā l:n = K −1 n k=l a k . The reason for choosing the window length is because the dynamics of the frequency comb, and thereby the matrix coefficients, may be time-dependent. The correlation matrix is directly obtainable from the covariance matrix (16) by performing normalization of the diagonal: ρ u,l:n = diag Σ u,l:n −1/2 Σ u,l:n diag Σ u,l:n where the diag(·) operator sets the anti-diagonal elements of a matrix to 0. The diagonal terms of Σ u,l:n contain the information about the variance of amplitude or phase noise associated with the different frequency lines. The off-diagonal coefficients of matrix, ρ u,l:n , contain the information about the linear correlation among amplitudes and phases, ranging from 1 (perfect correlation) to -1 (perfect anticorrelation).

Numerical results
In this section, we investigate the performance of the proposed framework to estimate the phase noise covariance and correlation matrix of an EO comb. The dominant type of noise, for these combs, is the phase noises originating from the continuous wave (CW) laser and the RF oscillator.
In the simulation set-up, we generate directly the down-converted comb signal, y k that contains frequency components at ∆ω m = 2π∆f RF + 2π(m − M/2 + 1)∆f , where ∆f RF = 2.5 GHz and ∆f = 100 MHz. We consider M = 49 comb lines. The CW laser and the RF oscillator phase noise are simulated as Wiener processes. The combined linewidths for the laser and the RF oscillator are assumed to be ∆ν L = 1 kHz and ∆ν RF =0.5 kHz, respectively. For an EO comb, the phase noise of the down-converted frequency comb line m is assumed to be line-dependent and generated as [12]: where m r = m − M/2 + 1 is the relative line index, · is the floor function, φ L k and φ RF k are phase noise sources associated with CW laser and RF oscillator phase noise. The sampling frequency is F s = 10 GHz. Measurement noise is added to signal y k , resulting in the following, In Fig. 2, the ability of the proposed framework to learn the system parameters is shown. We start with our guess (prior) on Q φ at iteration k = 0, i.e. Q (0) φ . For the considered case, we assume that Q (0) φ is a diagonal matrix. The true matrix denoted by Q true φ is not diagonal as shown in Fig. 2. In fact, since the EO comb phases are generated according to (18), the true matrix can be expressed as Q true φ = 2πT s [∆ν L + m r m r ∆ν RF ], where m r = [− M/2 , . . . , M/2 ] . Now, as the number of EM iteration is increased, the estimated Q φ approaches the true matrix Q true φ . This is also indicated by the evolution of LL − (Q) which decreases as the number of EM iterations is increased. Even though our initial assumption on the initial guess Q (0) φ was wrong, the framework was still able to converge to the correct Q φ . Next, we investigate the effectiveness of the framework to perform filtering of the measurement noise and operate at a wide range of SNRs. We also perform a benchmarking with a conventional digital phase noise measurement technique, shown in Fig. 3, that does not employ any filtering of the measurement noise [15]. The approach shown in Fig. 3 employs a bank of bandpass (BP) digital filters, with bandwidth f BW , to extract different frequency combs lines. Thereafter, for each filtered line, denoted by i m k , a Hilbert transform is applied to recover the corresponding orthogonal quadrature q m k . The phase information, per comb line, is estimated by taking the inverse trigonometric tangent of reconstructed field, i.e. tan −1 (i m k /q m k ). The resulting phase is processed through an unwrapping function that smooths every phase jump greater than π by adding multiples of 2π to obtain smoother transitions and to make the phase as a continuous In Fig. 4, we compare the accuracy of the conventional and the Bayesian filtering based approach for estimation of the correlation matrices when varying SNR avg . The true correlation matrix is denoted by ρ true φ,T and is depicted in left column in Fig. 4(a). The structure of ρ true φ,T is as expected for an EO comb. For the particular case, of the simulated EO comb, it can be seen that the lines closer to each other have a positive phase correlation, while the distant lines are anti-correlated. The estimated matrices obtained by the conventional and the Bayesian method are denoted by ρ conv φ,T and ρ BF φ,T , respectively. In Fig. 4(a) these are shown for the SNR avg of 10, 15 and 20 dB. The noise measurement bandwidth is 100 MHz. Figure 4(a) illustrates qualitatively that the estimated correlation matrices, ρ BF φ,T , are very similar to ρ true φ,T for the considered ranges of SNR avg . This is not the case for ρ conv φ,T . The difference between ρ conv φ,T and ρ true φ,T is especially pronounced for SNR avg of 15 and 10 dB. The reason for the discrepancy is the inability of the conventional method to provide accurate phase noise estimation for multiple comb lines. This is investigated in more details in Fig. 4(b) where the Root Mean Square Error (RMSE), between the estimated and the true phase noise traces computed for all M = 49 lines, is plotted as a function of SNR avg . The Bayesian approach has a significantly lower RMSE compared to the conventional one. This is due to the intrinsic property of the Bayesian filtering to filter out measurement noise when performing phase noise estimation, a property that is not part of the conventional phase estimation method. In addition, we also investigated how the Laser source linewidth ∆ν L impacts the system performance. If such linewidth is increased up to 1 MHz, it is possible to notice in Fig. 4(b) an increase in the RMSE. This shows that larger phase displacements cause performance loss in the phase tracking framework.
To investigate the benefits of the proposed framework in greater details, we compute the variance of the differential phase ∆φ m k = φ m k − φ c k , where c is the line index of the central line. For an EO comb, the variance of ∆φ m k is expected to be a quadratic function of the line index due to line-dependent phase noise generation process, namely σ 2 [16]. We consider a relatively high SNR avg of 25 dB, and calculate the variance on the whole signal duration T. Figure 5 shows that the variance of the differential phase obtained using Bayesian filtering practically overlaps with the true one. The differential phase variance obtained using the conventional approach shows an offset and also an erratic behaviour due to the effects of the measurement noise. Lowering the SNR avg , results in an even more inaccurate computation of the differential phase variance using the conventional method. An insight into the of sources of phase noise and their contributions to individual frequency comb lines can be obtained by performing the eigenvalue decomposition of the phase noise correlation matrix. More precisely, the eigenvectors will depict how the sources of phase noise are contributing to the phase noise of each frequency comb lines. The eigenvalues will be directly proportional to the variance of the phase noise sources. For the considered case, the CW laser and the RF oscillator are the main sources of phase noise and therefore only two eigenvalues (λ 1 , λ 2 ) and eigenvectors (v 1 , v 2 ) will be significant. The other M − 2 eigenvalues and the correpsonding eigenvectors will be close to zero.
The question is now to which phase noise source, the CW laser or the RF oscillator, (λ 1 , v 1 ) and (λ 2 , v 2 ) belong to. It can be shown that, using (18), the contribution from the CW laser will be equal for all frequency comb lines, while the contribution from RF oscillators will scale linearly with m. This is exactly what eigenvectors (v 1 ) and (v 2 ) are showing in Fig. 6. We can therefore conclude that the eigenvalue (λ 1 ) and eigenvector (v 1 ) are associated with the RF oscillator. Figure 7 shows the evolution of two main eigenvalues as a function of the observation time. We also plot the evolution of CW laser and RF oscillator phase noise variances. It is observed that the eigenvalues converge to the phase noise variance of the CW laser and RF oscillator with a proportionality constant equal to the norm of vector v 1 and v 2 . Since the phase noises of  For shorter observation times (τ obs <10 −8 s), there is a discrepancy between (λ 1 , λ 2 ) and (σ 2 RF v 1 2 , σ 2 L v 2 2 ). However, the discrepancy is more significant for (λ 1 , λ 2 ) obtained using the conventional phase noise estimation. Even for the Bayesian filtering approach, a relatively small discrepancy exists. The reason is that it is quite challenging to effectively filter the measurement noise at such short time-scales. The consequence of this is that there is uncertainty in phase estimation at shorter time scales and this uncertainty is significantly increased if the measurement noise is not filtered.

Experimental results
For the experimental validations, a set-up shown in Fig. 1(a) is used [15]. Two EO combs with line-spacings of 25 GHz and 25.05 GHz, respectively, are employed. It is assumed that the CUT and LO comb have equal and independent contributions to the overall phase noise and, therefore, the estimated phase noise covariance needs to be divided by a factor 2 [17].
After the balanced detection and the ADC, the down-converted comb consists of M =49 comb lines, spaced at ∆f = 50 MHz and centred at ∆f RF = 4.5 GHz. The power spectrum density (PSD) is shown in Fig. 1(b). The sampling frequency and the 3-dB bandwidth of the real-time sampling scope are F s = 50 GS/s and 16 GHz, respectively. The number of samples stored for processing and subsequent phase noise estimation is 12.5 · 10 6 . From the PSD, the extracted measurement noise variance is σ 2 n = 10 −3 corresponding to SNR avg = 19.6 dB. The noise measurement bandwidth is 50 MHz.
In Fig. 8, the estimated phase noise traces of frequency combs lines that have different powers (SNRs) are plotted. The evolution of the phase noise traces is important for gaining a qualitative insight into the phase stability and computation of the timing jitter as explained in [7]. We consider frequency comb lines that have 0 dB and -5 dBm of relative power. The first thing that should be noted is that the estimated phase noise traces for 0 and -5 dBm are very similar, as expected, when using the Bayesian filtering approach. This is clearly not the case when using the conventional method. The estimated phase noise traces using the conventional method show significantly more erratic behaviour, especially for -5dBm. This erratic behaviour is caused by the measurement noise. In Fig. 9, we plot the the variance of the differential phase noise defined as: where c is the line index of the central line. As an inset, we also show the estimated covariance matrices. Since the true variance, of the differential phase noise variance, and correlation matrix, are unknown, a quantitative computation of the error is not feasible. However, we do expect the variance of the differential phase noise to be a quadratic and smooth function of the line index as shown by numerical simulations in Fig. 5. Comparing Fig. 9 with the numerical one, a similar trend is observed. The curve obtained using the conventional phase noise estimation approach is on top of the curve obtained using Bayesian filtering. According to the numerical simulations, this is an indication of a bias with respect to the true curve. Moreover, we observe a more erratic behavior, when using the conventional approach, due to the measurement noise. Additionally, the erratic behaviour is also present in the correlation matrix estimated by the conventional approach. We can thereby, qualitatively, conclude that more accurate estimation of the differential phase noise variance and correlation matrix can be obtained using Bayesian filtering.
We also compared the extracted phase noise correlation matrices for different observation times τ obs , which correspond to the Fourier frequencies f obs = 1/τ obs . On Fig. 10, the correlation matrices calculated for both methods are shown for frequencies of f obs = 4 KHz and f obs = 5 GHz. It is possible to notice that in the low frequency regime the comb phases are all positively correlated, as at low frequency the combination of RF and Laser phase noise is the dominant source. In the high frequency regime instead, the correlation matrices appear close to diagonal, indicating that the phase noise of the comb is dominated by uncorrelated shot noise.
Next, we perform eigenvalue decomposition and plot the eigenvectors (v 1 , v 2 ) in Fig. 11. It is observed that v 1 is constant as a function line index indicating that (v 1 , λ 1 ) are associated with the CW laser. The eigenvector v 2 exhibits a linear evolution and is thereby associated with the RF oscillator. The behaviour of v 1 and v 2 is as expected for an EO comb. In Fig. 12, the eigenvalues λ 1 and λ 2 are plotted as a function of the observation time, τ obs . The evolution of λ 1 and λ 2 are proportional to the phase noise variances of the CW laser, σ 2 L,k , and the RF oscillator, σ 2 RF,k , respectively. For the simulations, the observation time, τ obs , down to 10 −9 s was considered, while for are experimental data we are limited to 2 · 10 −7 s. This is also one of the reasons why we do not see a big difference between the eigenvalues obtained using the Bayesian filtering and the conventional approach. For observation times below 10 −6 s, the evolution of λ 2 obtained by the Bayesian filtering approach is on top of the curve obtained using the conventional phase estimation approach. This observation is also observed for the numerical simulations, and we therefore tend to conclude that λ 2 obtained by the Bayesian filtering is closer to the true one.
It should be noticed that the evolution of eigenvalues λ 1 and λ 2 , as a function of τ obs , is different compared to the numerical simulation in Fig. 7. This implies that the phase noise statistics of the CW laser and the RF oscillator employed in the experiment are different compared to the simulations. An reason for this is likely that the Wiener process is not a good approximation of the CW laser and RF oscillator phase noise employed in the experiment.
To further investigate the nature of the phase noise processes, we employed the extracted eigenvectors (v 1 , v 2 ) to recover the laser and the RF oscillator phase noise from the detected phases of the comb lines. Let Φ be the matrix of dimension T × M containing the digitized phase noise samples of all the comb lines for the whole acquisition time. Using a generalization of (18), we can write where φ L 1:T , φ RF 1:T are the sources phase noise vector of dimension T × 1 and v 1 , v 2 are the corresponding eigenvectors of dimensions 1 × M. Exploiting the orthogonality of the eigenvectors, it is easy to show that the independent sources can be extracted from the whole phase noise matrix by using φ L 1:T = Φ · v 1 ||v 1 || −2 and φ RF 1:T = Φ · v 2 ||v 2 || −2 . In Fig. 13 we show the power spectral densities of the Laser and RF oscillator phase noise. As it can be seen from the spectra, these processes do not correspond to Wiener phase noise, as only the high frequency content follows the characteristic 1/f 2 power law. It can also be noticed that the Bayesian filtering technique is able to filter out the measurement noise, affecting the conventional noise characterization at high frequency.

Conclusions
A Bayesian filtering framework for joint amplitude and phase noise estimation of multiple frequency comb lines has been proposed and demonstrated for the first time. The framework is optimum in the mean square error sense, thereby providing the theoretically highest achievable accuracy. Significant improvements, compared to the state-of-the-art method, in terms of the accuracy of the estimated phase noise and correlation matrix as well as eigenvalues has been demonstrated numerically and experimentally. Most importantly, the proposed method provides accurate estimates over a wide range of SNRs. This is an important feature, as a wide range of frequency combs have large variations of SNR per comb lines. In summary, the flexibility and the optimally can promote the proposed method to become the new reference tool for comb noise characterization.