OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editor: Gregory W. Faris
  • Vol. 2, Iss. 5 — May. 17, 2007
« Show journal navigation

Highly accurate optical material parameter determination with THz time-domain spectroscopy

Ioachim Pupeza, Rafal Wilk, and Martin Koch  »View Author Affiliations


Optics Express, Vol. 15, Issue 7, pp. 4335-4350 (2007)
http://dx.doi.org/10.1364/OE.15.004335


View Full Text Article

Acrobat PDF (426 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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

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

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

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

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

, 5

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]

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

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]

] and Ref. [3

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]

] 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

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]

]). 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

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

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

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

] 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

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]

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

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]

].

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

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]

]. 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 . The sample signal is obtained by recording the photocurrent of the reference signal, additionally propagated through the sample material.

We shall denote by Erefex (ω) and Eexsamp(ω) 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:

Hexperiment(ω)=Esampex(ω)Erefex(ω).
(1)

Let Htheory(ω) be a theoretically derived transfer function of the sample material, which depends on the material parameters. We will derive an analytical expression for Htheory(ω) 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 Htheory(ω) is as close as possible to Hexperiment(ω). As a measure for how close these two functions are, we define the following error function Err:

M(ω)=Htheory(ω)Hexpriment(ω),
(2)
A(ω)=Htheory(ω)Hexpriment(ω),
(3)
Err=Σω(M(ω)+A(ω))
(4)

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:

    R12=n2˜n1˜n1˜+n2˜,
    (5)

    T12=2n˜1n˜1+n˜2.
    (6)

    Fig. 1. Ray propagation through a planar, homogenous sample.

    Furthermore, the propagation of an electric field Einit entering the material 1 is governed by the following relation:

    E(z)=EinitP1(z),
    (7)

    P1(z)=exp(n˜1zc0),
    (8)

    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.

    Fig. 1.Ray propagation through a planar, homogenous sample.
  • Standard notation used throughout this work

    For computations, we use the value cvacuum = 2.99796 ∙ 108 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 ñ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 Einit.

3.3. The analytical formula of the transfer function

Esamp=EinintP0T01P1T10[1+Σi=1δ(R102P12)i]
(9)

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:

Eref=EinitP0(x)
(10)

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

Htheory=EsampEref=P0()T01P1T10[1+Σi=1δ(R102P12)i]
(11)

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

    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]

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

    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]

    ].

  • Noise cancellation techniques

    In Ref. [4

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

    ,7

    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

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

    ], 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

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

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

Es,max=Er,maxexp(ωc0κ1),
(12)

where |Es,max| and |Er,max| 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 κ1:

κ1=1c0ωlogEs,maxEr,max.
(13)

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 1n 0. This delay can be expressed by the following equation:

Δt=Δnc0.
(14)

From equation (14) we obtain the following initial value for n 1:

n1=c0Δt+n0.
(15)

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 tmax 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 tmax. Inside the sample, the wave propagates with the speed c 0/n 1 along a path of length ℓ. Thus, We determine δ as the greatest integer, such that the following inequality is satisfied:

tmaxn1c0(1+2δ).
(16)

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.

4.3. Accurate determination of the sample thickness

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:

D(m)=n1(m1)n1(m)+κ1(m1)κ1(m),
(17)
TV()=mSD(m).
(18)

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 Htheory(ω) that lies as close as desired to Hexperiment(ω). However, the values of Hexperiment(ω) 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.

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 Htheory(ω) lies within the accepted confidence interval of Hexperiment(ω).

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 Erefex¯(ω)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 Esampex¯(ω). For simplicity, we introduce the following abbreviations:

  • HHexperiment(ω), such that:

    H=Esampex¯(ω)Erefex¯(ω).
    (19)

  • a{Esampex¯ (ω)}, b{Esampex¯ (ω)}such that:

    Esampex¯(ω)=a+jb
    (20)

  • c{Erefex¯ (ω)},d{Erefex¯(ω)}such that:

    Erefex¯(ω)=c+jd.
    (21)

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 ℜ{Htheory(ω)} and ℑ{Htheory(ω) of the modeled transfer function is acceptable if it lies in the intervals [ℜ{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:

H=a+jbc+jd=ac+bd+j(bcad)c2+d2.
(22)

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

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]

]. 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 ΔNx is frequency independent because of the assumption of white noise. In order to compute ΔNx, we consider a noise array 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:

mx=1kΣi=1kNi.
(23)

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

ΔNx=1k1Σi=1k(Nimx)2.
(24)

Optionally, one may define ΔNx 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 ΔNx having the same physical measuring unit as x.

Let ΔSx, 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 ΔSx = 0. If two or more data sets are available, we denote these arrays by x (1), x (2), ⋯, x (u). We compute the mean x (mean) of all measurements:

x(mean)=1uΣi=1ux(i).
(25)

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

ΔSx=1u1Σi=1u(x(i)x(mean))2.
(26)

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:

Δx=(ΔNx)2+(ΔSx)2.
(27)

As already discussed above, ΔNx is frequency independent in contrast to ΔSx, which is frequency dependent.

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

{H}=ac+bdc2+d2f(a,b,c,d),
(28)
{H}=bcadc2+d2g(a,b,c,d),
(29)

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

Δℜ{H}=(Δafa)2+(Δbfb)2+(Δcfc)2+(Δdfd)2,
(30)
Δℑ{H}=(Δaga)2+(Δbgb)2+(Δcgc)2+(Δdgd)2,
(31)

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

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

], 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):

fa=cc2+d2,
(32a)
fb=cc2+d2,
(32b)
fc=ad2ac22bcd(c2+d2)2,
(32c)
fd=bc2bd22acd(c2+d2)2,
(32d)
ga=dc2+d2,
(32e)
gb=cc2+d2,
(32f)
gc=bc2bd22acd(c2+d2)2,
(32g)
gd=ad2ac22bcd(c2+d2)2.
(32h)

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 Htheory with the new values of the material parameters.
  • For every point, it checks whether the real and imaginary parts of Htheory(ω) 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 spatially variant.
  • The new material parameter curves represent the input to the next iteration step.

4.4.5. Visualization

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 ℜ{Htheory(ω)} 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 {Htheory(ω)} at the respective point in an artificial way such that its value exceeds the allowed tolerances for ℜ{Htheory(ω)} and/or ℑ{Htheory(ω)}. Thus, the SVMAF does not alter such peaks, which is a striking advantage of the SVMAF over common noise cancellation techniques.

Fig. 2. Illustration of the mode of operation of the SVMAF

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

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]

] and [13

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

]. Please note that the sample used in Ref. [6

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]

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

Fig. 3. Material parameter extraction of a 543μm Si wafer
Fig. 4. Material parameter extraction of a 1252μm leu-gly-gly pellet
Fig. 5. Material parameter extraction of a 145μm artificial RNA pressing

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

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]

].

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

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]

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

OCIS Codes
(070.6020) Fourier optics and signal processing : Continuous optical signal processing
(300.6340) Spectroscopy : Spectroscopy, infrared

ToC Category:
Spectroscopy

History
Original Manuscript: February 16, 2007
Revised Manuscript: March 19, 2007
Manuscript Accepted: March 21, 2007
Published: April 2, 2007

Virtual Issues
Vol. 2, Iss. 5 Virtual Journal for Biomedical Optics

Citation
Ioachim Pupeza, Rafal Wilk, and Martin Koch, "Highly accurate optical material parameter determination with THz time-domain spectroscopy," Opt. Express 15, 4335-4350 (2007)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-15-7-4335


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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 Spectroscopy," 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 Spectroscopy," 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," Fluct. Noise Lett. 1, L65 - L69 (2001). [CrossRef]
  8. B. Ferguson and D. Abbott, "De-noising techniques for terahertz responses of biological samples," Microelectron. J. (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 II Proc.SPIE 4591, 197 - 209 (2001).
  12. DIN Deutsches Institut f¨ur 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 Tech. 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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited