Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Reducing errors in THz material parameter determination by model-based time-domain extraction methods

Open Access Open Access

Abstract

We demonstrate a reduction of THz material parameter estimation errors by a factor of 40, using a time-domain material-model-based extraction method, in comparison with the respective frequency-domain approach. Based on simulated and experimentally measured THz time-domain spectroscopy data with different levels of the signal-to-noise ratios (SNR), an improved robustness of the time-domain scheme with regard to noise is demonstrated. This approach thereby provides an additional tool for the extraction of material parameters, which is especially suited for the evaluation of spectral data of samples with strong resonant absorption features, even if the SNR at those features drops below 0 dB. The cause for the difference in performance between the frequency-domain and time-domain methods is identified and a formal description provided.

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

1. INTRODUCTION

In recent years, a large number of applications for THz systems have been demonstrated, ranging from 3D security and object imaging [1] to ultra-broadband device characterization [2], biosensing [3], and spectroscopic material analysis and recognition [4]. The advantages of THz radiation are its nonionizing nature, penetration capability through many dielectric materials (e.g.,  paper, plastics, cloth), and the presence of unique spectral characteristics of materials of interest, like explosives, drugs, pharmaceuticals, and chemicals [46]. It thereby provides a valuable tool for chemical compound analysis and process monitoring [7], sensing of gases [8], and security applications [4]. THz time-domain spectroscopy (THz-TDS) provides the most flexible access to this spectral information, as it covers a broad frequency range, exceeding those of electronic or tuneable continuous-wave laser systems [9]. In order to evaluate and compare the spectral characteristics of a substance, general material parameters, which are independent from the sample dimensions and the measurement system, have to be derived from the acquired measurement data, like the complex refractive index, the complex dielectric constant, and/or related material properties.

In the classical approach, these material parameters are determined individually for every frequency point. Such calculations are based on a frequency-space-based formulation of propagation models, i.e., the mathematical representations of the measured data in dependence of the sample material parameters. By comparing the measured spectral data to the modeled data, the value of the material parameters is extracted by numerical optimization, solving the posed inverse problem. This optimization approach is necessary, as the relationship between measured data and material data cannot be inverted analytically. Such calculations are then repeated for every frequency of the investigated spectrum in a point-by-point process, as shown in earlier work [10,11]. Based on these material parameters, further attributes of the sample, like relaxation-times, resonator strength, and resonant frequencies or plasma frequencies can be derived. Given sufficient a priori knowledge on a sample, like the general material class and number of spectral features, these material properties can be determined by assuming an analytic description for the physical material model (e.g.,  Lorentz oscillator model) to the previously extracted material parameters [6,12], which at the same time provides a condensed representation of those parameters over the frequency range. A schematic of such a standard frequency domain (FD) processing approach is displayed in the left section of Fig. 1.

 figure: Fig. 1.

Fig. 1. Schematics of the classical frequency-based approach (left, FD), evaluating each frequency point individually, and the material-model-based approaches in the frequency (center, FDm), and time domains (right, TDm), indicating the measured and determined data (rounded boxes) and the numerical processing steps.

Download Full Size | PDF

It has been demonstrated that the determination of material parameters is notably improved by merging these two processing steps (i.e., propagation model and physical material model) into one operation [13,14]. Since the noise-sensitive optimization of the individual frequency-dependent material parameters [15] is replaced by the optimization of the frequency-invariant global material properties of the physical models, the optimization process operates synchronously on the complete spectral data set with several hundred complex data points instead of on only two values (amplitude and phase) for each individual frequency point. This frequency-domain model-based analytic procedure (FDm) is schematically depicted in the central section of Fig. 1.

In our analytic procedure, we make use of the fact that, as the sample and its effect on the THz wave are described over the full spectral range by the parametric material model, the use of inverse Fourier transforms to derive time-domain model data becomes feasible, allowing the realization of a time-domain minimization and parameter extraction procedure [1618]. This provides an additional time-domain alternative for the determination of THz material parameters, which can be used as an extension of the frequency-domain-based processing. While inherent limitations of generally applicable frequency-domain methods can be circumvented, it is important to note that this time-domain approach can only be used for application cases in which some general a priori information on the material are available to select the appropriate parametric material model. This time-domain method (FDm) depicted in the right part of Fig. 1 is therefore particularly appealing for quality control, industrial analyses, or recognition tasks where the general analytic task is predefined.

In this work, we present a detailed, comparative study of the performance of material-model-based approaches operating in the time or frequency domains, used for the determination of THz material parameters from measured THz data. The methods are applied to simulated data sets, including artificial noise, allowing the generation of system-independent, fully controllable testing scenarios with designable data characteristics to allow precise quantification of the quality of the extracted material data. This investigation is complemented by the demonstration and analysis of parameters extracted from measured, experimental data. The capabilities of the methods in both testing scenarios are compared quantitatively, and the cause for differences in their performances is formally determined.

2. FORMAL DESCRIPTIONS

A. Propagation Model

To determine the complex refractive index of a sample from measured data, its effect on the THz radiation has to be formally expressed as the samples’ refractive-index-dependent relative amplitude transmission coefficient ${T_{\rm{samp}}}$ [10]. In the frequency domain, this can be described by the ratio between the field transmitted through the sample ${E_{\rm{samp}}}$ and the corresponding air reference ${E_{\rm{air}}}$ measured without a sample. These are equivalent to their respective absolute transfer functions ${H_{\rm{samp}}}$ and ${H_{\rm{air}}}$ multiplied by the impinging electric field ${E_{\rm{inp}}}$:

$${T_{\rm{samp}}}(f) = \frac{{{E_{\rm{samp}}}(f)}}{{{E_{\rm{air}}}(f)}} = \frac{{{H_{\rm{samp}}}(f) \cdot {E_{\rm{inp}}}(f)}}{{{H_{\rm{air}}}(f) \cdot {E_{\rm{inp}}}(f)}} = \frac{{{H_{\rm{samp}}}(f)}}{{{H_{\rm{air}}}(f)}}.$$

In the following, the transmission through an optically thick, single slab sample at normal incidence is considered. As the internal reflections in the sample can be excluded for an optically thick sample by temporal windowing, this allows the formulation of the most elementary propagation model. Following the description in [10], the transfer function is provided in this case by

$$\begin{split}{T_{\rm{samp\_calc}}}(f) &= \frac{{2{\eta _{\rm{air}}}(f)}}{{{\eta _{\rm{air}}}(f) + {\eta _{\rm{samp}}}(f)}} \cdot \exp\left({- i{\eta _{\rm{samp}}}(f)\frac{{2\pi\! f{l_{\rm{samp}}}}}{{{c_0}}}} \right)\\&\quad\cdot \frac{{2{\eta _{\rm{samp}}}(f)}}{{{\eta _{\rm{samp}}}(f) + {\eta _{\rm{air}}}(f)}} \cdot \exp{\left(\!{- i{\eta _{\rm{air}}}(f)\frac{{2\pi\! f{l_{\rm{samp}}}}}{{{c_0}}}} \!\right)^{- 1}}.\end{split}$$
Here, the frequency-dependent, complex refractive index of the sample
$${\eta _{\rm{samp}}}(f) = {n_{\rm{samp}}}(f) - i{\kappa _{\rm{samp}}}(f)$$
is composed of the refractive index ${n_{\rm{samp}}}$ and the extinction coefficient ${\kappa _{\rm{samp}}}$, the known complex refractive index ${\eta _{\rm{air}}}$ of the surrounding atmosphere and the samples thickness ${l_{\rm{samp}}}$.

With ${E_{\rm{samp\_meas}}}(f)$ and ${E_{\rm{air\_meas}}}(f)$ being the Fourier transforms of the measured transients, the experimental transfer function ${T_{\rm{samp\_meas}}}(f)$ can be calculated. The complex refractive index can then be derived from the measured data, by adapting the parameters of the calculated transmission in Eq. (2) to the measured data ${T_{\rm{samp\_meas}}}(f)$ and solving the resulting equation:

$${T_{\rm{samp\_meas}}}(f) = \frac{{{E_{\rm{samp\_meas}}}(f)}}{{{E_{\rm{air\_meas}}}(f)}} = {T_{\rm{samp\_calc}}}(f).$$
From Eq. (2), it is directly evident that no exact analytical solution is achievable and a numerical solution for the material parameter extraction becomes necessary [10], e.g.,  by minimizing the difference between the calculated and measured transmission magnitude
$$\Delta {T_R} = | {{T_{\rm{samp\_calc}}}(f)}| -| {{T_{\rm{samp\_meas}}}(f)} |$$
and phase
$$\Delta {T_{\rm{arg}}} = {\rm arg}({{T_{\rm{samp\_calc}}}(f)} ) - {\rm arg}({{T_{\rm{samp\_meas}}}(f)} ),$$
with ${\rm arg}(T)$ representing the unwrapped phase information. This unwrapping is required to remove the ambiguity of the discrete Fourier transform phase representation. This ambiguity would otherwise produce absolute errors for the real value of the complex refractive index and completely inhibit a correct parameter extraction.

A direct solution of Eq. (4) in the time domain in this form is not feasible, since this equation is defined for each frequency point individually by ${\eta _{\rm{samp}}}(f)$.

B. Material Models

As the information contained in THz-TDS data is not limited to the samples’ complex dielectric material values, more relevant attributes, like the properties of the underlying physical mechanisms that determine the samples interaction with the THz radiation, can also be derived. This is typically achieved by fitting a physical material model, i.e., a mathematical description of those mechanisms, to the complex refractive index data. Since the material models are reasonably based on the desired physical material properties, values like resonator strength, resonant frequencies, and relaxation times, for instance, are directly available as part of the fitting result. The selection of a suited material model is based on general a priori knowledge on the investigated material type, including, for instance, the overall number and type of spectral features.

As an example, a typical physical material model suited for many (nonmagnetic) electronically resonant materials in the THz range, like explosives, pharmaceutical substances, and chemicals, is the Lorentz oscillator model in Eq. (7) with multiple contributing oscillators [6,16]:

$$\begin{split}{\eta ^2}(f) &= \varepsilon (f) \\ &= {\varepsilon _\infty} + \sum\limits_a \frac{{{P_a}^2}}{{{{(2\pi \!{f_a})}^2} - {{(2\pi\! f)}^2} + i{\gamma _a}2\pi\! f}}.\end{split}$$
This physical material model describes the frequency-dependent dielectric constant $\varepsilon$ of a material, which is the square of its complex refractive index, based on the dielectric constant at the high-frequency limit ${\varepsilon _\infty}$ and the resonance frequency ${f_a}$, damping constant ${\gamma _a}$, and strength ${P_a}$ of the individual oscillator $a$.

Considering that the propagation model as well as the physical material model are dependent on the complex refractive index, a combination of both models can easily be realized in one step. By substituting the refractive index in Eq. (4) with Eq. (7). This step not only replaces all individual frequency-dependent unknowns by frequency-invariant global quantities but also massively reduces their number, since the typically hundreds of refractive index values at individual frequencies are now represented by the material model with a few physically relevant material quantities. Utilizing this combined model, the physical material properties of a sample can directly be determined in one optimization step as elucidated for the FDm and TDm case in Fig. 1, matching the models’ complex output to the measured transmission data simultaneously over the complete spectral range. The corresponding material parameters are calculated by the combined model and do not require an additional extraction step as depicted for the FD case.

Since the transfer function can now be described for the complete investigated spectral range by the frequency-invariant material properties, an inverse Fourier transform can easily be applied, and the optimization procedure can alternatively be performed in the time domain by comparing the convolution product of the acquired reference transient and the calculated transfer function with the acquired sample transient. This convolution product is determined by calculating the transfer function based on the previously described combined model in the frequency domain, multiplying the result with the Fourier transform of the reference transient and performing the inverse Fourier transformation. With

$$\begin{split}{E_{\rm{samp\_meas}}}(t) &= {T_{\rm{samp\_calc}}}(t)*{E_{\rm{ref\_meas}}}(t) \\& = {{\cal F}^{- 1}}({{T_{\rm{samp\_calc}}}(f) \cdot {E_{\rm{ref\_meas}}}(f)}),\end{split}$$
the optimization minimizes the difference
$$\Delta {E_{\rm{samp}}}(t) = {E_{\rm{samp\_meas}}}(t) - {{\cal F}^{- 1}}({{T_{\rm{samp\_calc}}}(f) \cdot {E_{\rm{ref\_meas}}}(f)}).$$

C. Noise Propagation and Material Parameter Error

In order to determine and compare the performance of the time-domain (TDm) and frequency-domain (FDm) extraction procedures, the impact of noise and resulting errors on the processed data is investigated. Since measurements are performed in the time domain, providing the transient $E({t_k})$ in its discrete form $E(k\tau)$, the standard deviation ${\sigma _E}(k\tau)$ can directly be derived as a measure of noise error in this domain, with the sampling interval $\tau$ and the time index $k$. For the frequency-domain case, the propagation of this standard deviation after Fourier transformation has to be addressed in order to calculate the corresponding standard deviation for the signal magnitude ${\sigma _{| E |}}(f)$ and phase ${\sigma _\vartheta}(f)$. Following the formalism derived in [19,20], this dependency of the standard deviations can be described, for statistically independent time samples, as

$$\begin{split}\sigma _{| E |}^2(f) &= \frac{1}{{{{| {E(f)} |}^2}}}\sum\limits_{k = 0}^{N - 1} ({[{{E_R}(f)\cos (2\pi\! fk\tau)} }\\& \quad - {E_I}(f)\sin (2\pi\! fk\tau) ]^2\sigma _E^2(k) )\end{split}$$
for the magnitude
$$| {E(f)}| = \sqrt {{E_R}{{(f)}^2} + {E_I}{{(f)}^2}}$$
and
$$\begin{split}\sigma _\vartheta ^2(f) &= \frac{1}{{{{| {E(f)} |}^4}}}\sum\limits_{k = 0}^{N - 1} ({[{{E_I}(f)\cos (2\pi\! fk\tau)}} \\&\quad + {E_R}(f)\sin (2\pi\! fk\tau) ]^2\sigma _E^2(k) )\end{split}$$
for the phase
$$\vartheta (f) = \arctan ({{E_I}(f)/{E_R}(f)}),$$
with ${E_R}$ and ${E_I}$ representing the real and imaginary parts of the discrete Fourier transform of the time-domain transient. The requirement for independent time samples can be easily realized in simulated scenarios (Section 3.A). Please note that, in the case of the experimentally measured data presented later (Section 3.B), the evaluation of the cross-correlation of the real and imaginary components of the Fourier-transformed data indicated that this independence requirement is not fulfilled in the evaluated data sets. Hence, the standard deviation in the frequency domain cannot be derived by the previously described concept but has to be calculated from the Fourier transforms of the individual measured transients [21].

The spectral signal-to-noise ratio can, in both cases, be defined as a ratio of the signal magnitude and standard deviation [21] as

$${\rm SNR}(f) = \frac{{| {E(f)} |}}{{{\sigma _{| E |}}(f)}}.$$
Thereby, the SNR and measure of the resulting phase error at a given frequency can directly be determined from the ratio of the measured data and its standard deviation.

The quality of the extracted parameters ${n_{\rm{samp}}}$ and ${\kappa _{\rm{samp}}}$ is evaluated by computing their root mean square error (RMSE), with respect to the exact reference values ${n_{\rm{ref}}}$ and ${\kappa _{\rm{ref}}}$ of the testing scenario. The RMSE for $N$ frequency points is defined as

$${\rm RMSE} = \sqrt {\sum\limits_{l = 1}^N \frac{{({{{[{{n_{\rm{extr}}}(l) - {n_{\rm{ref}}}(l)} ]}^2} + {{[{{\kappa _{\rm{extr}}}(l) - {\kappa _{\rm{ref}}}(l)} ]}^2}} )}}{{2N}}} .$$

This global measure of error is complemented by the investigation of the relative error of the parameters at every frequency point

$${\delta _n}(f) = \left| {1 - \frac{{{n_{\rm{samp}}}(f)}}{{{n_{\rm{ref}}}(f)}}} \right|$$
and
$${\delta _\kappa}(f) = \left| {1 - \frac{{{\kappa _{\rm{samp}}}(f)}}{{{\kappa _{\rm{ref}}}(f)}}} \right|,$$
providing direct information on spectral segments with an increased error, which can be related to the respective SNR and its effect on the extraction processes. The following analysis and discussion are focused on the SNR rather than the dynamic range (DR) because the SNR includes all noise effects overlaying the measured signal data (e.g.,  frequency- or power-dependent noise, contained in the standard deviation ${\sigma _{| E |}}(f)$ of the signal) and thereby all possible degradatory influences on material parameter extraction. An analysis with respect to DR and noise floor yields comparable results but would only partially include noise effects (for a detailed discussion, cf. [21]).

3. RESULTS AND DISCUSSION

A. Performance Evaluation by Simulated Data Sets

To evaluate the different methods for parameter extraction, two categories of testing scenarios are used. The first applies the methods to simulated noisy transient data sets and is presented in Section 3.A. The second, shown in Section 3.B, processes actually measured data at varying experimental signal-to-noise situations. All calculations (cf. Section 2) were implemented and performed using MATLAB and the lsqnonlin subroutine of its optimization toolbox.

The use of simulated data provides complete control over the pulse shape, the spectral characteristics of the sample, and the noise levels, independently from a particular measurement system and with full reproducibility. Since the material parameters are known by design, they can directly be quantitatively compared with the values extracted by the different, investigated methods. This provides an ideal direct measure of the errors of an analytic numerical approach and consequently its reliability for various SNR levels and sample properties.

 figure: Fig. 2.

Fig. 2. Simulated air reference and sample transients.

Download Full Size | PDF

The basis for the simulated data is a replica of an ideal THz transient, calculated from the derivative of a Gaussian function, similarly to the one shown in [22]. This transient is used as the air reference transient (Fig. 2). It has properties resembling those measured by the actual experimental system used later. The sample transient is obtained by convoluting this ideal air reference transient with the relative transfer function of the designed, simulated sample with four absorption resonances at 0.7, 1.25, 2.5, and 5.0 THz and a thickness of 900 µm. The corresponding calculated spectra are displayed in Fig. 3. In order to emulate the standard processing of measured data [10], typical preprocessing steps consisting of time-windowing the transient data, using a tapered cosine window, and applying a bandpass filter to the data in the frequency domain, according to the modeled and evaluated spectral range, are also performed. Normally distributed random numbers are used to generate the noise. White noise instead of field-amplitude-dependent noise [20,23] is used for the simulation, so that the effects of the noise can directly be interpreted as generally as possible, independently from a specific noise design. The amplitude noise is added to the simulated transients, which are scaled to different levels for the various investigated scenarios. These transients and their corresponding amplitude spectra are depicted in Figs. 4 and 5 along with the respective standard deviation, labeled according to their respective maximum SNR in the frequency domain.

 figure: Fig. 3.

Fig. 3. Simulated air reference and sample spectra with absorption resonances at 0.7, 1.25, 2.5, and 5.0 THz.

Download Full Size | PDF

 figure: Fig. 4.

Fig. 4. Simulated transients calculated at different signal amplitudes with noise maintained at a constant level. Inset shows an enlarged view of the ${{\rm{SNR}}_{{\max}}} = 42\,{\rm{dB}}$ case and the standard deviation in the time domain.

Download Full Size | PDF

 figure: Fig. 5.

Fig. 5. Amplitude spectra corresponding to the transients of different signal levels of Fig. 4 and the respective frequency-domain standard deviation.

Download Full Size | PDF

Figure 6 shows the spectral SNR of the generated data sets, with identical sample properties and different input signal levels. The general reduction of the SNR at the absorptive feature positions at 1.25 and 2.5 THz is clearly visible as well as the overall decline toward higher frequencies. In the ${{\rm{SNR}}_{{\max}}} = 82 \;{\rm{dB}}$ case, the SNR of the second feature at 2.5 THz drops just below 0 dB. With the further decrease of SNR, an almost 500 GHz wide segment drops below 0 dB.

 figure: Fig. 6.

Fig. 6. Spectral SNR of the simulated sample, labeled according to their maximum SNR (determined at the left dashed line). The right dashed line indicates the position of the minimum SNR, which first approaches 0 dB.

Download Full Size | PDF

In the following, the FDm and TDm material parameter extraction approaches are applied to the simulated data, as described in Section 2. The resulting RMSE of the parameters determined by the time- and frequency-domain material model-based methods are depicted in Fig. 7. The comparison of the RMSEs shows a clear overall superior quality of the parameters extracted by the time-domain approach. At the low SNR scenarios below ${{\rm{SNR}}_{{\max}}} = 62 \;{\rm{dB}}$, the frequency-domain method fails completely at the parameter extraction over a broad range of frequencies, indicated by the hollow points and the RMSE saturation at very high levels. The time-domain concept, on the other hand, still reliably determines the material parameters despite the calculating parameters from data with sub-0 dB SNR levels at the strongest resonance. Furthermore, it maintains an at least 40 times lower RMSE for all investigated noise cases, compared with its frequency-domain counterpart, as evident from the parallel lower RMSE values. With the minimum SNR (${{\rm{SNR}}_{{\min}}}$) at a spectral feature position close to 0 dB, as is noticeable in Fig. 6 for the ${{\rm{SNR}}_{{\max}}} = 82 \;{\rm{dB}}$ case at 2.5 THz, the time-domain method resulted in an RMSE of $4 \times {10^{- 5}}$ compared with the RMSE $2 \times {10^{- 3}}$ of the frequency-domain concept. Especially if the SNR of several segments of the spectral features lies below 0 dB, the time-domain processing maintains an RMSE better than $1 \times {10^{- 3}}$ and is able to provide correct data over the complete spectral range.

 figure: Fig. 7.

Fig. 7. RMSE of the parameters extracted at the SNR scenarios shown in Fig. 6, achieved by the material-model-based time-domain (TDm) and frequency-domain (FDm) methods, in black and red, respectively. The vertical dashed mark indicates the crossing of the SNR = 0 dB threshold by the strongest absorption feature.

Download Full Size | PDF

Figures 8 and 9 show the relative errors ${\delta _n}$ and ${\delta _\kappa}$ of the extracted refractive indices $n$ and extinction coefficients $\kappa$ of both approaches for the input data of Fig. 6 in reference to the ideal material parameters. It is evident by comparing the absolute values of the relative errors that the time-domain-extracted data are in significantly better agreement with the original data for all scenarios, including scenarios with low SNR. For instance, for the refractive index errors in Fig. 8 for the ${{\rm{SNR}}_{{\max}}} = 92 \;{\rm{dB}}$ case, the TDm has errors below $2 \times {10^{- 5}}$ for all frequency positions, while the FDm has errors reaching $1 \times {10^{- 3}}$. Furthermore, the results of the frequency-domain method show large deviations at low SNR. This is evident from a more or less constant spectral shape for the TDm as a function of decreasing SNR (only the error amplitude increases), while the spectral behavior of the FDm completely changes and rises up to $7 \times {10^{- 2}}$, where the spectral feature at 2.5 THz first drops below 0 dB, for the reduced ${{\rm{SNR}}_{{\max}}} = 62 \;{\rm{dB}}$ to ${{\rm{SNR}}_{{\max}}} = 42 \;{\rm{dB}}$ cases. In this situation, where the SNR of wide segments of the spectral feature lies below 0 dB, the relative errors, especially of the refractive index, show a broadband error and an offset like behavior of the extracted data following the spectral feature at 2.5 THz. The time-domain approach, on the other hand, maintains in the respective cases a maximum relative error of only $5 \times {10^{- 5}}$ in the ${{\rm{SNR}}_{{\max}}} = 82 \;{\rm{dB}}$ case and $5 \times {10^{- 3}}$ in the ${{\rm{SNR}}_{{\max}}} = 42 \;{\rm{dB}}$ case, and no broadband offset errors are occurring.

 figure: Fig. 8.

Fig. 8. Reference refractive index ${n_{\rm{ref}}}$ (top) and relative error of the refractive index ${\delta _n}$ determined by the material-model-based time-domain (TDm, center) and frequency-domain (FDm, bottom) methods. The spectral features positions are indicated by the red arrows.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Reference extinction coefficient ${\kappa _{\rm{ref}}}$ (top) and relative error of the extinction coefficient ${\delta _\kappa}$ by the material-model-based time-domain (TDm, center) and frequency-domain (FDm, bottom) methods. The spectral features positions are indicated by the red arrows.

Download Full Size | PDF

Since both methods process the same input data, same number of data points, and use identical settings and preprocessing steps, it stands to reason that the variation in performance, especially in the low SNR cases, stems from the different implementation of the optimization procedures. In the following, we will demonstrate that this difference results mainly from the handling of phase data in the frequency-domain approach. The standard deviation of the phase data is calculated with Eq. (12) and depicted for different SNR scenarios in Fig. 6. The increase in the phase standard deviation (or the phase error) with a decrease of SNR, especially at the regions of the spectral features, is clearly visible and reaching values close to $\pi$ and above for ${{\rm{SNR}}_{{\max}}} \le 82 \;{\rm{dB}}$. As a consequence, the phase values at the absorption features are highly affected by these effects. The inset in Fig. 10 depicts the ideal phase change between consecutive frequency points in the noise-free case. At the 2.5 THz feature, the change is ${-}{1.6}\;{\rm{rad}}$. For the ${{\rm{SNR}}_{{\max}}} = 82 \;{\rm{dB}}$ case, the standard deviation of the phase data at that point rises close to 3 rad. As the phase change of the noisy phase data is the sum of the noise-free phase change and the phase error, it is easy to see that, in this case, the ${\pm}\pi$ ambiguity threshold is exceeded. Consequently, an unwrapping based on this data cannot be performed properly, resulting in the introduction of an erroneous phase jump. Since several neighboring points also reach a critical standard deviation, several jumps can occur at this feature position. The number of those jumps drastically increases with the reduction of SNR. As those unwrapping errors cannot be distinguished from phase changes induced by the sample optical thickness, their removal by a manual or supervised unwrapping is mathematically not possible. Such errors, then, do not only affect the data locally at the spectral features, but, as they induce a phase offset to all succeeding phase points (due to the dependency of the phase on preceding spectral phase values in the unwrapping procedure), this error leads to global errors at all frequencies above these erroneous spectral positions. This completely inhibits correct phase handling in the optimization procedure and is the origin of the broadband errors observed earlier. As the material-model-based method in the frequency domain operates simultaneously over the complete investigated range, all determined material parameters are made erroneous by this effect.

 figure: Fig. 10.

Fig. 10. Standard deviation ${\sigma _\vartheta}$ of the spectral phase-information, corresponding to data sets of Fig. 6. Inset shows the ideal phase change between consecutive frequency points of the noise-free scenario.

Download Full Size | PDF

Since the time-domain method does not require any type of phase processing, as it directly performs the optimization step in the time domain, it is not affected by this limitation at all. Thereby, the effects of the rise of the error at the absorption features position remain spectrally localized and have only a limited degradatoriy effect on the overall extraction process. Broadband extraction of material properties thus becomes possible at limited SNR situations.

B. Performance Evaluation by Experimentally Measured Data Sets

The second section of the performance evaluation of the time- and frequency-domain methods is their application to actually measured, experimental transient data. Here, we investigate the data obtained from alpha-lactose samples pressed into pellets of 13.4 mm diameter and a thickness of 660 µm. The pellets were prepared without the use of a matrix material [24]. The measurements were performed with an in-house built TDS system using a biased GaAs photoconductive switch emitter and an electro-optic (GaP) differential detector in a lock-in measurement scheme, driven with a commercial 20 fs pulsed Ti:Sapphire laser with a repetition rate of 85 MHz (Femtolasers INTEGRAL 20). The THz pulses are guided by a set of off-axis parabolic mirrors. To emulate different SNR scenarios, the system was operated at various optical pump-power levels, resulting in a set of measurements with different THz field amplitudes. These transients, their amplitude spectrum, and the respective standard deviations are shown in Figs. 11 and 12. They are labeled according to their respective maximum SNR in the frequency domain. The calculations, as described in Section 2, are implemented and performed in the same fashion as for the simulated data, again using the lsqnonlin subroutine of MATLAB for optimization. The number of oscillators of the Lorentz model and their rough position is set accordingly to three absorption features directly discernible in the spectrum and a fourth feature visible as an absorption onset located outside the investigated frequency range.

 figure: Fig. 11.

Fig. 11. Measured transients of the lactose sample at different THz field levels (top) and their respective standard deviations (bottom).

Download Full Size | PDF

 figure: Fig. 12.

Fig. 12. Amplitude spectra and standard deviations corresponding to the transients in Fig. 11.

Download Full Size | PDF

In an analogous fashion to the evaluation of the simulated data sets, Fig. 13 depicts the spectral SNR of the lactose sample measurements. The inset shows the resulting RMSE values. To exclude effects of varying sample preparation, measurement system, and data handling [25], with regard to published lactose data, the RMSE is calculated in reference to the time-domain extracted results of the highest SNR case. Similarly to the evaluation of simulated data, the time-domain-approach TDm outperforms the frequency-domain FDm at a moderate minimum SNR with a RMSE of $2.2 \times {10^{- 3}}$ versus $1.4 \times {10^{- 2}}$ also in the real measurement case. Even for the ${{\rm{SNR}}_{{\max}}} = 25 \;{\rm{dB}}$ case with the SNR at the feature position close to 0 dB, the time-domain method provides reliable parameters with an RMSE of $7 \times {10^{- 3}}$, while the frequency-domain RMSE is 20 times higher and thereby insufficient for a proper parameter determination. The slightly reduced difference in performance of the measurement case compared with the simulation case, where a difference of 40 was observed, can most likely be attributed to noise effects in the measurement that are not fully represented by the simulated noise.

 figure: Fig. 13.

Fig. 13. Spectral SNR of the lactose sample measurements, labeled with the corresponding maximum SNR value. Inset shows the resulting RMSE of the respective extracted data.

Download Full Size | PDF

The extracted refractive indices are plotted for comparison in Fig. 14. For all time-domain cases and the high SNR frequency-domain case, all curves of the refractive index and the extinction coefficient agree closely around the ideal values. It is directly discernible that, for the low SNR case determined in the frequency domain, the spectral feature at 1.3 THz is missing in the extracted data, and an overall offset error is manifest.

 figure: Fig. 14.

Fig. 14. Complex refractive index of the lactose sample extracted by material-model-based time-domain (TDm) and frequency-domain (FDm) methods and the classical approach without the material model (FD) at the different SNR scenarios shown in Fig. 13.

Download Full Size | PDF

Following the same analysis of the simulated case, Fig. 15 shows the standard deviation of the measured phase data. The inset depicts the phase change between consecutive frequency points, which significantly peaks to a value of 2.9 rad at the sample feature position at 1.35 THz. A rise of the standard deviation at the absorption feature positions and an overall increase with decreasing SNR is again visible. For the ${{\rm{SNR}}_{{\max}}} = 38 \;{\rm{dB}}$ case, the ambiguity limit is approached due to the noise effects but not exceeded; therefore, the unwrapping can still be performed correctly, and an extraction including the feature in the frequency domain is possible. Still, the results are noticeably affected by the noise, resulting in a broadened extinction coefficient and a heightened change in the refractive index. For the low SNR case, the standard deviation at the feature position at 1.35 THz rises to 0.5 rad. Since this noise effect overlays with the 2.9 rad phase change caused by the samples’ properties, the overall phase change exceeds the ${\pm}\pi$ ambiguity limit, and the spectral feature is falsely removed from the phase data by the unwrapping procedure.

 figure: Fig. 15.

Fig. 15. Standard deviation ${\sigma _\vartheta}$ of the spectral phase information, corresponding to data sets of Fig. 13. Inset shows the phase change between consecutive frequency points of the high signal reference.

Download Full Size | PDF

As a consequence of the false phase representation, the samples’ spectral properties cannot be correctly reconstructed by the material model, as it describes the full investigated spectral range simultaneously and is thereby globally affected and shifted, resulting in complete failure of the optimization of the calculated and measured phase data. This inhibits proper parameter extraction by the frequency-domain method. These findings are consistent with those of the simulated scenarios, confirming the sensitivity of the frequency-domain method to phase errors, due to the resulting failure in phase-unwrapping and the superior robustness achieved by the use of the time-domain approach.

4. SUMMARY

We have shown that, by use of the material-model-based time-domain extraction method, material parameters can be extracted from THz transient data with an approximate factor of 20–40 lower root-mean-square error (RMSE), compared with the respective frequency-domain method. This was demonstrated by applying the data analysis procedures to simulated transient data and to experimentally measured data sets of alpha lactose at varying SNR conditions, utilizing only general a priori knowledge (visible from the spectra) on the number of rough positions of resonances for the adequate material model selection. The superiority of the time-domain concept is even more relevant in scenarios where the SNR is close or even below 0 dB at the absorption spectral feature position. While in such a case, time-domain methods still allow correct parameter extraction over the full measurement bandwidth, frequency-based approaches deliver erroneous results for all frequencies if one spectral position reaches a low SNR level.

In the demonstrated examples, the time-domain method extracted material parameters with a 40 times lower RMSE from the simulated data, compared with the results of the frequency-domain approach. In the case of the measured data, a 20 times lower RMSE was achieved by the time-domain concept. This improvement of noise robustness, on the one hand, further expands the capabilities of the THz-TDS analysis of samples with particularly strong resonant absorption features. On the other hand, it offers the option to speed up analytic procedures and measurements, since experimental averaging could be reduced at the cost of decreased SNR, while still maintaining a reasonable parameter extraction quality.

The cause for the performance difference was investigated by analyzing the noise propagation from the time to frequency domains. In the simulated and measured scenarios, a drastic rise in the standard deviation of the phase data at the resonant absorption features positions could be observed. In the context of the necessity of phase unwrapping in the frequency-domain method, this has a strong impact. If the ${\pm}\pi$ ambiguity limits are exceeded due to phase errors, the correct operation of the unwrapping procedure is inhibited, leading to a broadband offset error of the phase over a wide spectral range, which prevents proper parameter extraction by the frequency-domain approach. As the unwrapping is not required in the time-domain optimization concept, it remains reliable, even if the SNR locally decreases at specific frequencies due to an absorption feature, without adding further processing steps or weighting procedures.

Funding

Allianz Industrie Forschung (ZF4312305DF9).

Acknowledgment

We gratefully acknowledge useful discussions and experimental support by G. Schulte, T. A. Pham Tran, and Dr. R. Bornemann.

Disclosures

The authors declare no conflicts of interest.

REFERENCES

1. C. Jansen, S. Wietzke, O. Peters, M. Scheller, N. Vieweg, M. Salhi, N. Krumbholz, C. Jördens, T. Hochrein, and M. Koch, “Terahertz imaging: applications and perspectives,” Appl. Opt. 49, E48–E57 (2010). [CrossRef]  

2. T. Pfeifer, H. M. Heiliger, T. Löffler, C. Ohlhoff, C. Meyer, G. Lupke, H. G. Roskos, and H. Kurz, “Optoelectronic on-chip characterization of ultrafast electric devices: measurement techniques and applications,” IEEE J. Sel. Top. Quantum Electron. 2, 586–604 (1996). [CrossRef]  

3. P. H. Siegel, “Terahertz technology in biology and medicine,” IEEE Trans. Microw. Theory Tech. 52, 2438–2447 (2004). [CrossRef]  

4. M. C. Kemp, P. F. Taday, B. E. Cole, J. A. Cluff, A. J. Fitzgerald, and W. R. Tribe, “Security applications of terahertz technology,” Proc. SPIE 5070, 44–52 (2003). [CrossRef]  

5. K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue, “Non-destructive terahertz imaging of illicit drugs using spectral fingerprints,” Opt. Express 11, 2549–2554 (2003). [CrossRef]  

6. K. Yamamoto, M. Yamaguchi, F. Miyamaru, M. Tani, M. Hangyo, T. Ikeda, A. Matsushita, K. Koide, M. Tatsuno, and Y. Minami, “Noninvasive inspection of C-4 explosive in mails by terahertz time-domain spectroscopy,” Jpn. J. Appl. Phys. 43, L414–L417 (2004). [CrossRef]  

7. K. Lien Nguyen, T. Friscic, G. M. Day, L. F. Gladden, and W. Jones, “Terahertz time-domain spectroscopy and the quantitative monitoring of mechanochemical cocrystal formation,” Nat. Mater. 6, 206–209 (2007). [CrossRef]  

8. D. Mittleman, R. Jacobsen, R. Neelamani, R. Baraniuk, and M. Nuss, “Gas sensing using terahertz time-domain spectroscopy,” Appl. Phys. B 67, 379–390 (1998). [CrossRef]  

9. J. R. Freeman, H. E. Beere, and D. A. Ritchie, Generation and Detection of Terahertz Radiation (Springer, 2013), Chap. 1, pp. 1–28.

10. L. Duvillaret, F. Garet, and J. L. Coutaz, “A reliable method for extraction of material parameters in terahertz time-domain spectroscopy,” IEEE J. Sel. Top. Quantum Electron. 2, 739–746 (1996). [CrossRef]  

11. I. Pupeza, R. Wilk, and M. Koch, “Highly accurate optical material parameter determination with THz time-domain spectroscopy,” Opt. Express 15, 4335–4350 (2007). [CrossRef]  

12. M. Walther, B. Fischer, M. Schall, H. Helm, and P. Jepsen, “Far-infrared vibrational spectra of all-trans, 9-cis and 13-cis retinal measured by THz time-domain spectroscopy,” Chem. Phys. Lett. 332, 389–395 (2000). [CrossRef]  

13. O. S. Ahmed, M. A. Swillam, M. H. Bakr, and X. Li, “Efficient optimization approach for accurate parameter extraction with terahertz time-domain spectroscopy,” J. Lightwave Technol. 28, 1685–1692 (2010). [CrossRef]  

14. G. P. Kniffin and L. M. Zurk, “Model-based material parameter estimation for terahertz reflection spectroscopy,” IEEE Trans. Terahertz Sci. Technol. 2, 231–241 (2012). [CrossRef]  

15. P. U. Jepsen and B. M. Fischer, “Dynamic range in terahertz time-domain transmission and reflection spectroscopy,” Opt. Lett. 30, 29–31 (2005). [CrossRef]  

16. J. L. M. van Mechelen, A. B. Kuzmenko, and H. Merbold, “Stratified dispersive model for material characterization using terahertz time-domain spectroscopy,” Opt. Lett. 39, 3853–3856 (2014). [CrossRef]  

17. D. Stock and P. H. Bolívar, “Comparison of model-based material parameter extraction in frequency- and time-domain,” in 40th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz) (2015), pp. 1–2.

18. D. Stock and P. H. Bolívar, “Error analysis of model-based frequency- and time-domain methods for THz material parameter extraction,” in 41st International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz) (2016), pp. 1–2.

19. J. M. Fornies-Marquina, J. Letosa, M. Garcia-Gracia, and J. M. Artacho, “Error propagation for the transformation of time domain into frequency domain,” IEEE Trans. Magn. 33, 1456–1459 (1997). [CrossRef]  

20. W. Withayachumnankul, H. Lin, S. P. Mickan, B. M. Fischer, and D. Abbott, “Analysis of measurement uncertainty in THz-TDS,” Proc. SPIE 6593, 659326 (2007). [CrossRef]  

21. M. Naftaly and R. Dudley, “Methodologies for determining the dynamic ranges and signal-to-noise ratios of terahertz time-domain spectrometers,” Opt. Lett. 34, 1213–1215 (2009). [CrossRef]  

22. L. Duvillaret, F. Garet, J. F. Roux, and J. L. Coutaz, “Analytical modeling and optimization of terahertz time-domain spectroscopy experiments, using photoswitches as antennas,” IEEE J. Sel. Top. Quantum Electron. 7, 615–623 (2001). [CrossRef]  

23. L. Duvillaret, F. Garet, and J.-L. Coutaz, “Influence of noise on the characterization of materials by terahertz time-domain spectroscopy,” J. Opt. Soc. Am. B 17, 452–461 (2000). [CrossRef]  

24. B. M. Fischer, H. Helm, and P. U. Jepsen, “Chemical recognition with broadband THz spectroscopy,” Proc. IEEE 95, 1592–1604 (2007). [CrossRef]  

25. M. Naftaly and J. Molloy, “A multi-lab intercomparison study of THz time-domain spectrometers,” in 40th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz) (2015), pp. 1–2.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (15)

Fig. 1.
Fig. 1. Schematics of the classical frequency-based approach (left, FD), evaluating each frequency point individually, and the material-model-based approaches in the frequency (center, FDm), and time domains (right, TDm), indicating the measured and determined data (rounded boxes) and the numerical processing steps.
Fig. 2.
Fig. 2. Simulated air reference and sample transients.
Fig. 3.
Fig. 3. Simulated air reference and sample spectra with absorption resonances at 0.7, 1.25, 2.5, and 5.0 THz.
Fig. 4.
Fig. 4. Simulated transients calculated at different signal amplitudes with noise maintained at a constant level. Inset shows an enlarged view of the ${{\rm{SNR}}_{{\max}}} = 42\,{\rm{dB}}$ case and the standard deviation in the time domain.
Fig. 5.
Fig. 5. Amplitude spectra corresponding to the transients of different signal levels of Fig. 4 and the respective frequency-domain standard deviation.
Fig. 6.
Fig. 6. Spectral SNR of the simulated sample, labeled according to their maximum SNR (determined at the left dashed line). The right dashed line indicates the position of the minimum SNR, which first approaches 0 dB.
Fig. 7.
Fig. 7. RMSE of the parameters extracted at the SNR scenarios shown in Fig. 6, achieved by the material-model-based time-domain (TDm) and frequency-domain (FDm) methods, in black and red, respectively. The vertical dashed mark indicates the crossing of the SNR = 0 dB threshold by the strongest absorption feature.
Fig. 8.
Fig. 8. Reference refractive index ${n_{\rm{ref}}}$ (top) and relative error of the refractive index ${\delta _n}$ determined by the material-model-based time-domain (TDm, center) and frequency-domain (FDm, bottom) methods. The spectral features positions are indicated by the red arrows.
Fig. 9.
Fig. 9. Reference extinction coefficient ${\kappa _{\rm{ref}}}$ (top) and relative error of the extinction coefficient ${\delta _\kappa}$ by the material-model-based time-domain (TDm, center) and frequency-domain (FDm, bottom) methods. The spectral features positions are indicated by the red arrows.
Fig. 10.
Fig. 10. Standard deviation ${\sigma _\vartheta}$ of the spectral phase-information, corresponding to data sets of Fig. 6. Inset shows the ideal phase change between consecutive frequency points of the noise-free scenario.
Fig. 11.
Fig. 11. Measured transients of the lactose sample at different THz field levels (top) and their respective standard deviations (bottom).
Fig. 12.
Fig. 12. Amplitude spectra and standard deviations corresponding to the transients in Fig. 11.
Fig. 13.
Fig. 13. Spectral SNR of the lactose sample measurements, labeled with the corresponding maximum SNR value. Inset shows the resulting RMSE of the respective extracted data.
Fig. 14.
Fig. 14. Complex refractive index of the lactose sample extracted by material-model-based time-domain (TDm) and frequency-domain (FDm) methods and the classical approach without the material model (FD) at the different SNR scenarios shown in Fig. 13.
Fig. 15.
Fig. 15. Standard deviation ${\sigma _\vartheta}$ of the spectral phase information, corresponding to data sets of Fig. 13. Inset shows the phase change between consecutive frequency points of the high signal reference.

Equations (17)

Equations on this page are rendered with MathJax. Learn more.

T s a m p ( f ) = E s a m p ( f ) E a i r ( f ) = H s a m p ( f ) E i n p ( f ) H a i r ( f ) E i n p ( f ) = H s a m p ( f ) H a i r ( f ) .
T s a m p _ c a l c ( f ) = 2 η a i r ( f ) η a i r ( f ) + η s a m p ( f ) exp ( i η s a m p ( f ) 2 π f l s a m p c 0 ) 2 η s a m p ( f ) η s a m p ( f ) + η a i r ( f ) exp ( i η a i r ( f ) 2 π f l s a m p c 0 ) 1 .
η s a m p ( f ) = n s a m p ( f ) i κ s a m p ( f )
T s a m p _ m e a s ( f ) = E s a m p _ m e a s ( f ) E a i r _ m e a s ( f ) = T s a m p _ c a l c ( f ) .
Δ T R = | T s a m p _ c a l c ( f ) | | T s a m p _ m e a s ( f ) |
Δ T a r g = a r g ( T s a m p _ c a l c ( f ) ) a r g ( T s a m p _ m e a s ( f ) ) ,
η 2 ( f ) = ε ( f ) = ε + a P a 2 ( 2 π f a ) 2 ( 2 π f ) 2 + i γ a 2 π f .
E s a m p _ m e a s ( t ) = T s a m p _ c a l c ( t ) E r e f _ m e a s ( t ) = F 1 ( T s a m p _ c a l c ( f ) E r e f _ m e a s ( f ) ) ,
Δ E s a m p ( t ) = E s a m p _ m e a s ( t ) F 1 ( T s a m p _ c a l c ( f ) E r e f _ m e a s ( f ) ) .
σ | E | 2 ( f ) = 1 | E ( f ) | 2 k = 0 N 1 ( [ E R ( f ) cos ( 2 π f k τ ) E I ( f ) sin ( 2 π f k τ ) ] 2 σ E 2 ( k ) )
| E ( f ) | = E R ( f ) 2 + E I ( f ) 2
σ ϑ 2 ( f ) = 1 | E ( f ) | 4 k = 0 N 1 ( [ E I ( f ) cos ( 2 π f k τ ) + E R ( f ) sin ( 2 π f k τ ) ] 2 σ E 2 ( k ) )
ϑ ( f ) = arctan ( E I ( f ) / E R ( f ) ) ,
S N R ( f ) = | E ( f ) | σ | E | ( f ) .
R M S E = l = 1 N ( [ n e x t r ( l ) n r e f ( l ) ] 2 + [ κ e x t r ( l ) κ r e f ( l ) ] 2 ) 2 N .
δ n ( f ) = | 1 n s a m p ( f ) n r e f ( f ) |
δ κ ( f ) = | 1 κ s a m p ( f ) κ r e f ( f ) | ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.