## Abstract

Dual-comb spectroscopy is a rapidly developing spectroscopic technique that does not require any opto-mechanical moving parts and enables broadband and high-resolution measurements with microsecond time resolution. However, for high sensitivity measurements and extended averaging times, high mutual coherence of the comb-sources is essential. To date, most dual-comb systems employ coherent averaging schemes that require additional electro-optical components, which increase system complexity and cost. More recently, computational phase correction approaches that enables coherent averaging of spectra generated by free-running systems have gained increasing interest. Here, we propose such an all-computational solution that is compatible with real-time data acquisition architectures for free-running systems. The efficacy of our coherent averaging algorithm is demonstrated using dual-comb spectrometers based on quantum cascade lasers, interband cascade lasers, mode-locked lasers, and optically-pumped microresonators.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Dual-comb spectroscopy (DCS) [1–3] derives its strength from the massively parallel heterodyne down-conversion procedure that enables direct mapping of the information encoded in the optical domain to the radio-frequency (rf) domain, where high-speed analog-to-digital converters (ADCs) can be used to acquire the signals. This requires spatial overlap of the beams from two matched frequency comb sources on a fast photodetector, causing optical beating frequencies spread over the detector bandwidth. Unlike most broadband techniques, this gives nearly instantaneous access to the available optical bandwidth, with applications in environmental monitoring, industrial process control and characterization of transient chemical dynamics, where the combination of spectral coverage and time resolution can be effectively leveraged. However, a stringent requirement is placed on the mutual coherence of the two comb sources, and any significant drift or fluctuation in phase will degrade the DCS system performance over longer time-scales. This effect can be mitigated if the phase noise can be actively suppressed through high-bandwidth hardware feedback loop stabilization [4–6], or measured continuously to phase-shift the acquired signals and adaptively adjust the ADC sampling rate [7]. If implemented successfully, the time-dependent distribution of the spectral power will be aligned in the frequency domain enabling averaging without dispersing the power over a large rf bandwidth, a process referred to as coherent averaging, also known in radar systems as coherent integration [8]. It is important to note that for some applications a straightforward alignment of multiple short-time DCS magnitude spectra [5, 9, 10] or time-domain interferogram bursts [11] may be sufficient. However, there are some inherent disadvantages of this simplified approach: the linewidths of the beat notes are limited by the acquisition time intervals and the lack of repetition rate correction does not account for the cumulative broadening of the beat notes. In addition, only beat notes identifiable within a single frame can be used, which precludes weaker beat notes especially those at the extremes of the spectral coverage, thus limiting the operational bandwidth.

The conventional approach for coherent averaging in DCS systems incorporates additional electro-optical hardware, such as external cw reference lasers and photodetectors [12–14]. While such an approach is feasible in the near-infrared where the telecommunication components can be leveraged, at longer wavelengths (mid- to far-infrared) availability of electro-optical components is very limited. Moreover, hardware solutions can be complex, cause reliability issues, and contribute to an increase in instrument cost. For this reason, we focused on the development of a purely computational coherent averaging methodology that is compatible with real-time architectures and acts directly on the digitized photodetector signal. This allows for correction of relative phase drifts and fluctuations even when the rf beat notes are masked by excessive phase noise. Although the procedure can be applied to any free-running DCS system, e.g. such based on electro-optic modulators, optically pumped microresonators, or fiber lasers, the main motivation for this work is linked to the recent developments in chip-scale frequency combs based on electrically-pumped semiconductor sources such as quantum cascade lasers (QCLs) [9, 15–19] or interband cascade lasers (ICLs) [20–22]. These sources have reached a state of development where they can be considered as the basis for miniaturized, battery-operated, portable spectrometers with small footprints [23]. However free-running semiconductor laser combs exhibit fluctuations in repetition rate (*f*_{rep}) on the order of hundreds of hertz, which in comparison to typical repetition rates of several gigahertz correspond to a relative stability of ∼10^{−7}. This is several orders of magnitude worse than that seen in fiber-based combs, which routinely reach relative stabilities of <10^{−9}. In addition to the repetition rate, the comb offset frequency (*f*_{0}), commonly referred to as carrier-envelope offset (CEO) frequency, fluctuates as well. These two sources of phase noise are directly down-converted to the rf domain, as schematically shown in Fig. 1.

Even though such combs show high mode-to-mode coherence, they are corrupted by considerable amounts of phase noise. When such free-running combs are spatially overlapped on a fast photodetector, rf beat notes with linewidths approaching MHz-level are observed. In the worst case, instead of evenly spaced down-converted narrow rf comb lines, the rf power will be dispersed over a broad rf bandwidth, as shown in Fig. 2(a). Clearly, a scenario like this is incompatible with DCS, but the apparent loss of information can be remedied via proper phase correction.

A well-known fact about the Fourier Transform states that it is not suitable for analysis of non-stationary signals, i.e. such that change their statistics over the acquisition time [24]. Beat notes whose frequencies fluctuate randomly (often much more than the spacing between the lines) are undoubtedly examples of such non-stationarity [see Fig. 2(a)]. Additionally, due the types of phase noise presented in Fig. 1, each rf beat note has a unique frequency deviation and hence unique phase noise characteristics. However, by taking advantage of the underlying comb nature of the laser sources, the beat notes can be effectively recovered from the noise. Recent works on digital correction algorithms [25,26] clearly show this capability. In 2016, Burghoff et al. presented a purely computational correction of DCS signals [25] based on the Extended Kalman Filter (EKF), which reconstructed the amplitudes, phases, and frequencies of all RF beat notes in the DCS spectrum. Although effective, the computational complexity of the EKF caused predominantly by the inversion of a large matrix, limits its real-time application to systems with less than a few dozens of modes.

Another approach proposed a year later by Hébert et al. involves the use of the cross-ambiguity function developed initially for radar applications [26]. It is related to the cross-correlation technique capable of retrieving the time delay between consecutive bursts of the interferogram and hence the quasi-instantanous repetition rate, while in addition the frequency offset between two similar waveforms can be obtained. An impressive correction of approximately three thousand comb teeth initially dispersed in the RF domain due to *f*_{0} and *f*_{rep} fluctuations was demonstrated [26]. However, the algorithm requires the generation of a two-dimensional map with multiple points for each burst of the interferogram, followed by a search for the global maximum at each step, which can be computationally demanding.

Compared to these approaches, our computational coherent averaging (CoCoA) algorithm provides similar functionality in terms of phase correction, but has been developed to meet the demands of real-time processing and is therefore implemented in a structurally different way. First, it retrieves the instantaneous repetition rate through a simple single-step non-linear operation (digital difference frequency generation, DDFG) followed by coherent demodulation of any self-mixing harmonic using an IQ demodulator, which is easily realizable in programmable hardware. The retrieved signal acts as an adaptive clock to correct the repetition rate fluctuations and keep the rf comb lines equidistant. Once this is ensured, we employ a simple recursive frequency tracker optimized for real-time tracking, which corrects for the global offset phase noise. Consequently, we correct for both sources of instabilities in the down-converted rf comb. Most importantly, the linewidth limitation of a single frame resolution encountered in magnitude spectra alignment does not pertain here. The linewidth of the corrected rf comb lines is mainly limited by the total acquisition time, provided the sources show sufficient mutual stability at a 1/Δ*f*_{rep} timescale, and that the extraction of the correction signal is not corrupted by the noise of the photodetector. Furthermore, the output of the CoCoA algorithm is a phase-corrected time-domain signal, which allows the use of ordinary digital filters for analyzing data in swept spectroscopy mode [27]. The algorithm shows almost no degradation of correction efficacy for higher frequency rf beat notes, ascribable to the use of high-order harmonics of the repetition rate, and its computational complexity does not depend on the number of comb lines in the spectrum. Finally, even weak beat notes can be used for spectroscopy because the entire duration of the beating signal can be harnessed for the calculation of the Fourier Transform. It should be noted that the CoCoA algorithm works directly on the acquired interferograms and does not correct the phases of the optical modes interacting with the sample. It merely assures mutual coherence by manipulating the timing and phase of the acquired photodetector signal, but for absolute optical frequency calibration an optical frequency reference is required.

## 2. Short-time Fourier Transform analysis of radio-frequency spectra

Before we proceed to an actual description of the CoCoA algorithm, let us analyze the effect of acquisition time on the shape of the acquired DCS spectrum as illustrated by the data shown in Fig. 2(a) acquired with 10 ms averaging. In this experiment we optically beat two cryogenically-cooled terahertz QCL phase-noisy combs on a 7.5 GHz bandwidth superconducting hot electron bolometer (Scontel), and acquire rf data by IQ-demodulating at 1.75 GHz using a real-time spectrum analyzer (R&S FSW-43). This example is selected to illustrate the strengths of the CoCoA algorithm in presence of strong phase noise. Also, the spectrum contains a relatively low number of rf beat notes, which make the effects of the CoCoA algorithm clearly visible (compare the panels in Fig 2). However, the algorithm works equally well for larger sets of beat notes, of which examples are given in Section 4, 5, and 6.

The effects of averaging in the presence of phase noise are illustrated in Fig. 3. The top panel of Fig. 3 shows the frequency spectrum acquired at 1 μs, wherein a dozen of strong beat notes is easily distinguishable from the large-variance noise floor. The 3 dB linewidth is close to the Fourier limit defined by the acquisition time (∼1 MHz), which confirms that the two combs are stable at a 1/Δ*f*_{rep} timescale, unfortunately accompanied by the typical uncertainty of the peak amplitude on the order of tens of percent [6, 20]. In other words, a transmission DCS measurement using such a short acquisition time would be associated with large error bars. Naively, one would expect that by increasing the acquisition time ten times (10 μs), a significant increase in the signal-to-noise ratio would be observed. Unfortunately, the opposite is true, which stems from a broadening of the beat notes that take highly asymmetric shapes. A further increase by an additional factor of ten (100 μs) causes the beat notes to disperse over tens of megahertz and the retrieval of peak amplitudes for spectroscopic assessments with such non-uniform spectra is a non-trivial task.

One way to analyze non-stationary signals with the Fourier Transform is to calculate it over shorter time scales when the considered signal is locally stationary, commonly referred to as Short-time Fourier Transform (STFT). A plot of the squared magnitude of the STFT as a function time is known as the spectrogram, which is shown in Fig. 3(b). The spectrogram is generated for 1 ms of acquisition time with a temporal resolution of 1 μs. The oscillatory fluctuations of the beat notes are caused by the noisy environment of the vibrating cryostat, as well as the presence of optical feedback in the system corrupting stable comb operation. It is important to note that in general the fluctuations are much more complex without a periodic character as shown in this particular case. This example clearly demonstrates that ordinary averaging fails to be effective when applied to DCS signals affected by phase noise. Therefore, it is imperative to apply the appropriate corrections in order to suppress the noise and increase the signal-to-noise for low-uncertainty amplitude estimation.

## 3. Phase and timing correction for coherent averaging

At the core of the CoCoA algorithm is the assumption that the comb sources in the DCS system are coherent on a 1/Δ*f*_{rep} timescale, but on longer timescales they are affected by repetition rate and offset frequency noises as was illustrated in Fig. 1. Therefore, it is sufficient to resample the signal to compensate for a variable duration of the consecutive beating interferograms (correction of the fluctuating repetition rate difference), followed by a global tracking and correction of the difference in offset frequency. A block diagram of the proposed solution is shown in Fig. 4.

The algorithm consists of two estimation and two signal manipulation blocks, included in dashed, and solid boxes, respectively. We will provide a general description of each signal manipulation block in the following subsections and additional details related to the signal processing performed in each block are provided in the Appendix.

#### 3.1. Tracking of the instantaneous repetition rate

For simplicity, consider a complex dual-comb signal *s̃*(*t*) consisting of *N* unitary amplitude components with an rf offset frequency Δ*f*_{0} and a rf repetition rate Δ*f*_{rep}.

This formula, referred here to as DDFG, produces a signal that is the result of the beating between all possible rf beat notes, i.e. a superposition of beating products between neighboring beat notes is obtained (every other, third etc.). Obviously, higher order harmonics of the fluctuating repetition rate will be more dispersed in frequency than the lower ones due the cumulative broadening effect, which increases linearly with harmonic number.

Once the DDFG signal is obtained, the instantaneous frequency *k*Δ*f*_{rep}(*t*) of the *k*-th harmonic of the repetition rate can be tracked with a high signal-to-noise ratio (in practice ∼20 dB). The tracking procedure can be realized through computational IQ demodulation at the expected value of the repetition rate harmonic within a fixed bandwidth, by employing filtering and the Hilbert transform, or by implementing more sophisticated harmonic frequency trackers [28], as described in Section 3.3. Frequency trackers are capable of handling overlapping beat notes, which is often encountered in combs with more than a hundred spectral modes and small repetition rate difference Δ*f*_{rep}.

In an FPGA platform, the rf repetition frequency tracking procedure can be implemented in a typical digital down-conversion (DDC) scheme in real time. A digitally controlled oscillator (DCO) feeds a digital mixer to down-convert an arbitrary repetition rate harmonic, followed by low-pass filtering the mixed signal using a fixed-bandwidth digital filter. The phase of the down-converted narrow-bandwidth signal is used to adjust the frequency of the DCO (over longer time scales), thus establishing an adaptive digital frequency/phase locked loop.

The DDFG operation also serves as an important diagnostic tool. A DCS signal that shows no repetition rate harmonics is non-correctable using the procedure described here because it does not possess comb characteristics. This useful diagnostic feature is discussed in Section 4.

#### 3.2. Correction of the variable repetition rate

Having retrieved the *k*-th harmonic of the instantaneous rf comb spacing *k*Δ*f*_{rep}(*t*), the time-domain signal needs to be resampled onto a non-uniform grid to cancel the effects of the repetition rate fluctuations. The time axis used for resampling has a cumulative character, and is given in the continuous case by

*f*

_{rep}〉 denotes the expected value of the repetition rate. The resampling operation aims to ensure an equal duration of the beating events in the interferogram. In this example, we are using the most trivial linear interpolation between two sampled points. Note that in the case of semiconductor combs equipped with a microwave bias tee, one can electrically extract individual intermode beat notes, mix them and feed the difference frequency to the data acquisition reference clock, thereby performing hardware adaptive ADC sampling (similar to ref. [7]) instead of this numerical procedure.

#### 3.3. Tracking of the instantaneous offset frequency

After the Δ*f*_{rep} correction that ensures constant spacings in the down-converted rf comb, the global offset frequency fluctuations need to be extracted. In principle, all beat notes in the rf spectrum can be tracked using a modified Multiple Frequency Tracker [29] or the Kalman Filter [25]. However, after the repetition rate correction, the complexity of this task can be significantly reduced to tracking of only one beat note isolated using a digital band-pass filter. We have identified a few different schemes for this purpose with varying complexity, yet with nearly the same tracking performance. The simplest, proposed by R. Aarts [30] (here referred to as the Fast Frequency Tracker), utilizes a recursion approach for a real discretely-sampled signal where the formula

*x*with a forgetting factor

_{k}*γ*. Of course, indices

*k*− 1 and

*k*− 2 denote values from previous time instances. Since

*r̂*is an auxiliary signal that contains the tracked frequency in the argument of the cosine where

_{k}*T*=1/

_{s}*f*is the sampling period, to estimate the instantaneous frequency Δ

_{s}*f*

_{0}at

*k*-th iteration we need to calculate the arccosine of

*r̂*scaled by the sampling frequency and 1/2π

_{k}*r̂*uses only multiplications and additions, making it compatible with real-time correction. In addition, the arccosine function can be tabulated to facilitate real-time offset frequency tracking.

_{k}One of the drawbacks of the Fast Frequency Tracker is a degradation of its performance when *x _{k}* is composed of numerous sine waves, which occurs when the drift of the offset frequency exceeds the spacing between the rf beat notes. We have found two solutions to effectively solve this problem. The first employs the aforementioned Multiple Frequency Tracker [29], which is proven to be robust while having a very low tracking delay, at the expense of computational complexity due to its use of matrix inversions. The second is inspired by carrier recovery systems in wireless networks, which comprise a coarse and fine carrier frequency tracker, where the fine tracker is simply the previously discussed Fast Frequency Tracker. The coarse estimation keeps the band-pass filter for

*x*at the pertinent frequency and is realized through the Short-time Fourier Transform and cross-correlation detecting the shift in the position of the magnitude spectrum by the same token as in [5]. In short, we reduce the problem of drifts exceeding the spacing between the rf lines to tracking residual fluctuations around almost stationary beat notes. The use of the STFT is favored due to its real-time compatibility with modern signal processing platforms. Alternatively, one can employ an adaptive tracking algorithm similar to that proposed in the repetition rate tracking step, which is based on a digitally controlled oscillator (DCO), fixed bandwidth low-pass filter and carrier offset estimator. A review of different frequency trackers together with their derivation and implementation details is given in [31].

_{k}#### 3.4. Correction of the variable offset frequency

The final step in the phase correction is a complex multiplication to counteract the offset frequency fluctuations. The following formula multiplies the adaptively-resampled signal by the offset phase fluctuations in counterphase,

*f*

_{0}〉 term ensures that the corrected spectrum remains at a constant noise-free position, which would otherwise be unnecessarily shifted by the correction procedure.

#### 3.5. Flowchart of the CoCoA algorithm

To illustrate how the THz DCS spectrum introduced in Fig. 2 is corrected, a flowchart of the signals after each step are given in Fig. 5. Fig. 5(a) shows the uncorrected IQ-demodulated photodetector signal after applying the FFT. As can be seen, the beat notes are not resolvable for the acquisition time of 1 ms. The results of the DDFG procedure, described in Section 3.1, and the Appendix, are shown in Fig. 5(b). The 8^{th} harmonic (white arrow) is chosen for demodulation and the retrieved instantaneous differential repetition rate is displayed in Fig. 5(c). The repetition rate correction is applied according to Eq. (3) and the resulting spectrum is shown in Fig. 5(e). If the resampling is effective, the resampled signal should demonstrate an increase in signal-to-noise for the Δ*f*_{rep} harmonics after the DDFG procedure. This is shown in Fig. 5(d), which in comparison to Fig. 5(b) has an average amplitude increase of ∼20 dB and clear reduction of the spectral width of the higher harmonics.

The frequency tracking part of the algorithm is applied on the adaptively-sampled signal of Fig. 5(e), after which the instantaneous difference frequency offset is retrieved. The complex multiplication of Eq. (7) finally gives the corrected rf spectrum shown in Fig. 5(g). Fig. 5(h) shows the effect of applying only the frequency offset correction, which clearly corrects for the majority of the phase noise.

## 4. CoCoA as a diagnostic tool for dual-comb coherence detection

As mentioned earlier, the CoCoA algorithm can be used as a diagnostic tool to verify the mutual coherence of the two sources contributing to the interferograms. To illustrate this property, two interband cascade lasers [21] centered at ∼3.5 μm were used to generate DCS interferograms while operating in either the comb-regime or high phase-noise regime. Fig. 6(a) shows the rf spectrum acquired when operating one of the lasers in the high-phase noise regime, indicated by the broad linewidth of the intermode beat note shown in Fig. 6(d). The other laser showed a sub-kilohertz linewidth of the intermode beat note [see Fig. 6(c)], which indicates low-noise comb operation of this laser. Similar to Fig. 5(a), the rf spectrum of Fig. 6(a) shows no resolvable beat notes, but as discussed earlier this could be the result of a long acquisition time. However, unlike the previously presented data, the self-mixing spectrum shown in Fig. 6(b) lacks the characteristic Δ*f*_{rep}-harmonics seen whenever sufficient mutual coherence is present. Consequently, the CoCoA algorithm does not provide any improvement in the spectrum as the requirement on the input signal to be coherent on a 1/Δ*f*_{rep} timescale is not satisfied. In fact, the beat notes will become even more dispersed as there is no correlation between the frequency components of the signal. In contrast, when operating both lasers in the comb-regime, the self-mixing spectrum of Fig. 7(b) is obtained. The spectrum exhibits multiple Δ*f*_{rep}-harmonics, any of which can be IQ-demodulated for the adaptive sampling step described in Subsection 3.2. The result of the CoCoA algorithm is apparent when comparing Fig. 7(a) and (c). Again, the algorithm provides linewidths that are close to the Fourier limit imposed by the acquisition time [shown in Fig. 7(e) and (f)]. The diagnostics provided by the DDFG self-mixing spectrum can be a useful tool especially when access to the intermode beat note, either through optical detection or rf extraction, is not available.

## 5. Comparison to an analog phase-locked loop stabilization

It is interesting to note the similarities of the CoCoA algorithm to that of an optical phase-locked loop (OPLL), which actively stabilizes the frequency offset between the sources, thereby significantly reducing the difference in frequency offset. The effects of such a procedure implemented on a pair of long-wave infrared (LWIR) QCL combs [32] is shown in Fig. 8 and more details on the implementation of the phase-locking is given in ref. [6].

As seen in the figure, the OPLL produces resolvable beat notes especially close to the region where the phase-reference local oscillator is placed. However, there is no adaptive sampling of the interferograms and hence an increase in linewidth can be observed for higher rf frequency components.

The CoCoA algorithm clearly outperforms the analog OPLL across the spectrum and provides an improved beat note signal-to-noise ratio without introducing any of the complex instrumentation that is associated with a high-bandwidth locking system. The benefits of phase-locking can thereby be obtained even for portable systems, where requirements on size, weight and power are imposed.

## 6. Compatibility with other dual-comb platforms

To verify the source-independent character of CoCoA, we applied the algorithm to a stream of interferograms recorded in two additional, distinctly different, dual-comb platforms: one based on silicon microresonators [33], and the other based on passively mode-locked waveguide lasers [34].

#### 6.1. Silicon microresonators

For the microresonator-based dual-comb system, our algorithm was applied to the periodic time-domain signal of Fig. 9(a), which is obtained from the multiheterodyne beating between two mid-infrared (∼3 μm) microresonator combs with a mode spacing of 127 GHz and a repetition rate mismatch of ∼41.6 MHz [33]; data-set provided by courtesy of Alexander Gaeta at Columbia University. The uncorrected spectrum of Fig. 9(b) clearly shows a multitude of discrete peaks. However, a zoom of the beat notes located at higher frequencies [Fig. 9(d) and (e)] reveals the existence of a MHz-wide noise pedestal. As in the QCL and ICL case, the algorithm enables Fourier-limited rf beat notes linewidths with an effective cancellation of the noise pedestal [see Fig. 9(c) and (f)], albeit with a correction improvement not as impressive as before due to the large separation between the beat notes and relatively high initial stability of the rf comb lines. Also visible in panel (c) is a strong suppression of weak spurious interfering signals located around 500 MHz. Because these spurious components do not match the comb source assumption of the algorithm, they will be dispersed by phase excursions of the comb lines, and thus weakened after correction. Consequently, CoCoA minimizes the contamination of the rf spectrum and simplifies peak tracking procedures for further spectroscopic analysis.

For spectroscopic applications, particular attention needs to be focused on how phase manipulation, the underlying core of CoCoA, affects the amplitudes of the beat notes. For this purpose, we performed an analysis of peak amplitudes obtained using the Fast Fourier Transform at different integration times for two beat notes of different strengths, located at 415 MHz, and 1039 MHz respectively. The results of the analysis are plotted in Fig. 9(g), (h) and reveal that the algorithm is capable of reaching much lower uncertainty of the noise-free amplitude estimate for integration times exceeding 10 μs, which corresponds to approximately *n* =410 interferograms. For the strong beat note, a single per mille regime of precision becomes accessible, as opposed to the uncorrected case which considerably drifts with a global minimum in the percentage range. A similar behavior can be observed for the weak beat notes, although accompanied by a 10-fold decrease in the achievable precision due to the ∼20 dB difference in rf power. Unfortunately, in the single microsecond averaging time regime (averaging tens of interferograms), we observe a slight increase in uncertainty of the amplitude estimate. This is likely an effect of noise propagation introduced by the phase retrieval and correction procedure. This clearly implies that the algorithm is primarily beneficial at longer time scale acquisitions where the reduction of the frequency drifts is significantly larger than the increase in the random noise component.

#### 6.2. Passively mode-locked waveguide lasers

CoCoA was also tested on a dual-comb data-set obtained from a system based on shared-cavity passively-mode locked waveguide lasers (WGLs) operating in the near-infrared telecommunication range of 1.55 μm [34] provided by courtesy of Jérôme Genest of Université Laval. The DCS signal was recorded with a fast oscilloscope sampling dual-comb interferograms with a repetition rate Δ*f*_{rep}=57 kHz, which mapped the optical domain with a spectral point spacing *f*_{rep}=822.4 MHz. The span of the rf dual-comb spectrum corresponds to approximately 2.3 THz in the optical domain, hence the number of rf lines was approximately *n* =2800, which is significantly more than typically achievable in ICL and QCL-based systems.

In contrast to the previously described sources, the WGL-based DCS signal possesses clearly identifiable bursts caused by their pulsed, true mode-locked operation [Fig. 10(a)]. Although the temporal structure is significantly different, correction was possible without any modification to the algorithm. Due to the relatively high initial stability of the WGL-based spectroscopic system, the gain of the correction is only on the order of 5 dB. Nevertheless, we observe a significant improvement of the rf line widths and despite the large number of comb teeth in the rf spectrum, the correction efficacy does not considerably decrease at frequencies located far away from the center (carrier) frequency, as shown in Fig. 10(c) and (d), and the computational time does not suffer from the increase in the number of modes. This demonstrates the compatibility of CoCoA with DCS based on mode-locked sources, although further investigation may be needed to facilitate a more effective retrieval of repetition rate fluctuations for even larger number of modes.

## 7. Summary

A purely computational coherent averaging (CoCoA) algorithm for free-running dual-comb spectroscopy is presented. The algorithm is compatible with modern real-time architectures and can be implemented without any additional electro-optical components, which simplifies the opto-mechanics of the system while keeping the costs low. The CoCoA algorithm works directly on the acquired digitized photodetector signal and its efficacy is here demonstrated on spectra acquired with two terahertz quantum cascade laser frequency combs [16, 36] operated under conditions that introduce significant amounts of phase noise. The generality of the proposed algorithm is verified by application to other systems with varying degrees of phase noise: mid-intrared quantum cascade and interband cascade laser frequency combs, mid-infrared silicon microresonators, and near-infrared waveguide mode-locked lasers. In principle, the procedure can be implemented in any DCS system, but the main benefits will be seen in free-running on-chip semiconductor systems that exhibit intrinsically high degrees of phase fluctuations. If implemented, coherent averaging over long time-scales is feasible uncovering rf components otherwise hidden in noise.

## Appendix Derivation of the digital difference frequency generation (DDFG) formula

To describe how this procedure works mathematically, we will start with the resampling step first, which employs Digital Difference Frequency Generation (DDFG). Consider a radio-frequency *N*-component complex multiheterodyne signal in the form:

*f*

_{0}is the rf offset frequency, and Δ

*f*

_{rep}is the rf repetition frequency. For simplicity, we assume equal amplitudes of the

*N*frequency components. Our goal is to retrieve

*N*harmonics of Δ

*f*

_{rep}, while eliminating any frequency components containing Δ

*f*

_{0}. This problem closely resembles difference frequency generation (DFG) – a process well known from nonlinear optics [37]. In the optical domain, the second order nonlinear susceptibility χ

^{(2)}of a medium is responsible for multiple frequency mixing phenomena, such as sum frequency generation (SFG), second harmonic generation (SHG), and difference frequency generation (DFG), all dependent upon the squared electric field. In principle, this scheme can be easily applied here to obtain the difference frequencies just by squaring the input signal. However, for a complex input we will just double and sum the signal’s frequencies without creating any difference components. Using the well-known formula for the square of the sum of

*N*numbers, we can write:

Let us now analyze, what happens if we process the real and imaginary components separately. By squaring the real part of *s̃*(*t*), we obtain the same components as given by *χ*^{(2)} in nonlinear optics:

*x*)cos(

*y*) = cos(

*x*+

*y*) + cos(

*x*−

*y*), where

*x*= 2

*π*Δ

*f*

_{0}, and

*y*= 2

*πn*Δ

*f*

_{rep}. Since the complex character of the signal is lost, the frequency spectrum becomes symmetric with a significant DC term. The latter in nonlinear optics is known as optical rectification, and here the same process takes place: by squaring the signal, we rectify it thus producing DC. More importantly, the DFG term now appears in the equation, yielding harmonics of Δ

*f*

_{rep}. Unfortunately, they will often be overlaid by spectrally-broad low frequency SHG and SFG components containing Δ

*f*

_{0}. To make matters worse, higher frequency SFG and SHG products can spectrally alias if they appear at frequencies greater than half of the sampling frequency, thereby spectrally overlaying the DFG products as well. Both effects will introduce a cross-talk between the DFG products used to retrieve correction signals, and the unwanted SHG and SFG by-products. Clearly, a nonlinear operation that yields the DFG components alone without any second harmonic or sum products including Δ

*f*

_{0}would be preferable.

If we now square the imaginary part of *s̃*(*t*), the result is nearly the same, however the sign of the components containing Δ*f*_{0} is changed:

*f*

_{0}disappears, therefore we arrive at the final formula for multiple harmonics of the repetition rate of a multiheterodyne signal, referred here to as Digital Difference Frequency Generation (DDFG):

## Funding

DARPA SCOUT Program (W31P4Q161001); DARPA SIGMA+ Program (HR00111920006); Thorlabs Inc.; The Kosciuszko Foundation Research Grant; Foundation for Polish Science (FNP) (START 085.2018).

## Acknowledgments

The authors would like to thank Prof. Jérôme Faist at ETH and Dr. Mahmood Bagheri at JPL for providing the LWIR QCLs and the ICLs used in this work. Prof. Qing Hu at MIT is acknowledged for providing the THz QCL combs. Furthermore, the authors express sincere gratitude to Prof. Alexander Gaeta and Dr. Mengjie Yu at Columbia University for the microresonator DCS data, and Prof. Jérôme Genest and Nicolas Hébert at Université Laval in Canada for providing the passively mode-locked waveguide laser DCS data. The authors would also like to thank Dr. David Burghoff at University of Notre Dame for fruitful discussions and comments during the developments of the CoCoA algorithm.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **D. van der Weide and F. Keilmann, “Coherent periodically pulsed radiation spectrometer,” US Patent 5,748,309 (May 5, 1998).

**2. **S. Schiller, “Spectrometry with frequency combs,” Opt. letters **27**, 766–768 (2002). [CrossRef]

**3. **I. Coddington, N. Newbury, and W. Swann, “Dual-comb spectroscopy,” Optica **3**, 414 (2016). [CrossRef]

**4. **I. Coddington, W. C. Swann, and N. R. Newbury, “Coherent Multiheterodyne Spectroscopy Using Stabilized Optical Frequency Combs,” Phys. Rev. Lett. **100**, 013902 (2008). [CrossRef] [PubMed]

**5. **G. Villares, A. Hugi, S. Blaser, and J. Faist, “Dual-comb spectroscopy based on quantum-cascade-laser frequency combs,” Nat. Commun. **5**, 5192 (2014). [CrossRef] [PubMed]

**6. **J. Westberg, L. A. Sterczewski, and G. Wysocki, “Mid-infrared multiheterodyne spectroscopy with phase-locked quantum cascade lasers,” Appl. Phys. Lett. **110**, 141108 (2017). [CrossRef]

**7. **T. Ideguchi, A. Poisson, G. Guelachvili, N. Picqué, and T. W. Hänsch, “Adaptive real-time dual-comb spectroscopy,” Nat. Commun. **5**, 3375 (2014). [CrossRef] [PubMed]

**8. **R. G. Wiley, *ELINT: The Interception and Analysis of Radar Signals* (Artech House, 2006).

**9. **P. Jouy, J. M. Wolf, Y. Bidaux, P. Allmendinger, M. Mangold, M. Beck, and J. Faist, “Dual comb operation of λ ∼ 8.2 μm quantum cascade laser frequency comb with 1 W optical power,” Appl. Phys. Lett. **111**, 141102 (2017). [CrossRef]

**10. **O. Kara, L. Maidment, T. Gardiner, P. G. Schunemann, and D. T. Reid, “Dual-comb spectroscopy in the spectral fingerprint region using OPGaP optical parametric oscillators,” Opt. Express **25**, 32713 (2017). [CrossRef]

**11. **A. V. Muraviev, V. O. Smolski, Z. E. Loparo, and K. L. Vodopyanov, “Massively parallel sensing of trace molecules and their isotopologues with broadband subharmonic mid-infrared frequency combs,” Nat. Photonics **12**, 209–214 (2018). [CrossRef]

**12. **J.-D. Deschênes and J. Genest, “Frequency-noise removal and on-line calibration for accurate frequency comb interference spectroscopy of acetylene,” Appl. Opt. **53**, 731–735 (2014). [CrossRef] [PubMed]

**13. **J. Roy, J.-D. Deschênes, S. Potvin, and J. Genest, “Continuous real-time correction and averaging for frequency comb interferometry,” Opt. Express **20**, 21932 (2012). [CrossRef] [PubMed]

**14. **Z. Zhu, K. Ni, Q. Zhou, and G. Wu, “Digital correction method for realizing a phase-stable dual-comb interferometer,” Opt. Express **26**, 16813–16823 (2018). [CrossRef] [PubMed]

**15. **A. Hugi, G. Villares, S. Blaser, H. C. Liu, and J. Faist, “Mid-infrared frequency comb based on a quantum cascade laser,” Nature **492**, 229–233 (2012). [CrossRef] [PubMed]

**16. **D. Burghoff, T.-Y. Kao, N. Han, C. W. I. Chan, X. Cai, Y. Yang, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, “Terahertz laser frequency combs,” Nat. Photonics **8**, 462–467 (2014). [CrossRef]

**17. **Q. Lu, D. Wu, S. Slivken, and M. Razeghi, “High efficiency quantum cascade laser frequency comb,” Sci. Reports **7**, 43806 (2017). [CrossRef]

**18. **Y. Bidaux, I. Sergachev, W. Wuester, R. Maulini, T. Gresch, A. Bismuto, S. Blaser, A. Muller, and J. Faist, “Plasmon-enhanced waveguide for dispersion compensation in mid-infrared quantum cascade laser frequency combs,” Opt. Lett. **42**, 1604 (2017). [CrossRef] [PubMed]

**19. **Y. Yang, D. Burghoff, J. Reno, and Q. Hu, “Achieving comb formation over the entire lasing range of quantum cascade lasers,” Opt. Lett. **42**, 3888 (2017). [CrossRef] [PubMed]

**20. **L. A. Sterczewski, J. Westberg, L. Patrick, C. Soo Kim, M. Kim, C. L. Canedy, W. Bewley, C. D. Merritt, I. Vurgaftman, J. R. Meyer, and G. Wysocki, “Multiheterodyne spectroscopy using interband cascade lasers,” Opt. Eng. **75**, 011014(2017).

**21. **M. Bagheri, C. Frez, L. A. Sterczewski, I. Gruidin, M. Fradet, I. Vurgaftman, C. L. Canedy, W. W. Bewley, C. D. Merritt, C. S. Kim, M. Kim, and J. R. Meyer, “Passively mode-locked interband cascade optical frequency combs,” Sci. Reports **8**, 3322 (2018). [CrossRef]

**22. **L. A. Sterczewski, J. Westberg, M. Bagheri, C. Frez, I. Vurgaftman, C. L. Canedy, W. W. Bewley, C. D. Merritt, C. S. Kim, M. Kim, J. R. Meyer, and G. Wysocki, “Mid-infrared dual-comb spectroscopy with interband cascade lasers,” Opt. Lett. **44**, 2113 (2019). [CrossRef] [PubMed]

**23. **M. Rösch, G. Scalari, G. Villares, L. Bosco, M. Beck, and J. Faist, “On-chip, self-detected terahertz dual-comb source,” Appl. Phys. Lett. **108**, 171104 (2016). [CrossRef]

**24. **F. J. Theis and A. Meyer-Bäse, *Biomedical Signal Analysis: Contemporary Methods and Applications* (MIT Press, 2010).

**25. **D. Burghoff, Y. Yang, and Q. Hu, “Computational multiheterodyne spectroscopy,” Sci. Adv. **2**, e1601227 (2016). [CrossRef] [PubMed]

**26. **N. B. Hébert, J. Genest, J.-D. Deschênes, H. Bergeron, G. Y. Chen, C. Khurmi, and D. G. Lancaster, “Self-corrected chip-based dual-comb spectrometer,” Opt. Express **25**, 8168 (2017). [CrossRef] [PubMed]

**27. **L. A. Sterczewski, J. Westberg, and G. Wysocki, “Molecular dispersion spectroscopy based on Fabry–Perot quantum cascade lasers,” Opt. Lett. **42**, 243 (2017). [CrossRef] [PubMed]

**28. **V. F. Pisarenko, “The Retrieval of Harmonics from a Covariance Function,” Geophys. J. Royal Astron. Soc. **33**, 347–366 (2007). [CrossRef]

**29. **P. Tichavsky and P. Handel, “Recursive estimation of frequencies and frequency rates of multiple cisoids in noise,” Signal Process. **58**, 117–129 (1997). [CrossRef]

**30. **R. M. Aarts, “Low-complexity tracking and estimation of frequency and amplitude of sinusoids,” Digit. Signal Process. **14**, 372–378 (2004). [CrossRef]

**31. **L. A. Sterczewski, “Signal processing in terahertz and mid-infrared spectroscopy with frequency combs,” Ph.D. thesis, Wroclaw University of Science and Technology (2018).

**32. **J. Westberg, L. A. Sterczewski, F. Kapsalidis, Y. Bidaux, J. Wolf, M. Beck, J. Faist, and G. Wysocki, “Dual-comb spectroscopy using plasmon-enhanced waveguide dispersion compensated quantum cascade lasers,” Opt. Lett. **43**4522–4525 (2018). [CrossRef]

**33. **M. Yu, Y. Okawachi, A. G. Griffith, N. Picqué, M. Lipson, and A. L. Gaeta, “Silicon-chip-based mid-infrared dual-comb spectroscopy,” Nat. communications **9**, 1869 (2018). [CrossRef]

**34. **N. B. Hébert, D. G. Lancaster, V. Michaud-Belleau, G. Y. Chen, and J. Genest, “Highly coherent free-running dual-comb chip platform,” Opt. Lett. **43**, 1814–1817 (2018). [CrossRef] [PubMed]

**35. **P. Werle, R. Mücke, and F. Slemr, “The limits of signal averaging in atmospheric trace-gas monitoring by tunable diode-laser absorption spectroscopy (TDLAS),” Appl. Phys. B **57**, 131–139 (1993). [CrossRef]

**36. **Y. Yang, D. Burghoff, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, “Terahertz multiheterodyne spectroscopy using laser frequency combs,” Optica **3**, 499 (2016). [CrossRef]

**37. **C. Li, “Polarization Theory of Nonlinear Medium,” in *Nonlinear Optics*, (Springer, 2017), pp. 23–50. [CrossRef]

**38. **P. Tichavsky and P. Handel, “Recursive estimation of linearly or harmonically modulated frequencies of multiple cisoids in noise,” IEEE **3**, pp. 1925–1928 (1997). [CrossRef]

**39. **J. Van Zaen, “*Efficient schemes for adaptive frequency tracking and their relevance for EEG and ECG*,” Tech. rep., EPFL (2012).

**40. **L. Tan and J. Jiang, “Novel adaptive IIR filter for frequency estimation and tracking [DSP Tips Tricks],” IEEE Signal Process. Mag. **26**, 186–189 (2009). [CrossRef]