1. Introduction

At the heart of most applications of terahertz time-domain spectroscopy (THz-TDS) is a mathematical procedure that relates raw THz-TDS waveform measurements to parameters of scientific and technological interest [1–3]. Typically this analysis is formulated in the frequency domain, since it provides the most natural description of any linear, time-invariant system of interest. But since a THz-TDS waveform describes the electric field as a function of time, it must be Fourier-transformed for use in any frequency-domain analysis. The Fourier transform scrambles the THz-TDS measurement noise, which is more naturally represented in the time domain, and produces artifacts that can degrade the overall measurement quality and yield misleading results [4–13]. Furthermore, the standard approaches to parameter estimation in THz-TDS involve the ratio of two noisy waveforms, which further distorts the noise spectrum and can yield noise distributions that are far from Gaussian. Other approaches have emerged that improve on the standard procedures [14–16], but so far all of the existing approaches to THz-TDS analysis lack clear, model-independent statistical criteria for establishing parameter confidence intervals or for deciding whether a given physical model describes the data well or not. Here, we introduce a general framework to remedy this [17]. We describe a maximum-likelihood approach to THz-TDS analysis in which both the signal and the noise are treated explicitly in the time domain, together with a frequency-domain constraint between the input and the output signal. We show that this approach produces superior results to existing analysis methods.

2. Signal and noise in THz-TDS

A Monte Carlo simulation of an elementary THz-TDS analysis illustrates the basic problems that our method solves; it also allows us to develop notational conventions that we will use to describe our main results. Figure 1(a) shows an ideal, noiseless THz-TDS waveform $\mu (t)$ normalized to its peak value ${\mu _{\textrm{p}}}$ [18], together with a time-domain Gaussian noise amplitude $\sigma _\mu (t)$ given by

 

Fig. 1. (a) Simulated time-domain signal (black solid line) and noise amplitude (gray dashed line). (b) Power spectral density (relative to peak) obtained from ten simulated time-domain measurements, with one shown in red and the rest shown in gray.

Download Full Size | PDF

We simulate a single THz-TDS waveform measurement by computing

3. Transfer function estimation in THz-TDS

To measure the response of a system requires measurements of two THz-TDS waveforms: an input, $\mu (t)$, and an output, $\psi (t) = [h{*}\mu ](t)$, where $h(t)$ is the impulse response of the system and the ${*}$ symbol denotes continuous-time convolution. Fourier transforming the ideal relationship $\psi (t) = [h{*}\mu ](t)$ converts the time-domain convolution to a frequency-domain multiplication, $\tilde {\psi }_{\textrm{c}}(\omega ) = H(\omega )\tilde {\mu }_{\textrm{c}}(\omega )$, where $\tilde {\mu }_{\textrm{c}}(\omega )$ and $\tilde {\psi }_{\textrm{c}}(\omega )$ denote the continuous-time Fourier transforms of $\mu (t)$ and $\psi (t)$, respectively, and $H(\omega ) = \tilde {h}_{\textrm{c}}(\omega )$ is the continuous-time transfer function of the system. Proceeding as we did for a single waveform, we simulate input and output waveform measurements $x(t_k)$ and $y(t_k)$, respectively, by computing $\mu (t)$ and $\psi (t)$ at the discrete times $t_k$ and adding time-domain noise $\sigma _\mu (t_k)\epsilon (t_{k})$ and $\sigma _\psi (t_k)\delta (t_k)$, where each $\epsilon (t_k)$ and $\delta (t_k)$ is an independent random variable with a standard normal distribution. After applying the discrete Fourier transform to obtain $X(\omega _l) = \tilde {x}(\omega _l)$ and $Y(\omega _l) = \tilde {y}(\omega _l)$, we can determine the empirical transfer function estimate (ETFE) [19],

To illustrate these deficiencies, Fig. 2 shows 250 simulations of $\hat {H}_{\textrm{ETFE}}(\omega _l)$ with nominally identical input and output pulses, corresponding to $H(\omega ) = 1$. The red-highlighted curves show correlations similar to those shown in Fig. 1, but that extend to frequencies where the signal is much stronger. The average over all simulations reveals the bias in $\hat {H}_{\textrm{ETFE}}(\omega _l)$: in Fig. 2(a), $\textrm{Re} \{\hat {H}_{\textrm{ETFE}}\}$ falls from the correct value of 1 to the biased value of 0 as the frequency increases above 5 THz, where the signal reaches the noise floor in Fig. 1(b). Also, the noise distribution for $\hat {H}_{\textrm{ETFE}}(\omega _l)$ develops wide, non-Gaussian tails at frequencies where the signal is small, because noise fluctuation that cause $|X(\omega _l)|\rightarrow 0$ become more likely, and these in turn cause $\hat {H}_{\textrm{ETFE}}(\omega _l)$ to diverge [20]. If the noise is uncorrelated, such large fluctuations are only expected when the signal is small. But in the more typical case that the noise is correlated, they can influence frequencies well within the primary signal bandwidth, as indicated by the red-highlighted curves in Figs. 2(c) and 2(d).

 

Fig. 2. Gray dots show the real (a,c) and imaginary parts (b,d) of the empirical transfer function estimate $\hat {H}_{\textrm{ETFE}}$, Eq. (4), for 250 pairs of simulated time-domain measurements of the waveform shown in Fig. 1(a). One estimate is highlighted in red, with the dots connected by a thin red line. The solid black line shows the average over all simulations. Panels (a,b) show the full bandwidth and (c,d) show the same data over the primary signal bandwidth.

Download Full Size | PDF

Figure 3 shows how these problems distort standard measures of fit quality. It displays the results of weighted least-squares fits to the ETFE simulations in Fig. 2 with the model

 

Fig. 3. Measures of fit quality for ETFE fits with Eq. (5), obtained by minimizing Eq. (6). (a) Cumulative distribution of the GOF statistic, $S_{\textrm{ETFE}}$, for the Monte Carlo simulations shown in Fig. 2. The $\chi ^{2}$ distribution for the same number of degrees of freedom ($\nu = 254$) is shown for comparison. The red $\times$ marks a fit with $S_{\textrm{ETFE}} \approx 180$, which is just above the median value. The normalized residuals for this fit, $r_{\textrm{ETFE}}(\omega _l; {\boldsymbol {\hat {\theta }}})$, are shown in (b) as a function of frequency, and in the inset to (a) with a normal probability plot (black dots, which represent both the real and the imaginary parts of $r_{\textrm{ETFE}}$). The gray dash-dotted line in the inset represents the standard normal distribution.

Download Full Size | PDF

An alternative approach is to use a least-squares (LS) procedure that minimizes the sum of the squared differences between the output signal and the transformed input signal [14–16],

In principle, we should be able to account for time-dependent noise by assigning appropriate weights to the sum, but this is not as simple as it might seem. The problem is that Eq. (8) treats the input and output noise asymmetrically, which we can see more clearly if we express it explicitly in terms of the time-domain signal and noise sequences:

4. Maximum-likelihood estimation of a parameterized transfer function model

The likelihood function is equivalent to the probability density for obtaining the measured data under the assumptions of a given model, and forms the theoretical basis for the standard least-squares fit procedure. To define it in the THz-TDS context, we express the measured input and output signals in vector notation as ${\boldsymbol {x}} = [x(t_0), x(t_1), \ldots , x(t_{N-1})]^{\mathsf {T}}$ and ${\boldsymbol {y}} = [y(t_0), y(t_1), \ldots , y(t_{N-1})]^{\mathsf {T}}$, respectively, and assume that they represent noisy measurements of bandlimited but otherwise unknown ideal signal vectors ${\boldsymbol {\mu }} = [\mu (t_0), \mu (t_1), \ldots , \mu (t_{N-1})]^{\mathsf {T}}$ and ${\boldsymbol {\psi }} = [\psi (t_0), \psi (t_1), \ldots , \psi (t_{N-1})]^{\mathsf {T}}$ with noise covariance matrices ${\boldsymbol {\Sigma _x}}$ and ${\boldsymbol {\Sigma _y}}$, respectively. We further assume that ${\boldsymbol {\mu }}$ and ${\boldsymbol {\psi }}$ satisfy the linear constraint ${\boldsymbol {\psi }} = {\boldsymbol {h}}({\boldsymbol {\theta }}){\boldsymbol {\mu }}$, where ${\boldsymbol {h}}$ is a known matrix function of an unknown $p$-dimensional parameter vector ${\boldsymbol {\theta }}$. The likelihood function is then a product of two multivariate Gaussian distributions,

We can use discrete-time Fourier analysis to obtain explicit expressions for ${\boldsymbol {h}}({\boldsymbol {\theta }})$, ${\boldsymbol {\Sigma _x}}$, and ${\boldsymbol {\Sigma _y}}$. The time-domain transfer matrix ${\boldsymbol {h}}({\boldsymbol {\theta }})$ is

The maximum-likelihood (ML) estimate for the model parameters is given by the vectors ${\boldsymbol {\hat {\mu }}}$, ${\boldsymbol {\hat {\psi }}}$, and ${\boldsymbol {\hat {\theta }}}$ that maximize $L$ subject to the constraint ${\boldsymbol {\psi }} = {\boldsymbol {h}}({\boldsymbol {\theta }}){\boldsymbol {\mu }}$. Alternatively, we can minimize the negative-log likelihood,

Introducing an $N$-dimensional vector of Lagrange parameters ${\boldsymbol {\lambda }}$ to implement the constraint, we can define the modified cost function,

We can simplify these expressions further by defining

Substituting Eq. (26) in Eq. (25) yields

Figure 4 shows fit results for the same model and data shown in Fig. 3, but obtained by minimizing $Q_{\textrm{TLS}}({\boldsymbol {\theta }})$ in Eq. (30) instead of $Q_{\textrm{ETFE}}({\boldsymbol {\theta }})$ in Eq. (6). The statistical properties obtained with $Q_{\textrm{TLS}}({\boldsymbol {\theta }})$ are clearly superior to those with $Q_{\textrm{ETFE}}({\boldsymbol {\theta }})$. The GOF statistic, $S_{\textrm{TLS}} = Q_{\textrm{TLS}}({\boldsymbol {\hat {\theta }}})$, approximates a $\chi ^{2}$ distribution with $\nu = N-p$ degrees of freedom, as expected. The normalized residual vector ${\boldsymbol {r}}_{\textrm{TLS}}({\boldsymbol {\hat {\theta }}})$, shown for one of the fits, exhibits a standard normal distribution with no discernible correlations among neighboring time points. The distribution for $S_{\textrm{TLS}}$ is more than 7 times narrower than the distribution for $S_{\textrm{ETFE}}$, making it that much more sensitive to a poor fit. Similarly, the lack of structure in ${\boldsymbol {r}}_{\textrm{TLS}}({\boldsymbol {\hat {\theta }}})$ makes it more useful for residual analysis than $r_{\textrm{ETFE}}(\omega _l; {\boldsymbol {\hat {\theta }}})$, since structure associated with poor fits may be discerned more readily. Finally, unlike the ETFE fits, the TLS fits yield estimates for the parameter uncertainties that agree with the values observed over the set of Monte Carlo samples, $\sigma _{\theta _1} = 0.002$ and $\sigma _{\theta _2} = 0.4~\textrm{fs}$, which are both about a factor of 2 smaller than the values observed in the ETFE parameter estimates. In summary, when compared with the standard ETFE fit procedure, the TLS fit procedure offers better discrimination between good models and bad ones, more precise values for the parameters, and more accurate estimates of the parameter uncertainties.

 

Fig. 4. Measures of fit quality for TLS fits with Eq. (5), obtained by minimizing Eq. (30); compare with Fig. 3. (a) Cumulative distribution of the GOF statistic, $S_{\textrm{TLS}}$, for the Monte Carlo simulations shown in Fig. 2. The $\chi ^{2}$ distribution for the same number of degrees of freedom ($\nu = 254$) is shown for comparison. The red $\times$ marks a fit with $S_{\textrm{TLS}} \approx 252$, which is just above the median value. The normalized residuals for this fit, ${\boldsymbol {r}}_{\textrm{TLS}}({\boldsymbol {\hat {\theta }}})$, are shown in (b) as a function of time, and in the inset to (a) with a normal probability plot (black dots). The gray dash-dotted line in the inset represents the standard normal distribution.

Download Full Size | PDF

5. Maximum-likelihood estimation of the noise model

We have assumed so far that the noise amplitudes $\sigma _\alpha , \sigma _\beta ,$ and $\sigma _\tau$ in Eq. (13) are known, but in general they must also be determined experimentally. One way to do this is to measure the noise at three different points on the THz waveform: as we saw in Fig. 1, the $\sigma _\beta$ term dominates near the peak, the $\sigma _\tau$ term dominates where the signal crosses zero with a large slope, and the $\sigma _\alpha$ term dominates where both the signal and its slope are small. Another approach is to determine the variance as a function of time from a set of repeated waveform measurements [5], then obtain the noise parameters from a fit with Eq. (1). But both of these methods are susceptible to systematic errors from drift, which produces excess signal variability over the timescale of multiple measurements [8,9,23]. In this section, we describe a likelihood method for estimating the noise parameters that accounts for this drift.

We consider repeated measurements of an ideal primary waveform, $\mu (t)$, which has an amplitude and a delay that drift on a timescale longer than the waveform acquisition time. We can then associate each measurement with an ideal secondary waveform,

With these definitions, we can express the sampled secondary waveforms in vector form,

The likelihood function for observing ${\boldsymbol {X}}$ under these assumptions is

The resulting estimates for the noise parameters (${\sigma _\alpha }$, ${\sigma _\beta }$, and ${\sigma _\tau }$) are biased below their true values, which we can attribute to the presence of the nuisance parameters (${\boldsymbol {\mu }}$, ${\boldsymbol {A}}$, and ${\boldsymbol {\eta }}$) [21]. For example, in 1000 Monte Carlo simulations of $M=10$ repeated measurements using the waveforms described in Sec. 2, the ratios of the ML estimates to their true values are ${\hat {\sigma }_\alpha }/{\sigma _\alpha } = 0.95 \pm 0.02$, ${\hat {\sigma }_\beta }/{\sigma _\beta } = 0.94 \pm 0.03$, and ${\hat {\sigma }_\tau }/{\sigma _\tau } = 0.89 \pm 0.09$, all significantly below unity. Increasing the number of observations to $M=50$ reduces the bias to ${\hat {\sigma }_\alpha }/{\sigma _\alpha } = 0.990 \pm 0.007$, ${\hat {\sigma }_\beta }/{\sigma _\beta } = 0.98 \pm 0.01$, and ${\hat {\sigma }_\tau }/{\sigma _\tau } = 0.93 \pm 0.04$, but some bias remains in ${\hat {\sigma }_\beta }$ and ${\hat {\sigma }_\tau }$ even in the limit $M\rightarrow \infty$. This residual bias arises because the number of elements in ${\boldsymbol {A}}$ and ${\boldsymbol {\eta }}$ both grow with the number of observations, which allows us to account for drift but dilutes some of the information that the data provide about the noise [21,24]. In principle, we can resolve this problem by integrating out all of the nuisance parameters in Eq. (35) to obtain a marginal likelihood that depends only on ${\sigma _\alpha }$, ${\sigma _\beta }$, and ${\sigma _\tau }$ [21], but this involves computationally expensive integrations for a relatively small benefit. As we discuss below, we can remove much of the bias by simply rescaling the noise parameters by a suitable correction factor.

To determine the bias correction factor it is helpful to consider a simplified example in which there is no drift and only additive noise, so that $A_l = 1$ and $\eta _l = 0$ for all $l$ and ${\sigma _\beta } = {\sigma _\tau } = 0$. The ML cost function then reduces to

6. Experimental verification

In this section we present experimental results that verify our analysis. Figure 5(a) shows the raw data that we use to specify the noise model, ${\boldsymbol {X}}$, which comprises $M=50$ waveforms, each with $N=256$ points sampled every $T = 50~\textrm{fs}$. With this data as input, we minimize Eq. (37) to obtain the ML estimates, ${\hat {\sigma }_\alpha }$, ${\hat {\sigma }_\beta }$, ${\hat {\sigma }_\tau }$, $\hat {{\boldsymbol {\mu }}}$, $\hat {{\boldsymbol {A}}}$, and $\hat {{\boldsymbol {\eta }}}$ for all of the free parameters in the model. The resulting time-dependent noise amplitude, corrected for bias using Eq. (41), is

 

Fig. 5. Noise-model fits with laboratory measurements. (a) Set of 50 nominally identical waveform measurements that compose the data matrix ${\boldsymbol {X}}$. (b) Time-dependent noise amplitude obtained from a fit with Eq. (37) to ${\boldsymbol {X}}$. The solid line shows the ideal noise model, with the bias-corrected ML estimates ${\hat {\sigma }_\alpha }^{*} = (0.623\pm 0.005)~\textrm{pA}$, ${\hat {\sigma }_\beta }^{*} = (1.55\pm 0.03){\%}$, and ${\hat {\sigma }_\tau }^{*} = (1.8\pm 0.1)~\textrm{fs}$. Points show the standard deviation of the measurements after correcting for drift, as described in the text. (c) Normalized residuals for one of the waveform measurements shown in (a).

Download Full Size | PDF

Figure 6 shows fits to two sequential waveforms from this set. A fit with Eq. (5) in the time domain, obtained by minimizing $Q_{\textrm{TLS}}$ in Eq. (28), yields ${\hat {\theta }}^{\textrm{TLS}}_1 = 1.019\pm 0.003$ for the amplitude and ${\hat {\theta }}^{\textrm{TLS}}_2 = (-2.8\pm 0.5)~\textrm{fs}$ for the delay, with the normalized residuals shown in Fig. 6(b). A fit with the same model in the frequency domain, obtained by minimizing $Q_{\textrm{ETFE}}$ in Eq. (6), yields ${\hat {\theta }}^{\textrm{ETFE}}_1 = 1.020\pm 0.004$ and ${\hat {\theta }}^{\textrm{ETFE}}_2 = (-4.4\pm 0.6)~\textrm{fs}$, with the normalized residuals shown in Fig. 6(d). Both fit methods reveal significant differences from the ideal values of $\theta _1 = 1$ and $\theta _2 = 0$, reflecting the measurement drift discussed in Sec. 5. As we found for the residuals of the noise-model fit in Fig. 5(c), the residuals of the transfer-function fit in Fig. 6(b) show small but statistically significant correlations. But as we also found with the idealized Monte Carlo simulations, the residuals of the ETFE fit in the frequency domain are much more structured than the residuals of the TLS fit to the same data in the time domain.

 

Fig. 6. Transfer-function fits with laboratory measurements. (a) Two sequential waveform measurements, taken from the set shown in Fig. 5. (b) Normalized time-domain residuals, $\boldsymbol {r}_{\textrm{TLS}}({\boldsymbol {\hat {\theta }}}^{\textrm{TLS}})$, for the TLS fit with Eq. (5). (c) Real and imaginary parts of $\hat {H}_{\textrm{ETFE}}(\omega _l) - 1$ (dots), fitted with Eq. (5) (solid lines) and Eq. (45) (dotted lines) by minimizing $Q_{\textrm{ETFE}}$ in Eq. (6). (d) Normalized frequency-domain residuals, $r_{\textrm{ETFE}}(\omega _l; {\boldsymbol {\hat {\theta }}}^{\textrm{ETFE}})$, for the fit with Eq. (5).

Download Full Size | PDF

To illustrate the analysis problem that this raises, in Fig. 6(c) we compare an ETFE fit with Eq. (5) to an ETFE fit with a more flexible transfer-function model,

If we divide the 50 experimental waveforms into 25 sequential pairs and fit each pair with Eq. (5), the results are qualitatively consistent with the Monte Carlo simulations and support the conclusion that $S_{\textrm{TLS}}$ offers the better absolute measure of fit quality. The distribution of $S_{\textrm{ETFE}}$ over the experiments has a mean $\bar {S}_{\textrm{ETFE}} \approx 235$ and standard deviation $\sigma (S_{\textrm{ETFE}}) \approx 113$, while the distribution for $S_{\textrm{TLS}}$ has a mean $\bar {S}_{\textrm{TLS}} \approx 246$ and standard deviation $\sigma (S_{\textrm{TLS}}) \approx 39$. For the simulated distributions shown in Fig. 3(a) and Fig. 4(a), we obtain $\bar {S}_{\textrm{ETFE}} \approx 235$, $\sigma (S_{\textrm{ETFE}}) \approx 160$, $\bar {S}_{\textrm{TLS}} \approx 250$, and $\sigma (S_{\textrm{TLS}}) \approx 22$. As we discussed at the end of Sec. 4, a narrower GOF distribution provides better sensitivity to fit quality. And while the experimental distribution for $S_{\textrm{TLS}}$ is nearly twice as broad it is in the simulations, it is still nearly a factor of 3 narrower than the experimental distribution for $S_{\textrm{ETFE}}$. Despite the quantitative differences between the simulations and the experiment, the TLS method consistently shows better performance than the ETFE.

7. Conclusion

In summary, we have developed a maximum-likelihood approach to THz-TDS analysis and demonstrated that it has numerous advantages over existing methods. Starting from a few simple assumptions, we derived a method to determine the transfer function relationship between two THz-TDS measurements, together with a method to specify a parameterized noise model from a set of repeated measurements. In each case, we also derived expressions for fit residuals in the time domain that are properly normalized by the expected noise. Through a combination of Monte Carlo simulations and experimental measurements, we verified that these tools yield results that are accurate, precise, responsive to fit quality, and generally superior to the results of fits to the ETFE.

We focused on simple examples to emphasize the logical structure of the method, but we can readily apply the same approach to more complex problems. For example, we have successfully used the framework presented here to measure a weak frequency dependence in the Drude scattering rate of a metal, which is predicted by Fermi liquid theory; we have also used it to measure small variations in the THz-frequency refractive index of silicon with temperature [17]. In both of these applications we found that maximum-likelihood analysis in the time domain provided better performance than a least-squares analysis based on the ETFE in the frequency domain. We expect similar improvements in other applications, and provide MATLAB functions in Code Repository 1 (Ref. [25]) for others to try.