Operational Bandwidth Extension in Continuous-Wave Terahertz Systems by Digital Postprocessing

Increasing the signal-to-noise ratio in terahertz systems extends their operational bandwidth and improves their dynamic range. Traditionally, the measurement noise is reduced by averaging each measured sample over a certain integration time interval, which proportionally increases the total measurement time. In this article, we demonstrate a digital postprocessing technique to suppress the noise by as much as 20 dB with a fractional additional postprocessing time. The operational bandwidth of a commercially available continuous-wave terahertz system is extended by $\sim$30% and the data quality is comparable to a measurement taken with two orders of magnitude longer integration time. We further utilize digital postprocessing to recover spectral information at frequencies where the raw data are in the noise floor. We also discuss its applicability in material characterization techniques using a photonic terahertz vector network analyzer.


I. INTRODUCTION
F OR more than half a century, the terahertz spectral range (0.1-10 THz) has been thoroughly investigated because of its potential applications in the fields of spectroscopy, astronomy, nondestructive testing, etc. These innumerous applications are underpinned by commercially available, table-top, noncryogenic terahertz systems, often equipped with powerful electronic and optical sources, emitting up to milliwatts of terahertz power [1], [2], [3]. Combined with sensitive receivers, they reach a dynamic range (DNR) as high as 120-140 dB [4], [5], at the lower end of the terahertz spectrum. Continuous-wave (CW) photomixing systems feature bandwidths as high as 4 THz [6], with frequency resolution of the order of 1 MHz [7] with comparatively simple laser systems. A comprehensive list of wide-spread applications of terahertz CW systems can be found in the literature [8]. Typically, all CW systems suffer from roll-offs that decrease the dynamic range usually with a power law of frequency (e.g., f 4 ) [9], [10]. Consequently, at higher terahertz frequencies the detected signal at the receiver becomes indistinguishable from the noise floor, mostly generated from the detection mechanism but also from the postdetection electronics. This limits the maximum operational frequency. The lock-in detection technique is typically employed to suppress noise by averaging the measured signal over the lock-in integration time, τ . This process, however, increases the measurement time proportionally in a homodyne system and even quadratically when direct detectors are employed. For example, a 1 THz scan in a typical commercial homodyne CW system [7] with a 50 MHz frequency step and 3 s integration time takes around a minute. To conduct scans above 2 THz, an integration time of at least 300 ms is necessary to extract meaningful data and a scan of similar frequency coverage will take nearly 3 h. The lasers operating the commercial system show significant phase drifts in such long scans. In this article, we demonstrate that noise can be significantly reduced by postprocessing with digital bandpass filters (BPFs) [11], yielding similar data quality to the unfiltered measurements that take 100 times longer. Causal digital filters are usually realized by weighted averaging of measured data, and hence, a loss of spectral resolution is unavoidable for band-limited signals, such as spectral peaks. Nevertheless, the filtering techniques, as demonstrated in this article, may further enable high-speed THz systems that fulfill industrial needs, in combination with faster hardware [6]. We further demonstrate the applications of digital postprocessing in CW terahertz spectroscopy, conducted around water vapor absorption resonances between 1.6 and 2.5 THz and a free-space photonic vector network analyzer (VNA) developed in house [12]. Fig. 1 shows a typical CW homodyne photomixing setup. The transmitter Tx is driven by a pair of 1550-nm lasers that differ in frequency by the desired THz frequency. In order to enable lock-in detection, Tx is further modulated at an intermediate frequency f IF , e.g., by a bias modulation. The THz signal propagates through free space, usually guided by optical components, such as lenses and parabolic mirrors. The receiver Rx detects the THz signal and generates a current I Rx oscillating at f IF . The detector current is proportional to the incident terahertz field on the receiver [9]. I Rx is typically of the order of a few nanoamperes down to a few picoamperes and is, therefore, amplified by a low-noise transimpedance amplifier (TIA) before forwarding This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Fig. 1. Schematic of a typical configuration of a terahertz homodyne, CW system. In this work, we use a p-i-n-diode-based Tx from Toptica Photonics/Frauhofer Heinrich Hertz Institute and an ErAs:In(Al)GaAs-based Rx, developed in-house [13]. Both the Tx and Rx are antenna coupled and are attached to hyper hemispherical silicon lenses, facilitating their free-space operation. The source control module contains the LIA, whose ENBW is inversely proportional to τ .
it to a single channel lock-in amplifier (LIA) with a fixed phase, optimized for the electrical path length of the cabling, for lock-in detection. The Rx is driven by the set of same lasers that drive Tx ensuring perfect phase locking. The generated detector current undergoes an additional frequency-dependent (optical) phase change induced by the phase difference between the THz wave and the laser signal driving the receiver, corresponding to a path length difference of Δ(nd). The detector output current I det obtained at a THz frequency f reads [14] I where c 0 is the speed of light in vacuum and I N (f, f IF ) is an additional noise current composed of IF noise sources and noise/clutter from the THz path. The IF noise consists of mostly white noise originating from the resistance of the device, the TIA, and the LIA. It can be described by an equivalent noise resistance R N with a voltage noise of where T is the temperature, k B the Boltzmann constant, and ENBW = b/τ is the equivalent noise bandwidth, wherein τ is the lock-in time constant, referred henceforth as (lock-in) integration time. The factor b is a constant that depends on the roll-off of the low-pass filter of the lock in and is of the order of 0.5. As the detected current is proportional to the THz field, the THz noise power is proportional to the square of the noise current and, thus, scales linearly with the equivalent noise bandwidth (ENBW), i.e., ∝ 1/τ . Simply by increasing the lock-in integration time τ from 3 to 300 ms, the DNR improves by 20 dB, as shown for the measurements in Fig. 2(a). Generally speaking, since integration time is proportional to the total measurement time, each 10 dB DNR improvement require an order of magnitude increment in measurement time. The hard limit of noise reduction is given by systematic errors at extremely long measurement times. For integrated terahertz applications, such as optoelectronic VNA [15], which require both high system dynamic range and low measurement time, decreasing measured (a) Measured DNR of a commercial CW system with integration times 3, 30, and 300 ms, respectively, between 1.5 and 2.7 THz. A tenfold increment in integration time increases the system DNR by 10 dB. The CW system can be operated up to 1.8 and 2.4 THz with τ = 3 and 300 ms, respectively. (b) Unfiltered detector current is filtered using the equivalent window function generated by the IIR filter. Unwanted correlated reflections from the components in the measurement setup come both before and after the main signal. noise just by increasing τ becomes impractical. Noise and clutter sources originating from the reflections and stray radiation in the setup are usually colored, and hence do not scale inversely with integration time. In the following, we will show that appropriate filtering techniques not only allow for removing clutter but also drastically reduce the noise floor without increasing the measurement time.

II. THEORETICAL DESCRIPTION
The sinusoidal oscillations, induced by the phase difference in the homodyne setup, act as a carrier signal for the spectral information stored in its envelope. The envelope of the detected signal can be calculated using the Hilbert transformation [14], [16]. In the time domain, the envelope information translates to a pulse, occurring at a delay dictated by the homodyne oscillations, such as the path length difference Δ(nd). Reflections in the setup occur at different path lengths. These reflections can simply be filtered out by windowing in the time domain, corresponding to bandpass filtering of the frequency scan. The filtering operation is implemented using the "filtfilt" function in MATLAB [17], which is an implementation of anticausal zero-phase filtering, i.e., phase distortion does not occur during the filtering process [18]. Thus, both finite impulse response (FIR) and infinite impulse response (IIR) filters preserve the phase information of the signal. IIR filters always have lower order than its FIR variants, and thus, it is the filter of choice in this work. The designed filter bandwidth and the center frequency is dynamically calculated from the measured signal and a 40 dB stop-band attenuation is implemented. Fig. 2(b) shows a bandpass filtering with an exemplary normalized window size of 1% around the main peak A of the detector current from the LIA (I LIA ). It suppresses noise and clutter in the frequency domain by ≈20 dB. This improves the DNR by the same amount, and, as a result, increases the maximum operational frequency by a factor of ∼1.3. The BPF further removes the colored noise and clutter components, i.e., the unwanted standing waves and stray reflections, R1 and R2, occurring at various optical components in the table-top terahertz setup. The overseeable loss of resolution for bandlimited signals by the applied filtering techniques are usually not dramatic: homodyne spectrometers use frequency step sizes that must be small enough to resolve the homodyne fringes, which are usually chosen to have a period much smaller than a resonance under investigation. Typical step sizes are in the range of a few MHz to a few 10 MHz. The aforementioned filtering example by a factor of 100 in the time domain performs an equivalent weighted averaging of the frequency sweep data over ≈180 samples, 1 i.e., over a frequency span of ∼1-2 GHz. Any spectral feature narrower than this resolution will be broadened. Since the filter function and its Fourier transform are purely analytical, i.e., do not contain any statistic fluctuations, the broadening can be undone (theoretically) by a deconvolution. Regardless, the resulting broadened resolution is still almost an order of magnitude better than that of typical time domain systems.

III. RESULTS AND DISCUSSION
We use a commercial "Terascan" system from Toptica Photonics AG [7], with a p-i-n diode transmitter and a photoconductive receiver. We conducted frequency scans between 1.5 and 2.7 THz with integration times ranging between 3 and 300 ms. The yellow plot in Fig. 3(a) shows the inverse Fourier transformed detected current I LIA measured with an integration time of 3 ms, and Fig. 3(b) shows the corresponding envelope in yellow, which becomes indistinguishable from noise at ∼ 1.9 THz. The plots in red in both Fig. 3(a) and (b) show the bandpass-filtered signal, where the noise is significantly reduced and spectral features, such as water absorption lines become visible till ∼ 2.5 THz. For comparison, the water absorption data from the HITRAN database [19] are also added in Fig. 3(b). We note that above 2.5 THz, the signal power becomes comparable to the noise floor, and thus, water absorption spectra cannot be accurately distinguished from noise artefacts. BPFs used in postprocessing only reduce the out-of-band noise in the received signal, unlike the averaging process in the LIA, which reduces the total measured noise power (both in-band and out-of-band noise) by reducing the ENBW.  Fig. 4(a) compares the DNR of frequency sweeps conducted with τ = 3 ms (M1) and τ = 300 ms (M2) to a filtered version of M1. An IIR BPF of order ten is used in this case. The figure shows that the filtered data features ≈20 dB higher dynamic range and is comparable to M2, i.e., the data obtained using 100-times larger integration time, and consequently, 100-times longer measurement time. The quality factors (Q-factors) of the resonance peaks have not noticeably decreased. Fig. 4(b) shows a systematic reduction of measured noise current after filtering for datasets measured with seven different integration times ranging between 3 and 300 ms. The employed filter reduces the overall noise current by an order of magnitude, equivalent to a 20 dB improvement in DNR. However, we expect this DNR improvement will eventually saturate for very high integration times when the noise power in the unfiltered dataset is extremely low.
The additional postprocessing times required for the filter implementation and filtering process of a dataset of ≈41600 samples is of the order of tenths of a second with MATLAB running on a 64-bit windows 10 PC with 8 gigabytes of RAM and a 4-Core Intel Core i7-7700 CPU. Executing these filtering algorithms in programming languages, such as "C," should be faster. Hardware implementations (e.g., FPGAs) can further improve performance, however, at a cost of adaptability and flexibility.

A. Application in Spectroscopy
Spectroscopic measurements in terahertz domain involves the study of absorption spectra generated by various materials under investigation, such as polar gases, polymers, and crystalline materials [20], [21]. Traditionally, time-domain terahertz systems are employed for spectroscopy due to their broadband characterization capabilities at high-speed. However, these systems usually feature a spectral resolution in the order of a few gigahertz, limited by physical size of the employed delay stage. In contrast, CW frequency-domain Terahertz systems feature spectral resolutions as small as 1 MHz [7]. Furthermore, homodyne systems using Terahertz field detectors give access to the phase information of the spectral features, which, even though is frequently ignored in spectroscopy, can provide further information about the Q-factor of detected resonances [16], [22]. Datasets containing peaks (poles) can be better modeled as an autoregressive process, and hence, IIR filters better model their impulse response than their FIR counterpart. Further, spectroscopic fingerprints at room temperature and normal pressure typically have resonance peaks with widths in the few GHz range, i.e., 2-3 orders of magnitude larger than the resolution of CW photomixing systems. Thus, broadening of the resonance peaks caused by filtering is not dramatic. In addition, the phase information is often not used in spectroscopy, hence even nonzero-phase filtering using IIR filters can produce accurate results.
We investigate water absorption lines as a representative spectroscopic dataset, measured using the homodyne setup shown in Fig. 1 and compare the resonance peaks of an unfiltered dataset measured with τ = 300 ms (yellow) to a filtered dataset measured with τ = 3 ms (red). The total path length between the Tx and the Rx is 28 cm, where the terahertz beam propagates through ambient air. An IIR filter of order 10 and 1% normalized bandwidth is employed. Fig. 5(a)-(d) shows four different waterlines at 1.671, 1.797, 2.042, and 2.345 THz, respectively. To calculate the Q-factor of the resonance in the unfiltered data, we fit it to a Lorentzian function [23]. The filtered data almost perfectly replicate the fitted resonance features till at least ∼ 2 THz. In addition, the filtering technique allows to visualize resonances where even the 300 ms measurement fails as the resonances are hardly visible and hidden in the noise floor [cf., Fig. 5(d)]. An estimated envelope spectrum is calculated from the HITRAN database [19] for a laboratory temperature of 18 • C and a relative humidity of 30% at standard atmospheric pressure  Table I lists the measured Q-factors for unfiltered and filtered versions of measurements taken with τ = 3, 30, and 300 ms, respectively. The Q-factors of the filtered dataset with τ = 3 ms is comparable or higher than the unfiltered dataset measured with τ = 300 ms as they are less obscured by noise. As compared with the raw data, filtering increases the Q-factor and mostly gets it closer to the value of the HITRAN database, particularly where the raw data peaks are not well resolved as they are too close to the noise floor. Note that the Q-factors of filtered dataset measured with τ = 3 ms is usually similar or larger than those of the raw data measured with τ = 300 ms. This emphasizes that the loss of spectral resolution due to the employed BPF has minuscule effect on the shape of the resonances. The relationship between spectral broadening of the resonance features and DNRimprovement is detailed in Appendix B.
We further observe a slight shift in resonance frequency of about 1 GHz (i.e., ∼ 16% of the width of the resonance) upon filtering as compared with the fitted raw data. It probably occurs due to the noise-induced inaccuracies in the fitting algorithm of the unfiltered data. A further shift of ≈1.3 GHz, caused by a lag between set and actual frequency, has been calibrated out.
In order to further demonstrate the strength of the filtering technique, Fig. 6 compares two additional plots. First, a filtered dataset measured with τ = 300 ms and corrected for the 1.3 GHz frequency offset, and second, the estimated spectrum from the HITRAN database. The filtered THz spectrum measured with 300 ms matches excellently to the estimated ideal spectrum I HT , and we can observe waterlines up to the maximum measurable frequency of the Toptica Terascan system [7] of 2.75 THz. Insets show a zoom-in to the frequency ranges 1.65 − 1.75 THz and 2.3 − 2.5 THz. The noise floor of the filtered data is at ≈0.4 pA, i.e., about an order of magnitude lower than that of the thermal noise floor of the photoconductor.

B. Photonic Vector Network Analyzer
CW homodyne systems, with two transmitters and two receivers, can be utilized to setup a free-space optoelectronic VNA [12]. Apart from circuit characterization, VNAs are frequently used for high-precision spectroscopy and material characterization, i.e., extracting thickness, refractive index, and absorption coefficients of the material under test. Imaging of nanometric depositions on a highly resistive silicon wafer using a photonic system has also been demonstrated [14], enabled by accurate phase resolution of the transmission coefficient. VNAs are, in general, used to obtain complex scattering parameters of the devices under test, strictly requiring preservation of the phase information. IIR filters are no linear phase filters; however, implementation of the zero-phase filtering technique helps to preserve the phase characteristics of the actual signal.
We demonstrate an application example using a photonic terahertz VNA for characterization of a ∼ 220 μm thick quartz plate. Plots in Fig. 7(a) shows the measured field transmission and reflection coefficients. The filtered transmission and reflection coefficients excellently match the theoretical predictions. The plot in Fig. 2(b) shows the equivalent time-domain pulse of the unfiltered reflected current at the VNA port. The observed stray reflections R1 originate from the two polarizers employed in the setup [12], which appear before the reflected signal from the sample under test (e.g., quartz plate). Filtering can remove these unwanted, systematic sources of clutter. Reflections R2 are the remnants of multiple reflections occurring in the cavity between the sample and the polarizers, which appear after the reflected signal from the sample. These stray reflections are usually correlated to the signal components and constitute the colored noise component in the measured data, which cannot be suppressed by increasing the integration time. Bandpass filtering attenuates these reflections and, thus, is crucial for reflection coefficient measurement in photonic VNAs. Transmission measurements would only have stray reflections temporally trailing the main signal, and hence a low-pass filter would suffice. The residual noise power after bandpass filtering is proportional to the filter bandwidth. In other words, narrower filter bandwidths increase the system DNR. However, the appropriate filter bandwidth should be determined considering the type of application, frequency resolution, and scanning bandwidth. On the one hand, narrow bandwidths are crucial for spectroscopic measurements to remove any additional clutter arising in the measurement setup. On the other hand, too narrow bandwidth will suppress multiple reflections from samples under test (e.g., quartz plate and distributed Bragg reflectors), thereby not resolving the Fabry-Pérot interference resulting from the resonant cavity of the structures. Broader passband widths increase the filter order of IIR filters significantly and, thus, are more suited for spectroscopic measurements. The filter order of FIR filters depend on the transition width between the pass and stop bands and the stopband attenuation, and thus are better suited for transmission and reflection measurements in a VNA, where correlated reflections arising from the material under test are nontrivial and, thus, require broader pass-and transition bands.
The commonly employed noise reduction techniques by averaging over N samples is a discrete implementation of the integration-based low-pass filtering and is equivalent to a lowpass FIR filter, where all filter coefficients are equal. The resulting low-pass filter is, thus, a sinc function, with high passband ripple and the main sidelobe as high as 4% of the main peak, i.e., ≈26.5 dB suppression. Alternatively, a low-pass filter of the same order, but with varying filter coefficients, can be implemented to have a very small passband ripple and a sidelobe suppression as high as 40-60 dB, and consequently have a superior performance than averaging filters [24]. BPFs for homodyne setups reduce the out-of-band noise even further. Low-pass filtering with subsequent downsampling can also significantly reduce measured noise power, however at the cost of reduced spectral resolution. Wavelet-based filtering techniques notably improve the system DNR, but spectral features are lost when the signal power of the unfiltered data becomes comparable to the noise power, and thus are incapable of extending the operational bandwidth of the system. Alternatively, different window functions can be multiplied to the inverse Fourier-transformed time-domain data and, subsequently, can be Fourier transformed to obtain the frequency spectrum. A comprehensive study of effects of different windows shapes on time-domain spectroscopic data is done by Vázquez-Cabo et al. [25]. However, this process is more computationally expensive than directly using BPFs on the frequency-domain data.
We further note that the phase-zero filtering technique is important to remove phase distortion and delays induced by the filtering processing. For hardware implementation, however, standard filtering architecture (direct form I and II) are simpler, and in such cases, IIR filters should only be utilized where phase information is not critical.

IV. CONCLUSION
Digital filters in postprocessing can significantly improve system performance for homodyne CW terahertz systems. We have demonstrated a dataset with ∼ 41600 samples, measured over a time period of 2 min (additional ∼ 0.19 s of postprocessing time) with τ = 3 ms filtered to have comparable DNR and spectral coverage to a dataset measured with τ = 300 ms over a time period of 3 h. The operational bandwidth of the terahertz system at τ = 3 ms extends from 1.9 to, at least, 2.5 THz with filtering techniques, recovering spectral data even though the raw data seem to be dominated by noise. We have also shown application examples where digital filters can be used to significantly improve the data quality measured with a very low integration time. Measurement noise can be alternatively suppressed, even adaptively, using more complex Wiener filters [26]. The measurement dataset can be further improved using higher order autoregressive-moving-average models combining of FIR and IIR filters, or even using prediction filters to remove all unwanted noise components; however, these techniques require more complex computation and a-priori knowledge of the sample under test for proper modeling.

APPENDIX A INTEGRATION AT LIA
The mathematical equivalent of the (ideal) LIA output current, I LIA , integrated over an integration time τ is where I TIA is the input current to LIA after amplification with the TIA and r(t) = cos(ω mod t) is the normalized reference signal. While the received signal I det (t) as the major component of I TIA is also modulated at the reference frequency, it contains further frequency components originating from noise with an amplitude I N,ω . With I TIA,ω (t) = I det cos(ω mod t + ϕ ω ) + I N,ω cos(ωt + ϕ ω ), the integral can alternatively be written as As the integration time is chosen usually at least several times larger than the inverse modulation frequency, the second term essentially integrates to zero. The first term, however, turns into a constant for ω = ω mod , i.e., when the modulation frequency of the signal matches the reference frequency. If the signal frequency is slightly off with |ω − ω mod | −1 < τ, its oscillation is so slow that the integral returns a nonzero result, capturing also some noise components in the vicinity of ω mod . At larger difference frequencies, the integral decays to zero as sinc(ω − ω mod )τ . The sinc decay in the time domain turns into a boxcar decay in the frequency domain where Π(t) is the normalized boxcar function and I TIA (ω) is the frequency spectrum of the TIA current. The integration operation at LIA, thus, acts as a low-pass filter, with a bandwidth of π/τ [24], making the ENBW of the LIA inversely proportional to τ . In real LIAs, the filter function may differ from the ideal boxcar function with thus slightly different values for the ENBW, however, the same scaling with τ .

APPENDIX B DNR IMPROVEMENT AND SPECTRAL BROADENING
BPFs can effectively remove colored noise from the measured data; however, the residual white noise after lock-in detection cannot be completely eliminated. For the IIR filter, as shown in Fig. 2(b), the total unattenuated noise power (σ 2 N,filt ) passing through the filter can be written as where σ 2 N is the total noise power and ω p is the normalized 3-dB bandwidth of the filter. Thus, the estimated DNR improvement (ΔDNR) with such a filter can be approximately calculated to Δ DNR ≈ 10 log 10 σ 2 N σ 2 N,filt = −10 log 10 (ω p ).
Evidently, a smaller 3-dB bandwidth of the BPF ensures a larger DNR improvement, however, the choice of filter bandwidth depends on the type of application in question. In spectroscopic applications, the resonant features of materials under test are bandlimited, and consequently, bandpass filtering of such signals results in a loss of spectral resolution due to the weighted averaging of the spectral data resulting from convolution. This leads to broadening and smoothing of the spectral resonances, as observed in Fig. 6. Smaller ω p results in a higher order averaging, leading to a higher loss of spectral resolution. Numerically, we found that the 3-dB width of the envelope of the impulse response of the BPF should be at least three times smaller than the 3-dB bandwidth of the spectral resonance to avoid the broadening of such resonance features. Approximating the passband of the BPF to a boxcar function, a lower limit of ω p can be, thus, calculated to ω p ≥ 2.67 Δf reso Δf spec (7) where Δf reso is the resolution of the frequency scan and Δf spec is the 3-dB linewidth of the spectral resonance. For example, the 3-dB linewidth of waterline at 1.671 THz can be calculated from the HITRAN data in Table I to be 5.90 GHz. Thus, a maximum DNR improvement without broadening of this particular spectral feature, for a Δf reso = 50 MHz, can be calculated from (6) and (7) to be 16.45 dB. We, however, used a narrower filter to achieve ∼20 dB improvement in DNR in this work, which resulted in an evidently reduced Q-factor for the investigated spectral features. Applications, where Fabry-Pérot resonances are investigated, need an even wider passband to avoid attenuating the subsequent reflections originating at the resonant cavity.