Digital signal processing techniques for noise characterization of lasers and optical frequency combs: A tutorial

Optical frequency combs are indispensable tools in many scientific and engineering disciplines.^1,2 Over the past two decades, the field of optical frequency combs has witnessed significant advancements, leading to the development of various types of frequency combs with distinct characteristics and applications.^3–6 The development of broad and stable combs has led to their integration in various fields, notably playing vital roles in the establishment of ultra-stable frequency references for metrology applications, high-capacity coherent telecommunication, and high-precision spectroscopy.^7–10 

While extensive research efforts have been made to realize various types of optical frequency combs, the development of novel noise characterization techniques has not followed a similar trend. Performing noise characterization of optical frequency combs poses a challenge due to their complex phase noise dynamics and the corresponding scaling with comb-line number.^11–15 Additionally, different comb lines may have different signal-to-noise ratios (SNRs), requiring measurement techniques to estimate phase noise over a wide range of SNRs accurately.¹⁶ 

Understanding and characterizing the noise properties of optical frequency combs is paramount to harnessing their full potential. The magnitude of their phase noise is a crucial property, as its presence degrades the frequency stability of the individual comb lines, resulting in limited system performance.^1,17

Various phase noise measurement techniques have been proposed in the literature.^18–20 As these methods are implemented in the analog domain, they introduce additional complexity to the experimental setup. More importantly, analog methods usually rely on optical bandpass filters to isolate single comb lines from the spectrum. In the aforementioned approaches, the measurement accuracy can be limited by the SNR of single comb lines. This is because the combination of low-power comb lines and measurement noise may lead to low SNR, leading to an inaccurate phase noise estimate. While cross-correlation methods can reduce the impact of measurement noise on the phase noise estimate, they require additional calibration and introduce complexity.²¹ Additionally, optically filtering out single lines can be infeasible for very narrowly spaced combs.

Another drawback of performing phase noise characterization of frequency combs based on line-by-line measurement is the inability to extract the time-domain correlations between neighboring lines. To remedy this, analog measurement methods have been developed to measure multiple comb lines simultaneously.^22,23 However, these methods again introduce more experimental complexity while still being limited to measuring only a few lines.

The complementary approach of using balanced detection in combination with digital signal processing for phase noise estimation shifts most of the measurement complexity into the digital domain. The availability of ultra-wideband photodetectors and analog-to-digital (ADC) converters at the order of 100 GHz bandwidth allows the detection of many comb lines simultaneously in multi-heterodyne measurements. The sampled data are stored offline, allowing advanced digital signal processing (DSP) for joint phase noise estimation of multiple frequency comb lines to be efficiently implemented.^7,24

Applying DSP techniques thereby allows simultaneous phase noise characterization of multiple comb lines while retaining a minimal measurement setup.²⁵ The extracted phase noise can then be used to compute the phase noise correlation matrix, which can then be used to identify the underlying phase noise sources and their comb line number dependent scaling.^11,16

However, some of the aforementioned DSP-based phase noise characterization methods employ digital bandpass filtering of individual comb lines, thereby inheriting some of the limitations of analog techniques. This implies that the SNR of the individual lines limits the estimation accuracy of each comb line’s phase noise, and narrowly spaced combs still pose an issue.

Bayesian filtering, implemented as an extended Kalman filter (EKF), offers joint digital signal processing of multiple comb lines. It eliminates the shortcomings of digitally filtering individual comb lines while maintaining the reduced experimental complexity. In particular, the EKF can provide near-optimal phase noise estimation in the presence of measurement noise.²⁶ Given an analytical description (model) of the detected signal and the associated phase noise sources, the EKF can distinguish between phase and measurement noise. This distinction allows the digital suppression of measurement noise and an improved phase noise estimate over state-of-the-art conventional methods.¹⁶ The EKF inherently provides confidence intervals (CIs) with the phase noise estimates, allowing the estimate’s uncertainty to be quantified. In a laser phase noise characterization, the EKF was able to predict phase noise up to 7.6 dB below the measurement noise floor with 95% confidence.²⁷ 

In Ref. 16, a framework based on the EKF is presented to characterize the phase noise of an electro-optic (EO) comb. The EKF operates in the time domain using a joint model of the detected comb lines, so no prior bandpass filtering is required. The estimation accuracy is improved compared to the conventional DSP methods, and the phase noise correlations can be recovered. In this framework, however, the computational complexity scaling of the EKF is $O(N3)$⁠, where N is the number of detected comb lines. This scaling is unfavorable and practically disqualifies the method for combs with many lines.

In Ref. 28, the EKF was used to compensate for phase noise in a dual-comb spectroscopy setup. Here, the complexity scaling issue was avoided by estimating the carrier-envelope frequency and the repetition rate phase noise of the comb instead of estimating each line. The complete reconstruction of a noise-free comb spectrum from a free-running comb is demonstrated. However, this method is tailored toward applying noise compensation in dual-comb spectroscopy rather than characterization.

This tutorial introduces an improved EKF-based framework for phase noise characterization of optical frequency combs. The framework’s significant improvement compared to Ref. 28 is the inclusion of automatic differentiation (AD). As the model within the framework may contain many free parameters that need to be determined, AD allows for gradient-based static parameter optimization. This approach to model parameter optimization significantly reduces complexity and enables efficient optimization in high-dimensional search spaces. Moreover, AD allows for the introduction of a latent phase noise space, a compressed representation of comb line phase noise. The structure of the latent space is learned directly from the measurement data. The introduction of an optimizable latent space not only solves the complexity scaling problem but also enables the search for phase noise modes beyond the commonly agreed upon noise sources in combs.¹¹ 

We apply AD and adaptive gradient-based optimization to learn the structure of the latent space. AD allows the optimization of an analytical model that describes the signal. Critically, AD enables the calculation of all model parameter gradients within any differentiable EKF model expression, eliminating the need for almost all auxiliary parameter estimation techniques. In this joint optimization process, the models of the measured signal and the statistical properties of the noise sources are learned. This versatility makes the method very general and, in principle, applicable to any comb.

In this work, we demonstrate in detail how to use this framework on different types of combs in simulation and experimental data. The remainder of this tutorial is structured as follows: In Sec. II, we introduce the required experimental setup and derive analytical expressions to describe the detected signal. We supply the necessary theoretical background for applying the EKF and optimizing a model. The framework is applied to a synthetic laser beat note measurement to detail every step in the process. Additionally, the state-of-the-art conventional phase estimation method is explained and applied. In Sec. III, the framework is applied to an EO comb on simulated and experimental data. Here, we demonstrate the versatility of the optimization process, the accuracy of the phase noise estimation, and the extraction of phase noise correlations. In Sec. IV, we show how the framework can be used in dual-comb measurements. Using two mode-locked lasers, we present the phase noise estimation process based on an experimental measurement. Furthermore, we discuss how this framework can be applied for digital noise compensation, similar to the method presented in Ref. 28. Finally, we give some concluding remarks and a brief outlook in Sec. V.

This section describes the complete framework in detail. First, the experimental setup common to all multi-heterodyne measurements considered in this work is described. A conventional phase noise measurement method using simple digital phase noise estimators is briefly explained. Finally, the EKF-based method is detailed. We describe how AD is applied in an efficient and versatile optimization framework. At the end of the section, we step through the example of a laser phase noise characterization, as it is conceptually simpler than the frequency comb case while still showcasing the intricacies of the method.

The experimental setup described in this work is shown in Fig. 1(a). The same setup is also used for simulation studies. It is used to generate heterodyne beat note signals for use in frequency comb and laser characterizations. The setup requires minimal analog components, including an optical coupler, a BPD, and an ADC, making it easy to reproduce.

Throughout this article, the same experimental setup is used for all characterization studies, whether in simulation or experiment. All optical components in the setup are assumed to be fiber-bound, but a free-space setup would also work.

In the measurement process, the comb or laser under test (CUT/LUT) and the local oscillator (LO) are combined using a balanced fiber coupler. The combined signal is then detected by a balanced photodetector, which suppresses intensity noise and reduces shot noise contributions, resulting in an improved signal-to-noise ratio (SNR) of the signal.²⁹ 

The detection process produces one or multiple beat notes, depending on the number of optical tones in the light sources and their detuning with respect to each other. This setup is known as (multi-) heterodyne detection.

This general description can model any combination of CUT/LUT and LO. Mathematically, a LUT is equivalent to a CUT with only a single line. As this is a very general signal description, it can be unwieldy and impractical. We will, therefore, examine multiple special cases of this expression throughout this article and detail how this expression simplifies.

If a single photodiode were used instead of a BPD, the detected signal would be proportional to |E_a|². This would result in constant power offset terms and no amplitude noise suppression. Additionally, single photodiode detection has a 3 dB SNR penalty compared to BPD. While BPD is the superior detection method, single PD detection can be used to the same effect if amplitude noise and SNR are not a concern.

The quantities of interest for the phase noise characterization of the CUT are $ϕlCUT(tk)$⁠. However, as both $ϕlCUT(tk)$ and $ϕmLO(tk)$ are randomly fluctuating phase noise terms, it is not possible to separate their difference $ϕlmtotal(tk)$⁠. As we cannot rely on additional prior knowledge about their statistics, they become indistinguishable. To still be able to estimate $ϕlCUT(tk)$⁠, one of two approaches should be chosen. For the first option, a LO needs to be selected whose phase noise is negligible compared to the CUT, in which case $ϕlmtotal(tk)≈ϕlCUT(tk)$⁠. The other option is to use a reference with close to identical phase noise properties as the CUT. It is still impossible to separate their contributions in the time domain, as the time evolution of their phases is uncorrelated. However, both noise contributions will have the same statistical behavior, which leads to them having approximately identical phase noise power spectral densities (PSDs) $SϕCUT,l(f)≈SϕLO,l(f)$ per line l. As the total phase noise $ϕlmtotal(tk)$ is the difference between the individual noise sources, the total PSD can be expressed as $Sϕtotal,ll(f)=SϕCUT,l(f)+SϕLO,l(f)≈2SϕCUT,l(f)$⁠. Here it is immediately clear that in this specific case, a division by two yields the PSD of the lth line of the CUT. Notably, this argument is only effective when investigating the beat notes that originate from comb lines with the same properties, as indicated by the identical line index l. This scenario is useful when multiple copies of the CUT are available.

If neither of these options can be chosen due to a lack of suitable references, the measured phase noise will contain an unknown mixture of both sources. While this is not ideal for the characterization of the CUT, extracting the total phase noise per comb line can still be useful in certain applications. In dual-comb spectroscopy, for example, the amount of total phase noise per detected comb line can be a limitation.³⁰ Here, the exact composition of $ϕlmtotal(tk)$ from the CUT and LO is less important in practice.

In this section, a brief summary of the conventional digital phase noise estimation algorithm is given.^11,16,26 It is based on the discrete Hilbert transform $H$ to obtain the second quadrature of the detected multi-heterodyne beat, followed by an inverse tangent operation to obtain the phase. A schematic illustration of the algorithm is given in Fig. 1(c).

Starting from the detected multi-heterodyne signal as in Eq. (4), each beat note needs to be isolated. This is accomplished using a bank of digital bandpass filters. Each filter is centered around the beat note frequencies ω_lm. The bandwidth is determined by the frequency difference between the detected notes. Applying the filter bank produces a time-domain trace y_lm(t_k) for each beat note. Each y_lm(t_k) now only contains one beat note and represents one of the quadratures of the signal.

Notably, this method does not attempt to remove measurement noise present in the signal. The measurement noise term ξ(t_k) in Eq. (4) is still present in each y_lm(t_k) after the bandpass filtering operation. In high SNR regimes where measurement noise is negligible, this method has been shown to be an optimal estimator of $ϕlmtotal(tk)$⁠.²⁶ In moderate to low SNR regimes, however, measurement noise can dominate such that the estimator is no longer optimal. The resulting phase noise estimate then contains measurement noise contributions, as the estimator cannot distinguish between the two. The resulting phase noise estimate is biased and overestimates the phase noise magnitude. In regimes with low SNR and low phase noise, the method can even fail completely to produce a reasonable estimate.

Additionally, the bandpass filtering operation limits the maximum offset frequency at which the phase noise can be estimated. This frequency is half of the bandpass filter’s bandwidth, which is half the frequency distance of a given beat note to its neighboring beat note. This limitation becomes relevant for narrowly spaced combs. In dual-comb spectroscopy, the multi-heterodyne beat notes can have spacings of less than 100 Hz. This method would, therefore, only allow a phase noise characterization of up to 50 Hz.

Throughout this work, this method will be called the conventional method and serve as a benchmark for the EKF framework.

In signals with low SNR or low phase noise beat notes, a method that can reduce the impact of measurement noise and isolate phase noise is necessary. For this purpose, Bayesian filtering in the form of the extended Kalman filter (EKF) has been proposed. The EKF has been shown to provide near-optimal phase noise estimates in (multi-) heterodyne measurements in the presence of measurement noise.^16,26 After formulating a state-space model that describes the detected signal in Eq. (4), it can be used to track dynamic quantities inside the model over time that cannot be measured directly. These are the phase noise terms and are often called the “hidden” state of the state space model. Simultaneously, the EKF attempts to isolate the measurement noise present in the signal based on its statistical properties. By separating the measurement noise from the phase noise, an accurate phase noise estimate can be produced even in low SNR regimes.

If the analytic expression underlying the detected signal was linear, the regular Kalman filter can be shown to provide the optimal estimate with regard to the mean squared error.³¹ In our case, however, the expression in Eq. (4) is highly non-linear. This is known as the non-linear filtering problem, which does not have an analytical solution. This implies that no optimal filter exists for our purpose. The EKF provides an approximate solution by locally linearizing the state-space model. While the optimality property is lost, this first-order approximation is typically sufficient to produce accurate estimates. We, therefore, call the estimate “near-optimal.”

The EKF operates directly on the detected signal in the time domain given in Eq. (4). All beat notes in the signal are simultaneously part of the state-space model, which makes this a joint estimation approach. This means that no prior bandpass filtering around the individual combines is required as a joint estimation of the phase noise of each beat note line is performed. Consequently, the method does not exhibit the same limitations induced by bandpass filtering as the analog methods or the conventional method.

In the following subsections, the general theory of the EKF is presented. A schematic representation of the algorithm is presented in Fig. 1(b). After a rundown of the prerequisites, the EKF is applied to characterize the phase noise of a laser in Sec. II D.

The state-transition function f(x_k, θ) models how the hidden variables evolve over time. In addition to the deterministic function f, the state space model includes an additive stochastic noise term q_k. This so-called process noise is a random variable distributed as a zero-mean Gaussian with a covariance matrix Q. This allows the modeling of hidden variables based on their statistical properties. In the multi-heterodyne case, Eq. (6) models the phase noise evolution.

Once a state space model of the form given in Eqs. (6) and (7) that describes the detected signal has been defined, the EKF can be applied. Applying the EKF consists of iteratively applying the so-called filtering equations. Given an initial value for the phase noise x₀, applying the EKF to the first signal sample y₁ produces an estimate for the phase noise x₁ at the first sample. Recursively stepping through the equations produces a vector of phase noise estimates x_1:T for each time step, where T is the number of samples in the signal.

Furthermore, the EKF not only produces the point estimates x_1:T. The EKF estimate for each sample is a Gaussian distribution; therefore, a vector of covariances P_1:T is produced alongside the point estimates. For each time step t_k, the EKF produces a Gaussian probability distribution $N(xk,Pk)$⁠, which describes the likelihood for a certain value of the estimate. This means that the phase noise prediction inherently includes a quantification of the uncertainty of the estimate. At each t_k, the standard deviation of the Gaussian estimate can be interpreted as a confidence interval.

Iterating through these equations and saving the estimates x_k and their covariances P_k then produces a series of phase noise estimates, including uncertainties that are based on the presence of measurement noise. The equations are given here for completeness; for an in-depth introduction and derivations, see Ref. 31.

Thus far, we have assumed that the values of the static parameters θ are known. This is generally not the case, and their values need to be determined. This step is crucial, as accurate values for θ are a strict prerequisite for accurate phase noise estimation.

In the EKF framework, this is most commonly performed by maximizing the posterior distribution p(θ|y_1:T) of θ given the observations y_1:T. It represents the probability of a specific value of θ being responsible for producing the observations y_k. If there is no prior information about the optimal value of θ, the prior probability distribution p(θ) is constant. According to Bayes’s rule, maximizing the posterior distribution then becomes equivalent to maximizing the likelihood p(y_1:T|θ).

A common technique to solve the optimization problem in Eq. (11) is the expectation–maximization algorithm. It is an iterative optimization scheme in which convergence to a (local) optimum is theoretically guaranteed. For certain state-space models, including the one used in Ref. 16, Eq. (11) can even be solved using a closed-form analytic expression.³³ 

For state-space models with many static parameters θ, however, an analytic form to solve the optimization problem is not available. For high-dimensional θ, an efficient approach to solve Eq. (11) requires the explicit computation of the gradients ∂LL⁻/∂θ. While it is possible to calculate the gradients manually, it requires the symbolic propagation of the gradients through the filtering equations (9). This quickly becomes infeasible for state-space models with many parameters θ.¹⁶ To alleviate this complication, the use of automatic differentiation (AD) to calculate ∂LL⁻/∂θ has been proposed.³⁴ Using AD allows the accurate and efficient evaluation of partial derivatives of any differentiable function specified in a computer program.^35,36 It keeps track of all elementary operations inside the function to build a graph of the function composition. By defining the differentiation rules of these elementary operations, the gradient of a function can be calculated by propagating through the composition graph and recursively applying the chain rule of differentiation. The key advantage of AD is its ability to differentiate arbitrarily complex functions. Instead of the error-prone derivation of symbolic derivatives, AD handles the computation automatically.

To evaluate ∂LL⁻/∂θ, the AD is propagated through the recursive EKF equations (9) that are required to evaluate LL⁻. For our implementation of the EKF framework, we used the JAX library, which provides a fast and easy-to-use implementation of AD.³⁷ 

With a way to calculate ∂LL⁻/∂θ, we are able to use any gradient-based optimization method. A robust and efficient optimization algorithm is required that performs well in high-dimensional search spaces since θ can potentially have many free parameters. Therefore, we propose the use of modern adaptive gradient-based optimizers, which are commonly used in the optimization of large neural networks. These algorithms are designed for fast convergence in high-dimensional optimization problems. In particular, we will use the AdaBelief optimizer throughout this article, as it converged very quickly and accurately in testing.^37,38

The new parameters $θ̃$ that are exposed to the optimizer are unconstrained real numbers. Choosing appropriate transformation functions Λ then allows the actual model parameters θ to be scaled and constrained.

To illustrate this concept, we explain a possible choice of optimization transform function for an exemplary parameter. The variance of the measurement noise R in the measurement Eq. (7) of the state-space model needs to be strictly non-negative. In addition, its value can become so small that numerical accuracy may become an issue. Both of these issues are avoided if Λ is chosen to be the exponential function, which maps the real numbers to non-negative real numbers. On top of that, small values close to zero are mapped to large negative values. This effectively causes the optimization to be performed in logarithmic space. Λ can be a different transformation for every parameter in θ, as long as it is a differentiable operation.

To reduce the computational load during optimization, a truncated version $LLtrain−=LL−(y1:Ntrain,θ)$ of the full negative log-likelihood can be used, which only considers N_train < T samples of the measurement vector y_1:T. Under the assumption that the measurement and phase noise are approximately stationary, this truncation is a good approximation of the full negative log-likelihood. A balance between computational load and accuracy has to be found. Considering too few samples will result in a biased estimate of the model parameters, while using a large number of samples will result in unnecessary slow convergence. However, the long-term behavior of the system cannot be captured by the likelihood if too few samples are considered.

LL⁻ is generally a non-convex function, which means that the optimization problem potentially has multiple local minima. The optimization algorithm can, therefore, get stuck at a local minimum and cannot guarantee a global optimum. To mitigate this problem, we typically run the optimization multiple times with different initial guesses for the model parameters θ and choose the model parameters that yield the highest likelihood.

To test whether a good optimum was found, the parameters $θ̂MAP$ can be tested on a section of the measured signal y that was not used to calculate $LLtrain−$ in the optimization phase. To test a set of parameters after convergence in the optimization phase, we simply evaluate $LLtest−=LL−(yNtrain:Ntrain+Ntest,θ)$⁠. N_test is the number of signal samples used and is usually chosen as N_test = N_train.

If $LLtest−$ evaluates to a similar value as $LLtrain−$⁠, this is a good indication that a set of parameters was found that is general over the whole signal and a good optimum was reached. A similar value is also a good indication that the choice of N_train is reasonable. It means that the optimized parameters θ are generalized and represent the stationary equilibrium of the signal.

To finalize the methodology section, we will step through a sample case of laser phase noise characterization in detail. The setup in Fig. 1 will be simulated for the case where the LUT and LO both are single-frequency lasers. Subsequently, we will apply the conventional and EKF-based phase noise estimation methods and compare their accuracy.

It should be noted that in the heterodyne beat note case, where the signal as expressed in (16) is just a single sinusoidal, it is also possible to demodulate the signal experimentally. This would result in a measured signal y(t_k), which is linear in ϕ^LUT(t_k) and would, therefore, not require the full EKF formalism. Instead, the (linear) Kalman filter could be used, simplifying the DSP. However, in the general multi-heterodyne case, this is no longer possible, which is why we introduce the non-linear EKF formalism here.

To generate a time series realization of ϕ^LUT(t_k) that follows the model in Eq. (17), we use the method given in Ref. 41. In this study, we choose Ξ_flicker = 10⁶ Hz² and Ξ_Lorentz = 100/π Hz to generate ϕ^LUT(t_k). Substituting into Eq. (16) produces the simulated signal y(t_k), which is displayed in Fig. 2, together with its PSD.

Here, x_k is the hidden state of the EKF and corresponds to the scalar-valued phase noise ϕ^LUT(t_k). r_k represents the scalar-valued measurement noise ξ(t_k) with variance σ_R.

With the state model and the parameter transformation function defined, the optimization loop illustrated in Fig. 1 is executed using the AdaBelief optimizer. To minimize the chance of sub-optimal convergence by running into a local optimum, 100 sets of randomly drawn initial parameters θ₀ are chosen. Each parameter set is independently optimized until the change in parameters per iteration is negligible. The convergence behavior for a selected subset of parameters is illustrated in Fig. 3. The initial values are distributed over a wide range of multiple decades to demonstrate the optimization routine’s robustness. Some of the optimization runs for which the random initialization of θ₀ is particularly poor do not converge. Most runs, however, converge to one of two local optima, which can be seen in the two distinct final values for the variances in Figs. 3(b) and 3(d). Ultimately, only the final parameters θ of the best-performing run are used.

Using the optimal parameters $θ̂MAP$⁠, we can apply the EKF to estimate the time evolution of the phase noise as the final step. We receive a time series of estimates x_k with a corresponding time series of variances $σP,k2$⁠. Together, they form a time series of Gaussian distributions $N(xk,σP,k2)$ of near-optimal estimations of the phase noise time evolution. To put the estimation into perspective, we apply the conventional phase estimation algorithm described in Sec. II A to the same simulated signal. Both methods are compared to the true phase noise in Fig. 4. The time domain traces shown in (a) indicate that the long-term behavior of the phase is evaluated accurately by both methods. Zooming in, however, reveals that the EKF estimate is more accurate than the conventional method. The frequency noise PSD picture in Fig. 4(b) shows that the low-frequency dynamics up to 300 kHz are accurately captured by both methods. For higher offsets, the conventional method is limited by the measurement noise floor ξ(t_k), which causes the estimate to have a ∝f² frequency dependency. The EKF estimate follows the ground truth below the noise floor, demonstrating that the measurement noise is effectively filtered out. The limits of the estimation accuracy are quantified by the variances $σP,k2$ of the estimates. Using second-order error propagation, the time domain uncertainty can be propagated to the PSD estimate, allowing us to draw a CI around the PSD. A derivation is presented in the  Appendix.

In summary, the presented method using the EKF produces a more accurate estimate than the conventional method. The introduction of adaptive optimization algorithms and the use of AD make the method very versatile, as state space models of arbitrary form with arbitrarily many parameters can be efficiently optimized. In Sec. III, we will demonstrate how this framework can be used to estimate phase noise in frequency combs.

The beat note generation process described in Eq. (21) is illustrated in Fig. 5. The resulting signal is once again a comb with the same spacing as the original line spacing. To simultaneously detect as many lines as possible, we want to maximize the number of down-converted lines within the detection bandwidth. This is achieved by selecting a LO with a suitable frequency ω^LO that is close to the frequency of the highest power comb lines of the CUT. All lines for which |ω_l| < 2πf_BW can then be detected. The maximum number of detectable lines within the detection bandwidth f_BW is N_detected ≤ 2⌊2πf_BW/ω_rep⌋ + 1. While we can always detect at least one comb line, this relation clearly illustrates that this setup is only suitable for combs that have a repetition frequency significantly below the available bandwidth. Only then can we detect multiple lines at the same time and benefit from the reduced experimental complexity of the setup to infer information about the comb phase noise from only a single measurement.

This condition disqualifies this setup for types of combs that have a repetition frequency ω_rep in the tens of GHz and above, as even state-of-the-art photodiodes are limited to around 100 GHz of bandwidth.⁴² Even with these extreme bandwidths, only a few comb lines would be detectable. Therefore, combs with large spacing can effectively only be characterized in dual-comb setups, which are discussed in Sec. IV.

Furthermore, this setup only has limited usability for combs with many lines on the order of thousands, which includes femtosecond laser-based combs or generally octave-spanning combs. Their repetition rate is typically on the order of GHz or less, and beating with a single laser will produce many detectable beat notes. However, a single measurement will only ever cover a limited range of lines in the comb under test. For a full characterization of all lines, multiple multi-heterodyne measurements with a wavelength-tuneable LO would be required to access the phase noise of all lines. While not infeasible, a dual-comb setup with a suitable reference comb is preferable in this scenario.

For combs with ω_rep up to 10 GHz and few comb lines, however, utilizing this single-frequency laser setup is more advantageous compared to employing a dual-comb configuration. With a single-frequency laser acting as the local oscillator (LO), it is still possible to detect numerous lines while the experimental setup is simplified significantly. Furthermore, assessing the contribution of the LO toward the total phase noise per beat note $ϕltotal(tk)$ is more straightforward. As derived in Eq. (21), the LO phase noise ϕ^LO(t_k) is an additive contribution independent of the comb line index l. It, therefore, only adds a common phase noise offset that does not add any index-dependent dynamics. This distinction is generally not as easy in a dual-comb setup, as the LO itself is a comb and its phase noise is, therefore, dependent on the line index.

To show how the EKF-based method can be used to extract the phase noise term $ϕltotal(tk)$ in this scenario, we will first demonstrate its application on a simulated signal. An electro-optic (EO) comb is chosen as the CUT, as it is well-studied in the literature and has a known phase noise model. After showcasing the method on simulated data, the method’s application is also demonstrated on experimental data.

The third phase noise term that appears in Eq. (24) is hidden in V(t_k), the output of the RF source. As Eq. (23) states, ϕ^RF(t_k) is the phase noise in the RF signal. RF source phase noise generally has different and more complicated characteristics compared to laser phase noise. Therefore, we will use a different phenomenological model to generate ϕ^RF(t_k).

Note that this model and its parameters are phenomenologically chosen to emulate the phase noise of a specific device, which we have used as Ref. 45. The secondary purpose of this model is its explicit deviation from a simple phase noise random walk. As discussed in Sec. II D, the EKF state-transition function we use to approximate the phase noise behavior does not consider any deviations from a random walk. An intentional mismatch between the true and approximated model can, therefore, give an indication of whether the random walk state-space model can still capture complicated phase noise dynamics.

To simulate time series realizations ϕ^CW(t_k) and ϕ^RF(t_k) based on the PSD models in Eqs. (17) and (25), we again use the algorithm given in Ref. 41. Equation (24) can now be used to directly simulate the time domain signal of a down-converted EO comb. The parameters chosen for this numerical study are $2ηPCWPLO$ = 100 mV, KV₀ = 1.6, ω^CW = 2π × 1.1 GHz, and ω^RF = 2π × 500 MHz. Here, we have combined some individual parameters into effective combined quantities, as only total values, as opposed to individual contributions, matter. The spectrum of the resulting signal y(t_k) is shown in Fig. 6(b). The resulting comb is centered around the line with index l_center = 2 and features N_alias = 3 and N_reg. = 8 aliased and non-aliased comb lines, respectively.

To use the EKF-based phase noise estimation method, a state-space model needs to be defined. The measurement equation of the model, which describes the detected signal, is based on the general analytical expression in Eq. (24).

Using this state-space model for EKF-based phase noise tracking has been shown to produce accurate estimations that outperform the conventional method in an EO comb.¹⁶ The approach does, however, have multiple drawbacks that have to be addressed.

In this individual line tracking approach, the correlations between the phase noise terms $ϕltotal(tk)$ are found through the optimization of Q_N. This matrix represents the differential phase noise covariance, from which the differential correlations can be calculated. No prior information about the phase noise in the system is used in this model.

The state vector x_k has the same dimension as the number of detected comb lines N. This causes the associated covariance matrices P_k as well as the process noise covariance matrix Q_N of the state vector to be a matrix of size $N×N$⁠. The EKF filtering equations (9) include matrix operations that have an undesirable computational complexity scaling of $O(N3)$⁠. Applying the EKF to a large number of comb lines N can make this approach computationally infeasible. Additionally, the covariance matrix Q_N needs to be optimized, which requires the optimization of $O(N2)$ free parameters.

To solve the complexity scaling problem, we can exploit the fact that the individual comb line phase noise terms $ϕltotal$ is a combination of a few underlying phase noise sources. The expression for the signal in Eq. (24) shows that only three phase noise sources contribute. The N individual beat note phase noise terms $ϕltotal(tk)$ are a combination of the three sources, which means that they are correlated.

The function Λ_C that transforms the mixing matrix is based on Eq. (30). Only a single free parameter, l_center, is required to define the matrix. The latent phase noise space is constrained to only have a common mode and a differential mode phase noise contribution. We will call this the constrained latent space model, indicated by the subscript C_c..

The constrained model has |θ_c.| = 2N + 3L + 4 free parameters and L tracked variables. For large N, the linear scaling of the number of parameters causes it to have significantly fewer free parameters compared to the individual line approach. Notably, the number of parameters no longer scales quadratically with the number of comb lines. Furthermore, the EKF computational complexity reduces to $O(L3)$⁠, which is no longer dominated by scaling with the number of comb lines. This means that it becomes feasible to use this model to apply the EKF-based phase noise estimation method to combs with a large number of comb lines.

Unconstrained latent spaces While the latent phase noise space of an EO comb can be explicitly described by Eq. (30), a similar expression might not be available for all types of combs. In particular, evidence for latent phase noise sources beyond the common mode and differential mode has been reported for mode-locked lasers.^11,47

In this form, the model has |θ_uc.| = NL + 2N + 3L + 3 parameters, which still has a favorable scaling compared to the individual line tracking approach as N becomes large. In addition, the complexity of applying the EKF is still $O(L3)$⁠.

To evaluate the latent space phase noise model, we apply it to the simulated EO comb discussed in Sec. III A. We will compare two variants of the latent space model against ground truth phase noise.

The constrained latent state-space model assumes that the comb line phase noise is explained by a common and a differential phase noise mode. This model is derived from the known behavior of EO comb phase noise as expressed in Eq. (28). The state-space model is defined in (32) and (33). Together with the optimization transition function, which is defined in Eq. (34), the parameters θ_c are optimized, as detailed in Sec. II C 3, and used to estimate the latent phase noise sources $ψ̂c(tk)$⁠.

For comparison, we also test the unconstrained state-space model variant on the same simulated signal. This variant uses the same state-space model. The difference lies in the optimization transformation function, where we use the replacement detailed in Eq. (35). While the latent space is still kept two-dimensional, it is not constrained to common and differential noise modes. Using the optimization framework, we find the parameters θ_uc. and use them to estimate the latent phase noise $ψ̂muc.(tk)$⁠.

Figure 7 displays the results of the EKF-based estimation. In Fig. 7(a), the columns of the optimized mixing matrix C_c. are compared to the true matrix. The single free parameter C_c. in Eq. (34) that identifies the index of the central comb lines was optimized. Its final value, $lcenterMAP≈2$⁠, is very close to the true value. Therefore, the columns of C_c. align perfectly with the true values. In (b), the frequency noise PSDs of the latent phase noise sources are shown. As C_c. matches C_true, the latent noise sources map directly to the physical phase noises $ϕtruecomm(tk)=ϕtrueCW(tk)−ϕtrueLO(tk)$ and $ϕtrueRF(tk)$⁠. At high offset frequencies f > 100 MHz, the EKF estimates deviate from the true PSDs. The CI around the PSD, however, captures this estimation inaccuracy.

In Fig. 7(c), the optimized mixing matrix C_uc. of the unconstrained model is shown. The column vectors do not display a well-defined separation into common and differential noise modes. The latent phase noise sources will, therefore, not correspond to physical noise sources. This can be seen in (d), where their PSDs are shown. The shape of the PSDs is similar to the true ones, but they deviate considerably. In addition, both latent phase noise PSDs display strong spikes at high frequencies, which are not present in true phase noise sources. The origin of these high-frequency oscillations is unclear, as they are not part of the state-space model. It is likely that due to the sub-optimal set of parameters θ_uc., the EKF produces erratic phase estimates to compensate for the accordingly inaccurate state-space model. This showcases the importance of finding good model parameters since phase noise estimates cannot be trusted otherwise.

While it can be useful to isolate the physical noise sources in a comb, this is only possible if a prior model of the comb phase noise exists, as in Eq. (28). In all other cases, especially in the unconstrained mixing matrix case, the optimizer will find any C that can approximately decompose the comb line phase noise into latent noise sources.

Here, the latent space is not interpreted as a physical quantity. However, the observable of interest is the phase noise per comb line $ϕltotal(tk)$⁠. The latent space simply serves as a compressed representation. Using Eq. (31), $ϕltotal(tk)$ can easily be computed from the latent phase noise.

In Fig. 8, the estimation accuracy of both the constrained and unconstrained EKF-based approach and the conventional method is compared against the ground truth. The phase noise per comb line can also be estimated using the conventional method discussed in Sec. II A. The accuracy is estimated based on the estimated variance of the phase noise.

In Fig. 8(a), the integrated phase noise PSDs of each individual comb line for the different methods are shown. This measure is similar to the empirical phase noise time domain variance as it has the same unit, except that it considers a wide part of the spectrum instead. The phase noise PSD of each comb line $Sϕtotal,l(f)$ is integrated from f_min = 100 kHz to f_max = 100 MHz. The value of f_max is chosen since 100 MHz is the maximum offset frequency we can evaluate using the conventional method in this case.

The unconstrained variant with C_uc. produces accurate estimates for the central comb lines, which have a high SNR. For the outer lines with SNR < 15 dB, the estimated variance is larger than the true value, which indicates that the phase noise magnitude is overestimated. The conventional method estimates are similar, where the high SNR comb line phase noise is estimated accurately. For the outer lines, however, the method completely fails. The very low SNR indicates that the measurement noise dominates over the beat notes, a regime in which the conventional method cannot produce a good estimate.

In Figs. 8(b) and 8(c), two exemplary phase noise PSDs of the comb lines with indices l = −1 and l = 3 are shown. The former beat note has a low SNR, which is, however, still enough such that the conventional method does not fail completely. Beyond f = 2 MHz, its PSD estimate is dominated by measurement noise. Both EKF variants follow the true PSD, which is below the noise floor. The beat note with index l = 3 is one of the highest SNR notes. The conventional estimate is considerably closer to the true PSD, as the noise floor is lower than that of line l = −1.

In Fig. 9, the correlation matrix estimates for the different phase estimation methods are shown. Both correlation matrices from the EKF-based variants are similar to the true matrix. In the unconstrained case, however, the correlations for the line with index l = 7 deviate from the ground truth. This line has a very low SNR, which causes problems in the optimization process. The optimized parameters that correspond to this line are harder to find, which causes the inaccuracy.

The conventional phase estimations’ correlations match the true correlations for the central lines. As discussed before, the method fails for low SNR lines, resulting in completely uncorrelated estimates for the outer comb lines. In summary, the accuracy of the correlation matrix estimated seems to follow the phase noise variance estimate that is shown in Fig. 8. This is not surprising, as inaccurate phase noise estimates will lead to inaccurate estimates of their correlations.

In this section, we demonstrate the robustness of the method and its applicability to real-world experimental data. The EKF-based method will be applied to experimental traces of a down-converted EO comb with similar characteristics to the simulated comb discussed in Sec. III A.

The experimental setup for the measurement considered in this section is identical to the setup displayed in Fig. 6(a). A NKT Photonics Koheras BASIK E15 fiber laser is mixed with an EO comb. The EO comb is generated using a separate NKT Photonics Koheras BASIK E15 laser and an EO modulator, driven by an Agilent E8247C signal generator. The mixed signal is detected by a balanced photodetector and digitized using a sampling oscilloscope with f_BW = 13 GHz and f_sampling = 40 GHz. After digitization, the recorded signal is downsampled digitally by a factor 3, as no beat notes appear in the high-frequency part of the signal spectrum. This step reduces the bandwidth of the measurement, which increases the SNR. The spectrum of the resulting measured signal is displayed in Fig. 10(a). We identify 13 non-aliased down-converted comb lines and 13 aliased lines, which are centered around the line with index l_center = 6. Notably, all lines except lines four to eight have low powers and, therefore, low SNRs. Additionally, each aliased line happens to be only about 36 MHz away from a neighboring non-aliased line. Since most of the regular lines have a higher SNR, they dominate over the aliased lines, which makes them hard to characterize. This is a problem when using the conventional method since it relies on bandpass filtering around each line. Because of the narrow spacing, the maximum bandwidth of each filter is limited to 36 MHz. This limits the maximum measurable phase noise PSD offset frequency to half the bandwidth. For this measurement, this means that the conventional method can be used to characterize each comb line only up to 18 MHz. While in practice, higher bandwidths could be used, the resulting phase noise estimate would be biased since it would include an aliased line.

This issue could be mitigated in this particular experimental setup by adjusting the central frequency of the seed laser. Nevertheless, in frequency combs with low repetition rates, the narrow spacing between lines remains an inherent challenge. Consequently, this measurement highlights the benefit of employing the EKF, as it eliminates the need for bandpass filtering and, therefore, circumventing the problem.

As long as the state space model correctly describes the detected signal, it does not matter how close the lines are to each other in frequency. We can, therefore, separate the contribution of each line and, in principle, isolate each line from the signal. However, here we are interested in their phase noise properties.

To apply the EKF to the detected signal, we use the same steps as in Sec. III A. The same state-space model as defined in Eqs. (32) and (33) is used, which describe the latent space phase noise model. As shown in Sec. III C, the constrained model variant produces a more accurate phase noise estimate, especially for very low SNR comb lines. In this section, we will, therefore, only use this variant. The corresponding optimization transformation function Λ is defined in Eq. (34). After optimizing the model parameters θ, the EKF is applied to estimate the phase. As we used the constrained model for the mixing matrix C_c., the latent phase noise sources can be identified as the common and differential mode noise sources. Here, the common mode corresponds to the difference between the phase noise of the comb seed laser and the LO ϕ^CW(t_k) − ϕ^LO(t_k). Both contribute equally to all down-converted comb lines, as in Eq. (28). The differential mode is the phase noise originating from the RF source ϕ^RF(t_k).

Both the common and differential mode contributions are separately accessible in the experiment, so we can record independent phase noise reference measurements. The common mode phase noise is estimated by creating a single beat note signal between the seed laser and the LO, which has the same total phase noise as the common mode of the EO comb. For the RF source reference phase noise measurement, the source is connected directly to the ADC. The phase noise of both reference measurements is estimated using the conventional method.

The frequency noise PSDs of the estimated latent phase noise sources are shown in Fig. 10(b), alongside the reference measurements. The PSDs of the EKF estimates are accurate over a wide frequency range, as they match the reference measurements very well. The common mode latent phase noise ψ₀ matches the seed laser phase noise over the full frequency range. Up to 1 MHz offset frequency, the differential mode phase noise ψ₁ matches the RF source phase noise. Beyond 1 MHz, the reference measurement is limited by measurement noise, which can be seen by the ∝f² dependency of the frequency noise PSD. At the same frequency, the PSD of ψ₁ flattens out and its CI grows larger. Still, the interval predicts the PSD to be significantly below the noise floor of the reference measurement. A flat frequency noise PSD corresponds to the random walk behavior that is part of the EKF state space model in Eq. (32). It is, therefore, not surprising that the EKF PSD estimate takes on this behavior: Once the phase noise is limited by measurement noise, extrapolating from the state-space model is the only remaining information the EKF can use. The CI does not reflect this model assumption. Its calculation, as presented in the  Appendix, does not include uncertainty in the parameters of the model itself. It represents uncertainty originating only from the measurement noise while assuming complete certainty about the state-space model. While a random walk is a valid assumption for laser phase noise, this is not necessarily the case for an RF source, as discussed in Sec. III A. The latent phase noise PSDs below the measurement noise floor should only be trusted if the phase noise behavior is well-modeled in the state-space model. Since we cannot assume that the random walk is a good model for the RF source phase noise, we cannot trust the PSD of ψ₁ beyond 1 MHz.

Next, we investigate the estimation accuracy with respect to the phase noise of the individual comb lines. Just as for the simulated EO comb, the phase noise per comb line can be calculated using Eq. (31). They are compared to estimations using the conventional method.

Figure 11(a) shows the integrated phase noise variance per comb line for the EKF and conventional method. The phase noise variance per comb line is calculated by integrating the phase noise PSDs of each comb line between f_min = 4 kHz and f_max = 15 MHz. For the central comb lines with SNR > 15 dB, both methods agree on the phase noise variance. As we saw for the simulated EO comb in Fig. 8, both methods produce accurate estimates in high SNR regimes. Therefore, we can assume that both methods produce an accurate estimate for lines four to eight. Beyond these lines, the conventional method fails. Since we have used the constrained EKF model, its estimate follows the theoretical phase noise model of the EO comb in Eq. (28). Just as in the simulation, the variance follows the expected parabolic form. Since this is an experimental measurement, we do not have access to ground truth values. The simulation results from Sec. III C, however, indicate an accurate prediction.

In Figs. 11(b) and 11(c), the estimated correlation matrices of both methods are shown, which are calculated according to Eqs. (37) and (40). As expected, the EKF predicts strongly correlated comb line phase noise. The conventional method can only recover the correlations for the few central lines that exhibit enough SNR.

To summarize this section, we have shown how to use the EKF to characterize the phase noise of down-converted frequency combs using the example of an EO comb. Both in simulation and on experimental data, we demonstrate that the EKF outperforms the conventional DSP method. Especially for very low SNR comb lines, the EKF can still recover the phase noise while the conventional method fails. While the EO comb was used as a benchmark system, the latent space EKF method is much more general. Two variants of the latent space were introduced, which, in principle, allow for its application on any comb.

In this section, we investigate the dual-comb setup. We continue to utilize the same experimental setup as displayed in Fig. 1. Both the CUT and the LO are chosen to be frequency combs. In this case, the detected signal has the most general form as derived in Eq. (4). Each optical comb line of the CUT can, in principle, create a beat note with any of the optical comb lines of the RS. The number of beat notes in the detected signal depends on the repetition frequencies $ωrepCUT$ and $ωrepLO$ of the comb, as well as the available detection bandwidth.

The beat notes in dual-comb setups can be categorized into different orders. One line of each optical comb is arbitrarily chosen to have comb line index zero. The electrical beat note created by their mixing is a zeroth order beat note and is assigned a line index from each generating line. Its frequency will, therefore, be denoted as $ω0,0=ω0LO−ω0CUT$⁠. All beat notes generated from the neighboring comb lines of each comb with matching index are also considered zeroth order. The index difference of the generating lines determines the order number of a beat note. A beat note with frequency $ω0,1=ω0LO−ω1CUT$ would, therefore, be considered a beat note of order one.

In Fig. 12, the beat note generation process is illustrated. The spectrum of the zeroth-order beat notes is again a comb in the electrical domain with repetition frequency $Δωrep=ωrepCUT−ωrepLO$⁠.