Abstract
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
1. 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–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 electro-optic 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.
2. Theoretical framework for joint amplitude and phase noise tracking
2.1 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 $\Delta 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\Delta 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}_{3\textrm {dB}}^{\textrm {rec}}$, i.e. $M\Delta f \leq \textrm {min}(F_s/2,{BW}_{3\textrm {dB}}^{\textrm {rec}})$, where $F_s$ is the sampling frequency of the ADC.

Fig. 1. (a) System set-up for numerical and experimental investigations. Two optical frequency combs are heterodyned together, and the donwconverted frequency comb is sampled and digitized, obtaining the sequence $y_k$. (b) Typical power spectral density of the downconverted comb. (c) Bayesian filtering framework for joint estimation of static, $\mathbf {Q}$, and dynamic parameters, $[{\phi }_k^1,\ldots ,{\phi }_k^M]$ and $[{a}_k^1,\ldots , {a}_k^M]$
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, $\phi ^m_k$ is given by the phase difference: $\phi ^m_k = \phi ^{m,\textrm {CUT}}_k - \phi ^{m,\textrm {LO}}_k \:$ with $m=1\cdots M$. The $m$-th line amplitude $a^m_k$ is given by the product: $a^m_k= 2R\sqrt {P_k^{m,\textrm {CUT}} P_k^{m,\textrm {LO}}} \:$ . $R$ is the responsivity of the receiver photodiodes, $P_k^{m,\textrm {CUT}}$ and $P_k^{m,\textrm {LO}}$ are the powers of the CUT and LO frequency comb $m$-th lines, respectively.
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, $\boldsymbol {y}_{1:T} := [y_1,\ldots ,y_T]$, where $T$ is the length of the measurement data in samples, the objective is to perform joint tracking of the down-converted frequency combs’ phases $\boldsymbol {\phi }_k := [\phi _k^{1},\ldots ,\phi _k^M]^\top$ and amplitudes ${\boldsymbol {a}}_k := [a_k^{1},\ldots ,a_k^M]^\top$ for $k = 1\cdots T$, where $\top$ 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. $[\boldsymbol {\phi }_k,{\boldsymbol {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:
2.2 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 $\phi ^m_k$ is governed by (1) and (2), and the measurements depend on $a^m_k$ and $\phi ^m_k$ via (3). Strictly speaking, given $y_k$, we would like to estimate amplitudes and phases that minimize the following MSE problem: $\sum _k\mathbb {E}[(\boldsymbol {x}_k-\boldsymbol {x}_k^{\textrm {true}})^2]$, where $\boldsymbol {x}_k:= [\boldsymbol {\phi}_k;{\boldsymbol {a}}_k]^\top$, $\mathbb {E}[\cdot ]$ denotes the statistical expectation and $\boldsymbol {x}^{\textrm {true}}:=[\boldsymbol {\phi}^{\textrm {true}}_k;{\boldsymbol {a}}^{\textrm {true}}_k]^\top$ 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 $\boldsymbol {x}^{\textrm {opt}}_k := [\boldsymbol {\phi}^{\textrm {opt}}_k;{\boldsymbol {a}}^{\textrm {opt}}_k]^\top$. According to the Bayesian filtering theory, the optimal estimate $\boldsymbol {x}^{\textrm {opt}}_k$ are obtained given by computing the expectation of the filtering distribution [13]:
Implementing the system of equations in (6) requires initial conditions for $\boldsymbol {x}_0$. Typically, $\boldsymbol {x}_0=\textbf {0}$ and $\mathbf {\Sigma }_0$ is initialized as a diagonal matrix. Here, $\mathbf {h}_{\boldsymbol {x}}(\cdot ) = [\mathbf {h}_{{\phi }}(\cdot ),\mathbf {h}_{A}(\cdot )]$ is the gradient of the measurement function (3) with respect to the hidden states $\boldsymbol {x}_k$. $\mathbf {h}_{\boldsymbol {x}}(\cdot )$ contains the gradient with respect to the phases $\mathbf {h}_{{\phi }}(\cdot )$ and the amplitudes $\mathbf {h}_{A}(\cdot )$ expressed as:
2.2.1 Static parameter estimation
The frequency spacing $\Delta f$ and the measurement noise variance $\sigma _n^2$ can be extracted by inspection of the power spectral density of the photocurrent prior to the Bayesian filtering. The matrix $\mathbf {Q}$, however, needs to be jointly estimated from the observed data $\boldsymbol {y}_{1:T}$ together with $\boldsymbol {x}_{1:T}$. The Bayesian treatment for the optimal estimation of $\mathbf {Q}$, in the MSE sense, is obtained by maximizing the posterior distribution of $\mathbf {Q}$ given the observations $\boldsymbol {y}_{1:T}$:
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 $\textrm {LL}^-(\mathbf {Q}) := -\log {p(\boldsymbol {y}_{1:T}|\mathbf {Q})}$. Then, $\textrm {LL}^-(\mathbf {Q})$ is the cost function we want to minimize. Typically, for a general state`space model $\textrm {LL}^-(\mathbf {Q})$ can be approximated as [13]:
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 $\mathbf {Q}^{\textrm {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 $\mathbf {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 $\boldsymbol {x}^{\textrm {s}}_k := [\boldsymbol {\phi}^{\textrm {s}}_k;{\boldsymbol {a}}^{\textrm {s}}_k]^\top$. 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 $\boldsymbol {x}^{\textrm {s}}_{1:K}$ with associated smoothed covariance sequence $\mathbf {\Sigma }_{1:K}^{\textrm {s}}$. Here $K$ represent the number of data samples used for parameter estimation, and typically $N \ll T$. The EKS is initialized by setting $\boldsymbol {x}^{\textrm {s}}_K = \boldsymbol {x}^{\textrm {opt}}_{K}$ and $\mathbf {\Sigma }^{\textrm {s}}_{K} = \mathbf {\Sigma }^{\textrm {opt}}_{K}$. Then, for every $k = K-1\cdots 0$ the following equations are computed
2.3 Covariance and correlation matrix computations
Noise characterization of the frequency comb can the be performed by computing phase noise and RIN spectrum using $\boldsymbol {\phi}^{\textrm {opt}}_{1:T}$ and ${\boldsymbol {a}}^{\textrm {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$:
3. 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 $\Delta \omega _m = 2\pi \Delta f_{\textrm {RF}} + 2\pi (m - \lfloor {M/2}\rfloor +1)\Delta f$, where $\Delta f_{\textrm {RF}} = 2.5$ GHz and $\Delta 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 $\Delta \nu _L =$ 1 kHz and $\Delta \nu _{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 - \lfloor {M/2}\rfloor +1$ is the relative line index, $\lfloor \cdot \rfloor$ is the floor function, $\phi ^{\textrm {L}}_k$ and $\phi ^{\textrm {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, per comb line: $\textrm {SNR}_m=\frac {\bar {a^2}^m/4}{\sigma ^2_nT_sf_{BW}}$, where $i=1,\ldots ,M$, $\bar {a}^m$ is the average amplitude of the line $m$ and $f_{BW}$ is the bandwidth in which the SNR is measured. An average SNR is then defined as: $\textrm {SNR}_{\textrm {avg}} = (1/M)\sum _m \textrm {SNR}_m$.In Fig. 2, the ability of the proposed framework to learn the system parameters is shown. We start with our guess (prior) on $\mathbf {Q}_{\phi }$ at iteration $k=0$, i.e. $\mathbf {Q}^{(0)}_{\boldsymbol {\phi}}$. For the considered case, we assume that $\mathbf {Q}^{(0)}_{\boldsymbol {\phi}}$ is a diagonal matrix. The true matrix denoted by $\mathbf {Q}^{\textrm {true}}_{\boldsymbol {\phi}}$ 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 $\mathbf {Q}^{\textrm {true}}_{\boldsymbol {\phi}} = 2\pi T_s[\Delta \nu _{\textrm {L}} + \mathbf {m}_r\mathbf {m}_r^\top \Delta \nu _{\textrm {RF}}]$, where $\mathbf {m}_r = [- \lfloor {M/2}\rfloor ,\ldots ,\lfloor {M/2}\rfloor ]^\top$. Now, as the number of EM iteration is increased, the estimated $\mathbf {Q}_{\boldsymbol {\phi}}$ approaches the true matrix $\mathbf {Q}^{\textrm {true}}_{\boldsymbol {\phi}}$. This is also indicated by the evolution of $\textrm {LL}^-(\mathbf {Q})$ which decreases as the number of EM iterations is increased. Even though our initial assumption on the initial guess $\mathbf {Q}^{(0)}_{\boldsymbol {\phi}}$ was wrong, the framework was still able to converge to the correct $\mathbf {Q}_{\boldsymbol {\phi}}$.

Fig. 2. (Numerical) (a) Evolution of the estimated $\mathbf {Q}_{\boldsymbol {\phi}}$ matrix during the EM learning algorithm. Each matrix is estimated after the $n$-th EM iteration. (b) Negative log–likelihood of the matrix $\mathbf {Q}_{\boldsymbol {\phi}}$ plotted for every EM iteration number.
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. $\textrm {tan}^{-1}(i^m_k/q^m_k)$. The resulting phase is processed through an unwrapping function that smooths every phase jump greater than $\pi$ by adding multiples of $2\pi$ to obtain smoother transitions and to make the phase as a continuous function of the time index. Finally, a linear detrending of the traces will produce the desired phase noise of the lines $\phi ^{m,\textrm {conv}}_{k}$ for $m=1,\ldots ,M$.

Fig. 3. Illustration of the conventional method for phase noise estimation of the donwconverted. The resulting phases are denoted as ${\phi }^{1,\textrm {conv}}_k,\ldots ,{\phi }^{M,\textrm {conv}}_k$.
In Fig. 4, we compare the accuracy of the conventional and the Bayesian filtering based approach for estimation of the correlation matrices when varying $\textrm {SNR}_{\textrm {avg}}$. The true correlation matrix is denoted by $\mathbf {\rho }^{\textrm {true}}_{\boldsymbol {\phi},T}$ and is depicted in left column in Fig. 4(a). The structure of $\mathbf {\rho }^{\textrm {true}}_{\boldsymbol {\phi},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.

Fig. 4. (Numerical) (a) True correlation matrices (left) shown together with the the estimated matrices obtained by the conventional method (center) and the Bayesian filtering method (right). The true correlation matrix is the same for all the simulations. Each row is a simulation with different value of average $\textrm {SNR}_{\textrm {avg}}$, obtained by averaging the SNR per line over all lines. (b) The RMSE between the true and the estimated phase noise traces, for all $M=49$ comb lines, as function of the $\textrm {SNR}_{\textrm {avg}}$. The linewidth under consideration is 1 kHz and 1 MHz.
The estimated matrices obtained by the conventional and the Bayesian method are denoted by $\mathbf {\rho }^{\textrm {conv}}_{\boldsymbol {\phi},T}$ and $\mathbf {\rho }^{\textrm {BF}}_{\boldsymbol {\phi},T}$, respectively. In Fig. 4(a) these are shown for the SNR$_{\textrm {avg}}$ of 10, 15 and 20 dB. The noise measurement bandwidth is 100 MHz. Figure 4(a) illustrates qualitatively that the estimated correlation matrices, $\mathbf {\rho }^{\textrm {BF}}_{\boldsymbol {\phi},T}$, are very similar to $\mathbf {\rho }^{\textrm {true}}_{\boldsymbol {\phi},T}$ for the considered ranges of $\textrm {SNR}_{\textrm {avg}}$. This is not the case for $\mathbf {\rho }^{\textrm {conv}}_{\boldsymbol {\phi},T}$. The difference between $\mathbf {\rho }^{\textrm {conv}}_{\boldsymbol {\phi},T}$ and $\mathbf {\rho }^{\textrm {true}}_{\boldsymbol {\phi},T}$ is especially pronounced for $\textrm {SNR}_{\textrm {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 $\textrm {SNR}_{\textrm {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 $\Delta \nu _{\textrm {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 $\Delta \phi _k^m = \phi _k^m - \phi _k^c$, where $c$ is the line index of the central line. For an EO comb, the variance of $\Delta \phi _k^m$ is expected to be a quadratic function of the line index due to line-dependent phase noise generation process, namely $\sigma _{\Delta \boldsymbol {\phi}_k}^2 = (m-\lfloor {M/2}\rfloor )^2\sigma _{\textrm {RF},k}^2$ [16]. We consider a relatively high $\textrm {SNR}_{\textrm {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 $\textrm {SNR}_{\textrm {avg}}$, results in an even more inaccurate computation of the differential phase variance using the conventional method.

Fig. 5. (Numerical) Variance of the differential phase, $\sigma ^2_{\Delta \phi _T^m}$, using Bayesian filtering (blue curve) and a conventional method (black curve). The red curve represents the ground truth.
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 $(\lambda _1, \lambda _2)$ and eigenvectors $(\boldsymbol {v}_1, \boldsymbol {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, $(\lambda _1,\boldsymbol {v}_1)$ and $(\lambda _2,\boldsymbol {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 ($\boldsymbol {v}_1$) and ($\boldsymbol {v}_2$) are showing in Fig. 6. We can therefore conclude that the eigenvalue $(\lambda _1)$ and eigenvector $(\boldsymbol {v}_1)$ are associated with the RF oscillator.

Fig. 6. (Numerical) Normalized eigenvectors corresponding to the two principal eigenvalues for the longest observation time $k=T$. The eigenvectors are extracted using a conventional, (Black curve), and a Bayesian filtering, (Blue curve), method. The true eigenvectors that associated with the CW laser and the RF phase noise are depicted in red.
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 $\boldsymbol {v}_1$ and $\boldsymbol {v}_2$. Since the phase noises of both laser and RF oscillator are modelled as Wiener processes, the variances increase with the observation time.

Fig. 7. (Numerical) Evolution of two main eigenvalues as a function of the observation time ($\tau _{\textrm {obs}}$). We also show the evolution of the variance of the CW laser and the RF oscillator phase noise, $\sigma ^2_{\textrm {L},\tau _{\textrm {obs}}}$ and $\sigma ^2_{\textrm {RF},\tau _{\textrm {obs}}}$, respectively.
For shorter observation times $(\tau _{\textrm {obs}}< 10^{-8} s)$, there is a discrepancy between $(\lambda _1, \lambda _2)$ and $(\sigma ^2_{\textrm {RF}}\left \|\mathbf {v}_1 \right \|^2,\sigma ^2_{\textrm {L}}\left \|\mathbf {v}_2 \right \|^2)$. However, the discrepancy is more significant for $(\lambda _1, \lambda _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.
4. 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 $\Delta f=$ 50 MHz and centred at $\Delta f_{\textrm {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\cdot 10^6$. From the PSD, the extracted measurement noise variance is $\sigma _n^2 = 10^{-3}$ corresponding to $\textrm {SNR}_{\textrm {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.

Fig. 8. (Experimental) Evolution of the estimated phase noise as a function of time using for comb-lines with different relative powers. (a) Comb line number $m=3$ with 0 dBm of relative power. (b) Comb line number $m=7$ with -5 dBm of relative power.
In Fig. 9, we plot the the variance of the differential phase noise defined as: $\Delta \phi _k^m = \phi _k^m - \phi _k^c$, 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.

Fig. 9. (Experimental) Variance of the differential phase, $\sigma ^2_{\Delta \phi _T^m}$, using Bayesian filtering (blue curve) and a conventional method (black curve). (b) estimated correlation matrices.
We also compared the extracted phase noise correlation matrices for different observation times $\tau _{\textrm {obs}}$, which correspond to the Fourier frequencies $f_{\textrm {obs}} = 1/\tau _{\textrm {obs}}$. On Fig. 10, the correlation matrices calculated for both methods are shown for frequencies of $f_{\textrm {obs}} =$ 4 KHz and $f_{\textrm {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.

Fig. 10. (Experimental) Correlation matrices of the phase noise traces estimated for different Fourier frequencies. (a) and (c) show the correlation corresponding to $f_{\textrm {obs}} = 5$ GHz for the conventional and Bayesian filtering framework respectively. (b) and (d) show the correlation corresponding to $f_{\textrm {obs}} = 4$ kHz for the conventional and Bayesian filtering framework respectively.
Next, we perform eigenvalue decomposition and plot the eigenvectors $(\boldsymbol {v}_1, \boldsymbol {v}_2)$ in Fig. 11. It is observed that $\boldsymbol {v}_1$ is constant as a function line index indicating that $(\boldsymbol {v}_1, \lambda _1)$ are associated with the CW laser. The eigenvector $\boldsymbol {v}_2$ exhibits a linear evolution and is thereby associated with the RF oscillator. The behaviour of $\boldsymbol {v}_1$ and $\boldsymbol {v}_2$ is as expected for an EO comb.

Fig. 11. (Experimental) Normalized eigenvectors corresponding to the two principal eigenvalues for the longest observation time $k=T$. The eigenvectors are extracted using a conventional, (Black curve), and a Bayesian filtering, (Blue curve), method.
In Fig. 12, the eigenvalues $\lambda _1$ and $\lambda _2$ are plotted as a function of the observation time, $\tau _{\textrm {obs}}$. The evolution of $\lambda _1$ and $\lambda _2$ are proportional to the phase noise variances of the CW laser, $\sigma ^2_{\textrm {L},k}$, and the RF oscillator, $\sigma ^2_{\textrm {RF},k}$, respectively. For the simulations, the observation time, $\tau _{\textrm {obs}}$, down to $10^{-9}s$ was considered, while for are experimental data we are limited to $2\cdot 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.

Fig. 12. (Experimental) Evolution of two main eigenvalues as a function of the observation time $\tau _{\textrm {obs}}$.
For observation times below $10^{-6}s$, the evolution of $\lambda _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 $\lambda _2$ obtained by the Bayesian filtering is closer to the true one.
It should be noticed that the evolution of eigenvalues $\lambda _1$ and $\lambda _2$, as a function of $\tau _{\textrm {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 $(\boldsymbol {v}_1, \boldsymbol {v}_2)$ to recover the laser and the RF oscillator phase noise from the detected phases of the comb lines. Let $\mathbf {\Phi }$ be the matrix of dimension $T \times 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

Fig. 13. (Experimental) Power spectral density of the extracted phase noise traces using the Eigenvector projection with a conventional method (black curve) and with the Bayesian filtering framework (blue curve). (a) shows the extracted laser source phase noise and (b) the RF oscillator phase noise.
5. 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.
Funding
European Research Council (771410, 771878); Swedish Research Council (2016-06077).
Acknowledgments
This work is supported by the European Research Council through the ERC-CoG FRECOM and DarkoComb project (grant agreement no. 771878 and 771410), and Swedish Research Council (2016-06077).
Disclosures
The authors declare no conflicts of interest.
References
1. J. Kim and Y. Song, “Ultralow-noise mode-locked fiber lasers and frequency combs: principles, status, and applications,” Adv. Opt. Photonics 8(3), 465–540 (2016). [CrossRef]
2. R. Schmeissner, J. Roslund, C. Fabre, and N. Treps, “Spectral noise correlations of an ultrafast frequency comb,” Phys. Rev. Lett. 113(26), 263906 (2014). [CrossRef]
3. X. Xie, R. Bouchand, D. Nicolodi, M. Giunta, W. Hänsel, M. Lezius, A. Joshi, S. Datta, C. Alexandre, M. Lours, P. A. Tremblin, G. Santarelli, R. Holzwarth, and Y. Le Coq, “Photonic microwave signals with zeptosecond-level absolute timing noise,” Nat. Photonics 11(1), 44–47 (2017). [CrossRef]
4. P. J. Winzer, D. T. Neilson, and A. R. Chraplyvy, “Fiber-optic transmission and networking: the previous 20 and the next 20 years,” Opt. Express 26(18), 24190–24239 (2018). [CrossRef]
5. L. Lundberg, M. Karlsson, A. Lorences-Riesgo, M. Mazur, V. Torres-Company, J. Schröder, and P. Andrekson, “Frequency comb-based wdm transmission systems enabling joint signal processing,” Appl. Sci. 8(5), 718 (2018). [CrossRef]
6. L. Lundberg, M. Mazur, A. Mirani, B. Foo, J. Schröder, V. Torres-Company, M. Karlsson, and P. A. Andrekson, “Phase-coherent lightwave communications with frequency combs,” Nat. Commun. 11, 201 (2020). [CrossRef]
7. A. Schlatter, S. Zeller, R. Paschotta, and U. Keller, “Simultaneous measurement of the phase noise on all optical modes of a mode-locked laser,” Appl. Phys. B 88(3), 385–391 (2007). [CrossRef]
8. D. Zibar, H.-M. Chin, Y. Tong, N. Jain, J. Guo, L. Chang, T. Gehring, J. Bowers, and U. Andersen, “Highly–sensitive phase and frequency noise measurement technique using bayesian filtering,” IEEE Photonics Technol. Lett. 31(23), 1866–1869 (2019). [CrossRef]
9. I. Coddington, N. Newbury, and W. Swann, “Dual-comb spectroscopy,” Optica 3(4), 414–426 (2016). [CrossRef]
10. V. Durán, S. Tainta, and V. Torres-Company, “Ultrafast electrooptic dual-comb interferometry,” Opt. Express 23(23), 30557–30569 (2015). [CrossRef]
11. G. Brajato, L. Lundberg, V. Torres-Company, and D. Zibar, “Optical frequency comb noise characterization using machine learning,” Presented at the European Conference on Optical Communication (ECOC), Dublin, Ireland, 22–26 Sept. 2019.
12. A. Ishizawa, T. Nishikawa, A. Mizutori, H. Takara, A. Takada, T. Sogawa, and M. Koga, “Phase-noise characteristics of a 25-GHz-spaced optical frequency comb based on a phase-and intensity-modulated laser,” Opt. Express 21(24), 29186–29194 (2013). [CrossRef]
13. S. Särkkä, Bayesian Filtering and Smoothing, vol. 3 (Cambridge University, 2013).
14. C. J. Wu, “On the convergence properties of the em algorithm,” Ann. Statist. 11(1), 95–103 (1983). [CrossRef]
15. L. Lundberg, M. Mazur, A. Fülöp, V. Torres-Company, M. Karlsson, and P. A. Andrekson, “Phase correlation between lines of electro-optical frequency combs,” in Conference on Lasers and Electro-Optics, (Optical Society of America, 2018), p. JW2A.149.
16. D. R. Carlson, D. D. Hickstein, W. Zhang, A. J. Metcalf, F. Quinlan, S. A. Diddams, and S. B. Papp, “Ultrafast electro-optic light with subcycle control,” Science 361(6409), 1358–1363 (2018). [CrossRef]
17. F. Riehle, Frequency Standards: Basics and Applications (John Wiley & Sons, 2006).