## Abstract

We improve the existing data extraction algorithms for THz time-domain spectroscopy (THz TDS) in two aspects. On the one hand, we merge the up-to-date knowledge of THz TDS signal processing into a single powerful optical material parameter extraction algorithm. On the other hand, we introduce a novel iterative algorithm that further enhances the accuracy of the parameter extraction. In contrast to most of the published experiments, we are able to reliably investigate samples with thicknesses as small as 100μm, samples with low indexes of refraction, i.e. close to 1, as well as samples with sharp peaks in the material parameter curves.

©2007 Optical Society of America

## 1. Introduction

*Terahertz time-domain spectroscopy* (THz TDS) is a technique for measuring the complex refractive index of materials over a wide range of frequencies, usually spanning from a few tens of GHz to several THz. In the last decade several authors presented material parameter extraction algorithms to determine the complex refractive indexes of nearly homogenous solid samples with THz TDS, see Ref. [1, 2, 3, 4, 5]. Commonly, for material parameter extraction, a THz pulse which has propagated through a sample is compared to a THz pulse obtained without the sample in its propagation path.

The sample material acts as a Fabry-Perot resonator for the THz pulse. Thus, the measured signal exhibits multiple reflections of the THz pulse. A simplistic method for approximating the optical material parameters considers a narrow time window for the sample waveform, containing only the first pulse propagated through the sample. However, the isolation of the first pulse may not be possible in the cases of low refractive indexes and/or thin samples. Another major drawback of this method is the necessity of a very accurate a priori knowledge of the sample thickness.

The recent papers Ref. [2] and Ref. [3] present algorithms leading to much more accurate results than the simple method mentioned above. Additionally, they allow for an accurate determination of the sample thickness. However, even though very large, the spectrum of materials suited for investigation with these algorithms is constrained by certain boundaries, such as relatively high refractive indexes and/or relatively thick samples (compare with the limits of the method in Ref. [3]). Furthermore, the accuracy of the existing algorithms is restricted by the uncertainty of the measured signals.

In this paper we present two important advancements to THz TDS. On the one hand we improve the existing state of the art material parameter extraction algorithms. We merge the up-to-date knowledge on THz TDS signal processing presented in Ref. [2, 3, 4] into a single powerful material parameter extraction algorithm. On the other hand we introduce a novel algorithm that enhances the accuracy of material parameter extraction by using the uncertainty of the measurements in a constructive way. We define a confidence interval for the theoretically derived transfer function, which allows for the formulation of additional physical boundary conditions to the investigated sample. We introduce the *Spatially Variant Moving Average Filter* (SVMAF), which constitutes the core of an iterative algorithm for finding the optimal material parameter curves within the allowed confidence intervals considering the additional physical boundaries.

With our algorithm we reach a very high accuracy of the extracted complex indexes of refraction of materials. In contrast to other published results, our algorithm performs very accurately even in the difficult cases of thin samples, with thicknesses as small as 100μm, and samples with low refraction indexes. Hence, it has a wide range of applicability.

## 2. The experimental setup

Our THz TDS setup is very similar to the setup described in Ref. [6]. The heart of our setup is a femtosecond titanium sapphire laser operated by a diode pumped solid state laser. It produces 30fs laser pulses with a repetition rate of 80MHz. The laser pulses are used to excite two photoconductive dipole antennas. One part of the optical beam is focused onto the emitter antenna which consists of GaAs onto which two metallic striplines are deposited. The laser pulses generate free carriers in between the striplines which are accelerated by an applied bias.

The resulting time-dependent photocurrent has a subpicosecond rise time and acts as a source for the subpicosecond THz pulses. These pulses are radiated from the antenna and are then collimated, focused and guided by four off-axis paraboloidal mirrors to from an intermediate focus in which the sample is placed. Another antenna is used to detect the THz radiation. In this antenna carriers are excited by the second part of the optical laser pulse which has passed a variable delay line. The electric field of the THz pulse is then used to bias the detector. The induced photocurrent is proportional to this electric field and the transient carrier population. The time-dependent photocurrent is measured by scanning the delay line. The time resolution of this measurement is limited mainly by the carrier lifetime in the silicon on sapphire (SOS) photoconductive antenna detector. For further details see Ref. [6].

## 3. Analytical model of the transfer function

In this section we derive an analytical formula for the transfer function. This is done in close analogy to Ref. [3]. We include this derivation here since it is of capital importance for the understanding of the subsequent sections. Throughout this section we highlight the aspects that we consider particularly important.

#### 3.1. The inverse electromagnetic problem

The method used by THz TDS to compute material parameters is based on measuring the effect caused on the electric field of a THz wave by the pass through the investigated material. Since the parameters are not computed directly, but from their effect on the electric field, this electromagnetic problem is called an *inverse* problem.

The temporal profile of the THz pulse propagated through a material represents the convolution in time domain of the emitted THz waveform with the pulse response of the material. The technique widely used in THz TDS is the comparison of two signals in frequency domain: a *reference signal E*
_{ref} (*ω*) and a *sample signal E*
_{samp}(*ω*), where *ω* denotes the angular frequency. Both electric fields in the frequency domain are directly proportional to the Fourier transforms of the respective photocurrents multiplied by *jω*. The sample signal is obtained by recording the photocurrent of the reference signal, additionally propagated through the sample material.

We shall denote by *E _{ref}^{ex}* (

*ω*) and

*E*(

^{ex}_{samp}*ω*) the recorded reference and sample signals, respectively. The frequency dependent experimentally determined transfer function of the sample material is given by their quotient, which corresponds to the deconvolution in the time domain:

Let *H _{theory}*(

*ω*) be a theoretically derived transfer function of the sample material, which depends on the material parameters. We will derive an analytical expression for

*H*(

_{theory}*ω*) in Section 3.3.

Let *n*
_{1}(ω) be the real refractive index and α_{1}(ω) be the absorption coefficient of a material we indicate with the index 1. As usual, we denote the complex refractive index of the material 1 by ñ_{1}(ω) = *n*
_{1}(ω) − *j*κ_{1}(ω), where κ_{1}(ω) is the extinction coefficient indicating the absorption. The relation between the absorption coefficient α_{1}(ω) and the extinction coefficient κ_{1}(ω) is κ_{1}(ω) = α_{1}(2ω)*c*
_{0}(2ω), where *c*
_{0} is the speed of light in air.

The inverse electromagnetic problem consists in finding the material parameters of the sample material, such that for every frequency ω in the considered frequency range, the value *H _{theory}*(ω) is as close as possible to

*H*(ω). As a measure for how close these two functions are, we define the following error function

_{experiment}*Err*:

Note that all computations in this work are made in discrete time and frequency domains. This explains the sum over all frequencies in equation (4).

Solving the inverse electromagnetic problem for a material 1 with unknown material parameters is therefore equivalent to minimizing the value of *Err* with respect to the functions *n*
_{1}(ω) and κ_{1}(ω). Since *ω* has discrete values, equation (4) implies that this problem is equivalent to minimizing the expression |*M*(ω)| + |*A*(ω)| with respect to *n*
_{1}(ω) and κ_{1}(ω) for each frequency ω separately.

#### 3.2. Assumptions and notation

In the subsequent section we derive an analytical formula for the transfer function. The assumptions made about the experimental setup are very important for the validity of this theoretical model. In the following, we give a list of conditions which we consider fulfilled throughout this paper as well as the standard notation that we use.

*Homogeneity of the investigated sample*This is tantamount with its material parameters being spatially and directionally invariant. However, the frequency dependence of the material parameters may be arbitrary. In other words, absorption peaks and other variations of the material parameters are detected without losses of the quality of the results.

*Flat and parallel surfaces of sample materials and no scattering*We neglect roughness and curvature of the irradiated materials, as well as scattering effects at the surfaces and inside the material. In our analytical transfer function derivation, we assume that all surfaces are perfectly flat and parallel.

*Dry atmosphere in the environment of the experiment*We assume that the experimental setup is purged with dry Nitrogen gas to avoid complications due to water vapor.

*Orthogonal incidence of the THz rays*For simplicity we presume that the THz rays impinge orthogonally on the irradiated samples. The sample is placed in the spectrometer with an inclination error smaller than 5°, which leads to deviations of the refractive index and the extinction coefficient smaller than 0.002. Thus, we neglect this error.

*Fresnel coefficients notation*With

*R*_{12}and*T*_{12}we denote the Fresnel reflection and transmission coefficients, respectively. They describe the reflection and transmission of a THz wave at the interface of material 1 and material 2 in the time domain:Furthermore, the propagation of an electric field

*E*entering the material 1 is governed by the following relation:_{init}$${P}_{1}\left(z\right)=\mathrm{exp}\left(\frac{-{\mathrm{j\omega}\tilde{n}}_{1}z}{{c}_{0}}\right),$$where

*E*(*z*) is the electric field after the geometrical propagation length*z*and*P*_{1}(*z*) is the frequency dependent propagation coefficient of the material 1.*Standard notation used throughout this work*For computations, we use the value

*c*= 2.99796 ∙ 10_{vacuum}^{8}m/s for the speed of light in vacuum. We approximate the speed of light*c*_{0}in air by cvacuum and for the complex refractive index of air we use the value*n*̃_{0}= 1.00027 −*j*_{0}. We always indicate the medium air with the index 0. With*x*we denote the geometrical length of the optical path between the emitter and receiver antenna and by ℓ the thickness of the investigated sample. The frequency dependent initial electric field leaving the emitter antenna is denoted by*E*._{init}

#### 3.3. The analytical formula of the transfer function

Figure 1 shows the ray propagation through the investigated sample in seven different stages. First, the initial electric field *E _{init}* impinges on the front (left) surface of the sample. In the second stage, we have depicted the two rays that emerge from

*E*being split at the air-sample interface into a transmitted and a reflected portion. The reflected portion ends with the symbol ⊖, which indicates that this portion of the ray has no influence on

_{init}*E*and thus will be discarded. The transmitted portion is equal to

_{samp}*E*

_{init}T_{01}and propagates through the sample. When it reaches the second (right) surface of the sample, the electric field of the wave is equal to

*E*

_{init}T_{01}

*P*

_{1}(ℓ). The next sample-air interface causes the beam to split again. This process is depicted in the third stage of the figure. The transmitted portion, which is equal to

*E*

_{init}T_{01}

*P*

_{1}(ℓ)

*T*

_{10}is a part of

*E*. This fact is indicated by the symbol ⊕ in the third stage of the figure. Subsequently, the reflected part of the beam leads to a periodical reflection pattern, which we refer to as a

_{samp}*Fabry-Perot reflections*. The three dots at the bottom of Fig. 1 indicate the infinite continuation of this reflection pattern. For simplicity, let

*P*

_{1}denote

*P*

_{1}(ℓ). Also, the propagation of the electromagnetic wave through the air between the two antennas has to be considered. This implies the multiplication by a factor

*P*

_{0}(

*x*− ℓ), which we abbreviate by

*P*0. Summing up all components of the electric field transmitted through the sample yields the following formula:

The integer δ in equation (9) represents the number of Fabry-Perot reflections taken into account. Theoretically, this number is infinitely large. Practically, this number is upper bounded: on the one hand, the energy of the electric field is attenuated exponentially with every pass through the sample and on the other hand, the time window of our measurement is finite. In Section 4.2 we present an algorithm for the optimal choice of δ.

The reference THz pulse propagates through air, thus the following equation holds:

Since *P*
_{0} denotes *P*
_{0}(*x* − ℓ), the equalities (9) and (10) together with (8) yield the following
theoretical transfer function:

## 4. Material parameter extraction

For a structured presentation of our material parameter extraction algorithms, we have divided this section into four subsections. Each subsection contains a stand-alone step of the extraction process. In the subsections 4.1 to 4.3, we present several improvements of the existing extraction algorithms. In 4.4 we introduce the Spatially Variable Moving Average Filter, which leads to highly accurate results.

#### 4.1. Preliminary signal processing

Before proceeding to the actual material parameter extraction, we perform the following preliminary signal processing:

*Linear offset compensation of raw time domain data*We assume the existence of a linear offset in the measured photocurrent, which is an often observed experimental artifact. For its compensation, we subtract the linear offset computed from the mean of several data points before the incidence of the pulse and the mean of the last few data points.

*Determination of the reliable frequency range*The reliable frequency range of our experiment is determined as the frequency band, in which the magnitude of the sample signal is greater than or equal to the noise floor of the experiment. The noise floor is obtained by acquiring data from the spectrometer after blocking the THz beam path, and represents the frequency independent detector noise, as discussed in Ref. [5]. In this fashion, we prevent the misinterpretation of absorption coefficients in regions of the spectrum, where there is a detectable reference signal but the entire radiation is absorbed by the sample. Therefore, our method yields similar results to the analytical technique described in Ref. [5].

*Noise cancellation techniques*In Ref. [4,7, 8], noise cancellation techniques such as wavelet de-noising or Wiener filter are presented for a similar experimental setup. However, such techniques may alter the signal shapes in undesired ways. For example, using different wavelets for de-noising leads to different material parameter curves of the same sample. For a comprehensive description of the wavelet de-noising technique by soft thresholding, consult Ref. [9]. The optimal wavelet for de-noising depends on the THz pulse shape, which may vary for different antennas. Using sub-optimal wavelets may lead to suppression of important peaks in the material parameter curves. Therefore, we discard noise cancellation techniques in our work. Instead, we introduce the novel SVMAF technique described in Section 4.4. While respecting the investigated material’s physical properties, our iterative technique suppresses noise on the one hand and statistical errors on the other hand, yielding highly accurate material parameter curves.

#### 4.2. The core of the extraction algorithm

In contrast to Ref. [3], for the minimization of the error Err expressed in equation (4) with respect to the material parameters *n*
_{1}(ω) and κ_{1}(ω) we use the Nelder-Mead algorithm, see Ref. [10]. A similar procedure has been reported in Ref. [11] for the THz differential TDS method. In the following, we derive the initial values which we use for the minimum search algorithm. By measuring the sample thickness with a micrometer screw, we obtain an initial value for the variable ℓ. Furthermore, we assume that the attenuation of the maximum of the magnitude of the THz pulse in the sample is given by the following relation:

where |*E _{s,max}*| and |

*E*| represent the maximal magnitudes of the sample and reference signals, respectively, in the frequency domain. From equation (12) we obtain the following initial value for κ

_{r,max}_{1}:

An initial value for the real refractive index is also easily found. While passing the sample, the maximum of the magnitude of the reference signal in the time domain is delayed by the time Δ*t* due to a difference of the real refractive index Δ*n* = *n*
_{1} − *n*
_{0}. This delay can be expressed by the following equation:

From equation (14) we obtain the following initial value for *n*
_{1}:

In order to proceed with the minimum search, we only need to determine δ in equation (11). For this purpose, we use the rough approximation of *n*
_{1} from equation (15). Let *t _{max}* be the length of the time window measured from the incidence of the reference pulse until the end of the measurement. The integer δ indicates the number of copies of the main pulse that fit in the considered time window of length

*t*. Inside the sample, the wave propagates with the speed

_{max}*c*

_{0}/

*n*

_{1}along a path of length ℓ. Thus, We determine δ as the greatest integer, such that the following inequality is satisfied:

Note, that the factor 2 arises from the fact that each reflection at the sample’s exit interface traverses the sample back and forth before reaching the detector.

The possibility of setting the termination threshold of the Nelder-Mead search allows us to find values for the material parameters, such that *H _{theory}* is arbitrarily close to

*H*. A typical value we used for this threshold is 10

_{experiment}^{-4}.

#### 4.3. Accurate determination of the sample thickness

We employ the total variation technique described in Ref. [3], but, in contrast to Ref. [3], we combine the total variation technique with the Nelder-Mead search. The multiple Fabry-Perot reflections of the ultrashort THz pulse cause an oscillation of the sample signal and, thus, of the transfer function, in the frequency domain. A closer look at equations (5) to (8) and (11) reveals that this oscillation is determined by the sample thickness ℓ and the complex refractive index *ñ*_{1}(ω). Thus, while fitting *H _{theory}* to

*H*by searching the minimum of the error

_{experiment}*Err*, any uncertainty of

*H*caused by an uncertainty of ℓ has to be compensated by the search variables

_{theory}*n*

_{1}(ω) and κ

^{1}(ω). In other words, an uncertainty of ℓ is noticeable in the extracted material parameter curves

*n*

_{1}(ω) and κ

_{1}(ω) as an oscillation due to the Fabry-Perot effect.

We extract material parameter curves for a range of values of the sample thickness ℓ. By implementing a measure of smoothness for these curves, the smoothest curve pair can be chosen from the set of all computed curve pairs. Its corresponding sample thickness indicates the value, which is the closest one to reality among the values in the considered range. In order to define a measure of smoothness we consider the set S of all evaluation points in the frequency domain and we let the integer *m* run through *S*. For a given sample thickness ℓ, the total variation of degree one *TV*(ℓ) is given by the equations:

In other words, determining the sample thickness is equivalent to choosing ℓ such that *TV*(ℓ) is minimal. As mentioned above we use the Nelder-Mead algorithm to find the minimum of the error function *Err* as described in the previous section. This algorithm provides the material parameters, such that *H _{theory}* is arbitrarily close to

*H*. Therefore, uncertainties of

_{experiment}*H*induced by uncertainties of ℓ have no other possibility to be compensated, but by oscillations of the material parameter curves. For this reason, the total variation technique combined with the Nelder-Mead minimum search yield very accurate results for the sample thickness determination.

_{theory}#### 4.4. Accuracy improvement with the Spatially Variant Moving Average Filter

### 4.4.1. The principle

Using the algorithms described in Sections 4.1 to 4.3 for a single pair of reference and sample measurements, we are able to compute a transfer function *H _{theory}*(ω) that lies as close as desired to

*H*(ω). However, the values of

_{experiment}*H*(ω) vary slightly for each measurement of the same sample. Hence, one has to accept a confidence interval for the transfer functions. The confidence interval can be modeled on the basis of a set of measurements.

_{experiment}We use this confidence interval for setting additional physical boundary conditions to the extracted material parameter curves. We suppose that a material parameter curve is accurate if it satisfies the additional physical boundaries and if the modeled transfer function *H _{theory}*(ω) lies within the accepted confidence interval of

*H*(ω).

_{experiment}### 4.4.2. The additional physical boundary condition

Several factors, such as neglected scattering processes, the inhomogeneity of the geometry of the sample and additive noise, lead to a rippled shape of the material parameter curves. The additional physical boundary condition that we set to our model is that the material parameters should not vary strongly at consecutive frequency values, i.e., the curves should have a “smooth” shape rather than a rippled one. As we will see in the discussion of the method, this condition does not affect the detection of absorption peaks which investigated samples may exhibit.

### 4.4.3. The confidence interval

In the following, we assume that several measurements of the reference and sample signals are available. For each considered frequency, we compute an average reference signal $\overline{\stackrel{}{{E}_{\mathrm{ref}}^{\mathrm{ex}}}}(\omega )$having the magnitude and angle equal to the arithmetical means of the magnitudes and of the angles of all reference signals, respectively. Similarly, we obtain the average sample signal $\overline{{E}_{\mathrm{samp}}^{\mathrm{ex}}}(\omega )$. For simplicity, we introduce the following abbreviations:

*H*≔*H*(ω), such that:_{experiment}- $a\u2254\Re \stackrel{\xaf}{{\{E}_{\mathrm{samp}}^{\mathrm{ex}}}(\omega )\}$, $b\u2254\Im \stackrel{\xaf}{{\{E}_{\mathrm{samp}}^{\mathrm{ex}}}(\omega )\}$such that:
- $c\u2254\Re \stackrel{\xaf}{{\{E}_{\mathrm{ref}}^{\mathrm{ex}}}(\omega )\}$,$d\u2254\Im \overline{{\{E}_{\mathrm{ref}}^{\mathrm{ex}}}(\omega )\}$such that:

The goal of this subsection is to determine the frequency dependent confidence intervals Δℜ{*H*} and Δℑ{*H*}, such that for each frequency, every pair of values ℜ{*H _{theory}*(ω)} and ℑ{

*H*(ω) of the modeled transfer function is acceptable if it lies in the intervals [ℜ{

_{theory}*H*} − Δℜ{

*H*}, ℜ{

*H*} + Δℜ{

*H*}] and [ℑ{

*H*] − Δℑ{

*H*}, ℑ{

*H*} + Δℑ{

*H*}], respectively. Note, that an equivalent procedure can be derived using the magnitude and angle instead of the real and imaginary parts of the transfer function.

From the equalities (19) to (21) we derive the following frequency dependent equation:

The quantities *a*,*b*,*c* and *d* are measured values. They are affected by additive noise. In order to characterize the detector noise of the experimental setup, we performed several measurements without THz radiation, i.e. by blocking the THz beam, which have shown that the noise level is frequency independent. For this reason, we assume that the measured quantities a,b,c and d are affected by additive white noise. Furthermore, we observe variations among different measurements of the same sample materials, that are not correlated to the detector noise. The reason for these variations are fluctuations of the THz signal, as mentioned in Ref. [5]. Laser intensity instabilities or slight changes in the electromagnetic environment lead to such variations. We summarize these effects under the class of statistical errors and we assume that they are independent from the error due to the portion of additive white noise discussed above. In conclusion, we assume that the quantities *a*,*b*,*c* and *d* are affected by two independent error sources: *additive white noise* and a *statistical error*.

In the following, we derive expressions for the intervals Δ*x*, with *x* ∈ {*a*,*b*,*c*,*d*}, considering both error sources. Subsequently, we use these intervals to derive Δℜ{*H*} and Δℑ{*H*}.

Let Δ_{N}
*x*, with *x* ∈ {*a*,*b*,*c*,*d*}, be the confidence interval of the variable *x* due to white noise. Note, that Δ* _{N}x* is frequency independent because of the assumption of white noise. In order to compute Δ

*, we consider a noise array*

_{N}x*N*that we obtain from the array

*x*by isolating some last points, such that no signal is detectible in this array. For this purpose, a single pair of reference-sample measurements suffices. Let

*k*be the length of the array

*N*. First, we compute the mean of the noise array:

We estimate the confidence interval of the variable *x* due to white noise as the standard deviation of the noise array:

Optionally, one may define Δ* _{N}x* in a different manner. For example, choosing the variance, i.e. the squared standard deviation, may lead to a narrower confidence interval. However, the standard deviation approach delivers a confidence interval Δ

*having the same physical measuring unit as*

_{N}x*x*.

Let Δ* _{S}x*, with

*x*∈ {

*a*,

*b*,

*c*,

*d*} be the confidence interval of the variable

*x*due to the statistical error. Let

*u*be the number of measurements available for the variable

*x*. If only one reference or sample signal is available, no assertion can be made about the statistical error. In this case, we set Δ

*= 0. If two or more data sets are available, we denote these arrays by*

_{S}x*x*

^{(1)},

*x*

^{(2)}, ⋯,

*x*

^{(u)}. We compute the mean

*x*

^{(mean)}of all measurements:

We estimate the confidence interval of the variable *x* due to the statistical error as the standard deviation of the measurements:

Note that the expressions in (25) and (26) have to be computed for every frequency separately. Since the two error sources are statistically independent, we assume that the total confidence interval is equal to:

As already discussed above, Δ* _{N}x* is frequency independent in contrast to Δ

*, which is frequency dependent.*

_{S}xNow we can turn our attention to determining the confidence intervals Δℜ{*H*} and Δℑ{*H*}. From equation (22) we get:

Since Δℜ{*H*} and Δℑ{*H*} are functions of *a*,*b*,*c* and *d*, the Gaussian error propagation rule yields:

where *f* and *g* are the functions defined in equations (28) and (29), respectively. The confidence intervals Δ*x* with *x* ∈ {*a*,*b*,*c*,*d*} are given by equation (27). For simplicity, we discard the covariance between the variables *a* and *b*, which in principle are correlated. According to Ref. [12], this is equivalent to neglecting the uncertainty of one of them. The same holds for the variables *c* and *d*. However, this approximation leads to a narrower confidence interval. For this reason, the physical boundary condition which we pose here is stronger than if considering the covariances. The partial derivatives of *f* and *g* can easily be calculated from equations (28) and (29):

### 4.4.4. The algorithm

We introduce the Spatially Moving Average Filter (SVMAF) as an extension of the moving average filter, which is also known as adjacent averaging. For material parameter extraction, we have implemented an iterative algorithm using the SVMAF. It provides satisfactory results usually after up to 5 iterations. After describing in this section the actions performed in each iteration step, we illustrate the mode of operation of the algorithm by means of a graphical model in the subsequent section.

- For every point of the material parameter curves, it replaces the initial value of this point with the mean of the initial values of itself and its two neighbors. This is the reason for the name
*moving average*. The first and the last point of the curves keep their initial values. - It computes the transfer function
*H*with the new values of the material parameters._{theory} - For every point, it checks whether the real and imaginary parts of
*H*(ω) lie inside the confidence intervals or not. If both parts lie in the respective ranges, the new values of the material parameters are accepted. If not, the material parameters at the respective frequency are reset to the initial values. Since the action of the filter depends on every point, we called it_{theory}*spatially variant*. - The new material parameter curves represent the input to the next iteration step.

### 4.4.5. Visualization

Figure 2 depicts an example of one iteration step using the SVMAF. The upper half of the figure, i.e. the pictures 1 and 3, shows the progression of the real refractive index *n* while the lower half, i.e. the pictures 2 and 4, shows the real part of the transfer functions and its confidence interval. For simplicity, we only consider one of the material parameters as well as only the real part of the transfer function in this example. In the first pair of figures we see the initial state: the material parameter curve plotted with a continuous line connecting full dots, represents the initial values of *n*(ω). The transfer function plotted with a continuous line always represents ℜ{*H _{experiment}*(ω)}. However, at the beginning of the iteration, the function ℜ{

*H*(ω)} is very close to ℜ{

_{theory}*H*(ω)} since the Nelder-Mead search yields very accurate results. The dashed material parameter curve connecting empty dots, represents the result of the smoothing process in the first step of the algorithm described above and the dashed transfer function represents ℜ{

_{experiment}*H*(ω)} computed with the new material parameters. We see that the value of ℜ{

_{theory}*H*(ω)} at the first point lies outside the confidence interval. Thus, the value of n at this point is reset to its initial value. All other values of ℜ{

_{theory}*H*(ω)} lie inside the confidence interval. Therefore, for these points, the new values are being accepted. The curve consisting of the initial value of the first point and the new values of all other points is the continuous curve in the second plot of

_{theory}*n*. It represents the input to the next iteration step.

When computing the new averaged values, depicted dashed on the second *n*-plot, we observe that the value of the first point is closer to the initial value than its averaged value on the first plot. This is due to the updated values of its neighboring points. Also, we observe that the new averaged value causes the updated ℜ{*H _{theory}*(ω)} to entirely lie within the confidence interval. The three dots on the right hand side of the figure indicate that this procedure can be repeated for an arbitrary number of times.

### 4.4.6. Properties of the SVMAF

Because of the repeated averaging process, the iteratively obtained theoretical material parameter curves are smoother than the initial theoretical curves. Within the experimental error modeled by the confidence intervals, they agree with the experimental curves. The smoothing of peaks in the material parameter curves that physically make sense, modifies the value of {*H _{theory}*(ω)} at the respective point in an artificial way such that its value exceeds the allowed tolerances for ℜ{

*H*(ω)} and/or ℑ{

_{theory}*H*(ω)}. Thus, the SVMAF does not alter such peaks, which is a striking advantage of the SVMAF over common noise cancellation techniques.

_{theory}## 5. Experimental results

In this section we give some experimental results that exemplify the performance of the material parameter algorithm presented in this paper. For each parameter extraction, we use two reference and two sample measurements. Note that using more measurements leads to a better description of the statistical error and, thus, of the confidence interval. However, the confidence interval is also defined for a single reference and a single sample measurement. In this case, it describes only the additive white noise portion, since no assertion about the statistical error can be made. We use a micrometer screw to measure an initial values for the sample thickness determination algorithm described in Section 4.3. We assume an error of the micrometer screw of ±10μm for thick samples. For this purpose we evaluate the total variation of the material parameter curves for a thickness range around the measured value, in which the difference between two consecutive values is 1μm.

#### 5.1. High resistivity silicon wafer

In our first example, we extract the material parameters of a silicon wafer. The thickness measured with the micrometer screw is 539μm. Figure 3(a) shows the temporal profiles of the time domain reference and sample signals. Figure 3(b) depicts the total variation over the sample thickness. The minimum of this curve is located at 543μm, which is the value that we use for the material parameter extraction. The material parameters are depicted in Figures 3(c) and 3(d). The refractive index is in excellent agreement with the one published in Refs. [6] and [13]. Please note that the sample used in Ref. [6] was over 20mm thick, while our sample is nearly 40 times thinner. We observe that the SVMAF leads to much smoother curves than the curves obtained before the SVMAF was applied.

#### 5.2. Pellet consisting of 25mg leu-gly-gly tripeptide and 125mg polyethylene powder

Both material parameter curves of the Si wafer presented in the first example exhibit almost constant shapes. The smoother curves gained due to the SVMAF may lead to the false impression that the SVMAF algorithm only averages neighboring points. Our second example should exclude such an impression by examining a material that shows a clear frequency dependency of the material parameter curves. We investigate a pellet consisting of 25mg leu-gly-gly tripep-tide and 125mg polyethylene powder. The material is pressed into a pellet under 380MPa for 10 minutes.

The time domain signals are depicted in Figure 4(a). Even though the Fabry-Perot effect is almost invisible to the eye in the time domain curves, the total variation exhibits a clear minimum at 1252μm, as shown in Figure 4(b). This is due to the very accurate Nelder-Mead algorithm used to determine the error between the modeled and the experimental transfer functions. This fact points out the applicability of our thickness extraction algorithm even in cases with a low Fabry-Perot effect. Figures 4(c) and 4(d) show the material parameters. We observe the strong frequency dependency of both material parameters. Note that the SVMAF algorithm does not affect the peaks of the curves that make sense physically. This is one of the major qualities of the SVMAF technique.

#### 5.3. Thin artificial RNA pressing

In the following we investigate a pressing of 22mg of artificial RNA powder. Our goal in this section is to exemplify the performance of our parameter extraction algorithm for very thin samples. The sample under test lies far outside the ranges of experiments documented in several papers, compare to Ref. [3].

With a micrometer screw we measured a thickness of 150μm. However, this measurement may have been inaccurate due to the softness of the the sample. With the total variation technique we determine the actual sample thickness equal to 145μm.

The computed material parameter curves are shown in Figures 5(a) and 5(b). Obviously, the extraction algorithm works very well for this thin sample with a relatively small real refractive index. Since the sample is very thin, the multiple Fabry-Perot reflections overlap in the time domain. Thus, time windowing the first pulse passing the probe is impossible. In Figures 5(c) and 5(d) we show the material parameters extracted with an algorithm that does not consider the overlapping Fabry-Perot reflections. The experimental data, from which the material parameters were extracted are the same as the data presented in Ref. [14]. We can clearly observe artifacts in the region around and below 0.5THz in the curves caused by the Fabry-Perot effect. Furthermore, Figures 5(c) and 5(d) show that towards higher frequencies, the ripple of both curves increases. This happens due to poor SNR at higher frequencies. Figures 5(a) and 5(b) demonstrate that the SVMAF counters this effect.

## 6. Conclusions

There are two main contributions of this work. First, we have improved the existing material parameter extraction techniques for THz TDS by merging the up-to-date knowledge into a single algorithm. Our method provides highly accurate curves of the complex refractive index of a wide range of samples. It works well for samples as thin as 100μm, where multiple reflections overlap, samples with low indexes of refraction and samples with sharp peaks in the material parameter curves. The second main contribution is the introduction of the Spatially Variant Moving Average Filter which constitutes the core of an iterative algorithm yielding the physically most accurate modeled function out of a set of functions that have to satisfy certain boundaries. In our case, we use THz TDS to determine the material parameters of a homogenous sample by fitting a modeled transfer function to a measured one. Measured transfer functions vary for every measurement. We use these variations to define confidence intervals for the transfer functions which allow for the formulation of supplemental physical boundaries to the material parameter curves. Our iterative algorithm converges to a solution within the confidence intervals. We perceive this solution as the most accurate characterization of the optical parameters of the investigated material to date. However, the idea behind this algorithm may be useful for several other experimental techniques, where a modeled function is fitted to an experimentally determined one.

## Acknowledgments

We thank Radu Cernat from the Institute of High Voltage and EMC, University of Dortmund for the support and useful discussion. Rafal Wilk gratefully acknowledges the financial support from “Gottlieb Daimler and Karl Benz” foundation.

## References and links

**1. **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]

**2. **L. Duvillaret, F. Garet, and J.-L. Coutaz, “
Highly precise determination of Optical Constants and sample thickness in Terahertz Time-Domain Spectroscopy,” Appl. Opt. **38**,409 –415 (1999). [CrossRef]

**3. **T. Dorney, R. Baraniuk, and D. Mittleman, “Material parameter estimation with Terahertz Time-Domain Spec-troscopy,” J. Opt. Soc. Am. A **18**,1562 –1571 (2001). [CrossRef]

**4. **B. Ferguson and D. Abbott
, “Signal processing for T-Ray Biosensor Systems,” Smart Electronics and MEMS II **4236**,157 –169 (2001).

**5. **P. U. Jepsen and B. M. Fischer, “Dynamic range in Terahertz Time-Domain Transmission and Reflection Spec-troscopy,” Opt. Lett. **30**,29 –31 (2005). [CrossRef] [PubMed]

**6. **D. Grischkowsky, S. Keiding, M. van Exter, and C. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B **7**,2006 –2015 (1990). [CrossRef]

**7. **B. Ferguson and D. Abbott, “Wavelet de-noising of optical terahertz pulse imaging data,” J. Fluct. Noise Lett. **1**,L65 –L69 (2001). [CrossRef]

**8. **B. Ferguson and D. Abbott, “De-noising techniques for terahertz responses of biological samples,” Microelectronics Journal (Elsevier) **32**,943 –953 (2001).

**9. **D. L. Donoho, “De-noising by Soft-Thresholding,” IEEE Trans. Inf. Theory **41**,613 –627 (1995). [CrossRef]

**10. **J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions,”SIAM J. Optim. **9**,112 –147 (1998). [CrossRef]

**11. **S. Mickan, K.-S. Lee, T.-M. Lu, J. Munch, X.-C. Zhang, and D. Abbott, “Thin Film Characterization using Terahertz Differential Time-Domain Spectroscopy and Double Modulation,” Electronics and Structures for MEMS IIz Proc. SPIE **4591**,197 –209 (2001).

**12. **DIN Deutsches Institut für Normung, “DIN 1319 Fundamentals of metrology - Part 4: Evaluation of measurements; uncertainty of measurement,” Beuth Verlag GmbH (1999).

**13. **M. N. Afsar, “Dielectric Measurements of Millimeter-Wave Materials,” IEEE Trans. Microwave Theory and Techniques **12**,1598 –1609 (1984). [CrossRef]

**14. **B. M. Fischer, M. Hoffmann, R. Wilk, F. Rutz, T. Kleine-Ostmann, M. Koch, and P. Uhd Jepsen, “Terahertz
time-domain spectroscopy and imaging of artificial RNA,” Opt. Express **13**,5205 –5215 (2005). [CrossRef] [PubMed]