## Abstract

In rapid scan Fourier transform spectrometry, we show that the noise in the wavelet coefficients resulting from the filter bank decomposition of the complex insertion loss function is linearly related to the noise power in the sample interferogram by a noise amplification factor. By maximizing an objective function composed of the power of the wavelet coefficients divided by the noise amplification factor, optimal feature extraction in the wavelet domain is performed. The performance of a classifier based on the output of a filter bank is shown to be considerably better than that of an Euclidean distance classifier in the original spectral domain. An optimization procedure results in a further improvement of the wavelet classifier. The procedure is suitable for enhancing the contrast or classifying spectra acquired by either continuous wave or THz transient spectrometers as well as for increasing the dynamic range of THz imaging systems.

© 2003 Optical Society of America

## 1. Introduction

In recent years, considerable progress has been made towards the development of imaging systems utilizing THz transient spectrometers [1–5]. It is widely accepted that progress in new application areas such as biomedical imaging [6–15] where regions of high water content dramatically attenuate the signal will require efficient signal processing algorithms to improve the contrast of THz spectra. In an earlier paper [16], we described the advantages of performing filtering, regression and classification of THz spectra using a range of linear transforms. In addition, we also showed that it is possible to optimize the integration time in Fourier transform spectrometry after analyzing on-line the statistical properties of the wavelet transform coefficients [17]. Our work complements that in other laboratories which have been performing filtering in either the time or frequency domain using wavelet transforms [18–21]. In these papers, however, no attempt has been made to describe analytically the propagation of noise as these linear transformations are performed. Moreover, no work has been done on the optimization of the wavelet transform for the treatment of specific THz data sets. These are the subjects of the current paper.

In what follows we show that when the signal in the background interferogram is large compared to the noise, the noise in the sample interferograms, propagates linearly in the frequency domain. We characterize the noise in the wavelet coefficients resulting from the filter bank decomposition of the complex insertion loss after the standard operations of apodization, Fourier transformation and ratioing against the background. We show that the noise power in each wavelet coefficient can be related to the noise power in the sample interferogram by a noise amplification factor. Thus, under the hypothesis that the variance in the noise remains constant across different time instances in the interferogram acquisition process, the signal-to-noise ratio of that wavelet coefficient may be improved by maximizing an objective function composed of the power of the wavelet coefficients divided by the noise amplification factor. The formulation ensures that optimal feature extraction in the wavelet domain is performed. Examples using noisy lycra and leather transmission spectra show that the performance of a classifier based on the output of a single level filter bank using the db4 mother wavelet is considerably better than that of an Euclidean distance classifier in the original spectral domain. Furthermore, the proposed optimization procedure results in a further improvement in the signal-to-noise ratio of the wavelet coefficients used for classification.

## 2. Propagation of noise in the wavelet domain

The analysis starts by assuming that the noise in the background interferogram *b*(*t*) is negligible when compared to the noise in the sample interferograms *x*(*t*). A vector *x*(*t*) of length 2*J*, (*t*=0, 1, …, 2*J*-1), can be assumed to be a stochastic process composed of a signal term *x _{m}*(

*t*) and a zero mean noise term

*n*(

*t*) so that:

where *x _{m}*(

*t*)=

*E*[

*x*(

*t*)],

*E*[

*n*(

*t*)]=0, ∀

*t*,

*E*[

*n*(

*t*

_{1})

*n*(

*t*

_{2})]=0, ∀

*t*

_{1}≠

*t*

_{2}and

*E*[

*n*

^{2}(

*t*)]=

*σ*

^{2}, ∀

*t*with

*E*denoting the expectation operator and

*σ*the noise standard deviation. For an apodization window

*w*(

*t*), the apodized sample interferogram

*x*(

_{a}*t*) is

*x*(

_{a}*t*)=

*x*(

*t*)

*w*(

*t*) which can also be divided into signal and noise terms as:

The Discrete Fourier Transform (DFT) *X _{a}*(

*ω*) of the apodized interferogram can be written as:

where

The apodized noise in terms of its discrete Fourier transform can be written as ${N}_{a}\left(\omega \right)=\sum _{t}n\left(t\right)w\left(t\right){e}^{-j\omega t}$ and thus

In addition, the autocorrelation of the apodized noise in the frequency domain is:

$$=E\left[\left(\sum _{{t}_{1}\ne {t}_{2}}n\left({t}_{1}\right)w\left({t}_{1}\right){e}^{-j{\omega}_{1}{t}_{1}}n\left({t}_{2}\right)w\left({t}_{2}\right){e}^{+j{\omega}_{2}{t}_{2}}\right)+\left(\sum _{t}{n}^{2}\left(t\right){w}^{2}\left(t\right){e}^{-j\left({\omega}_{1}-{\omega}_{2}\right)t}\right)\right]$$

$$=\left(\sum _{{t}_{1}\ne {t}_{2}}\underset{0}{\underbrace{E\left[n\left({t}_{1}\right)n\left({t}_{2}\right)\right]}}w\left({t}_{1}\right){e}^{-j{\omega}_{1}{t}_{1}}w\left({t}_{2}\right){e}^{+j{\phantom{\rule{.2em}{0ex}}\omega}_{2}{t}_{2}}\right)+\left(\sum _{t}\underset{{\sigma}^{2}}{\underbrace{E\left[{n}^{2}\left(t\right)\right]}}{w}^{2}\left(t\right){e}^{-j\left({\omega}_{1}-{\omega}_{2}\right)t}\right)$$

where * denotes complex-conjugate. By letting *z*(*t*)=*w*
^{2}(*t*), it follows that:

where *Z*(*ω*) is the DFT of *z*(*t*). Let *L*(*ω*) be the complex insertion loss of the sample, that is,
$L\left(\omega \right)=\frac{{X}_{a}\left(\omega \right)}{{B}_{a}\left(\omega \right)}$
where *B _{a}*(

*ω*) is the DFT of the apodized background interferogram. From Eq. (3) it follows that:

where *N _{ab}*(

*ω*) is the apodized noise in the background. Under the hypothesis that this noise in the background is negligible compared to the background signal, 8a becomes:

Setting
${L}_{m}\left(\omega \right)=\frac{{X}_{\mathit{ma}}\left(\omega \right)}{{B}_{a}\left(\omega \right)}$
and
$M\left(\omega \right)=\frac{{N}_{a}\left(\omega \right)}{{B}_{a}\left(\omega \right)}$
it follows that *L*(*ω*)=*L _{m}*(

*ω*)+

*M*(

*ω*). The noise term

*M*(

*ω*) has the following properties:

and

## 3. Wavelet transform of the complex insertion loss function

The wavelet transform of the discrete signal **L**=[*L*(*ω*
_{0}) *L*(_{ω}
*1*) … *L*(*ω*
_{J-1})], obtained after eliminating the right half of the DFT result, can be written as:

where *ω _{n}* is the

*n*

^{th}frequency point in the DFT and the basis elements

*ψ*(wavelets) are built from a mother wavelet function

_{a,b}*ψ*as:

Subscript *a*∊R* is the “Scale” (or “Dilation”) which determines the width of the wavelet and subscript *b*∊R defines the position of the wavelet with respect to signal **L**. Equations (11) and (12) define a continuous wavelet transform, in the sense that the real-valued parameters *a* and *b* are allowed to vary continuously. The discrete wavelet transform is obtained by taking parameters *a* and *b* from a dyadic grid defined in the following manner:

where *s* and *r* are integer numbers. In this case, the wavelet transform results in a set of coefficients indexed by *s* (scale level) and *r* (translation index). These coefficients can be obtained in a computationally efficient manner by using the tree algorithm [22,23] depicted in Fig. 1. In this scheme, **c _{0}** is the result of convolving

**L**with a “scaling function”, which is the low-pass counterpart of the mother wavelet. This can be regarded as a low-pass pre-filtering procedure. The wavelet transform coefficients at scale level

*s*are stored in sequence

**d**

_{s}.

In Fig. 1, *H* and *G* represent respectively low-pass and high-pass filters associated with the mother wavelet adopted. Symbol (↓2) denotes the down-sampling operation, which consists of eliminating every other coefficient of a sequence. Down-sampling ensures that the number of points remains the same after the wavelet transform. In each scale level s two sets of coefficients are generated: **d _{s}** (which are also called “detail” coefficients in this context) and

**c**(“approximation” coefficients). The tree algorithm is iterative in the sense that coefficients

_{s}**c**

_{s}are further decomposed in order to generate

**c**and

_{s+1}**d**. The filtering and down-sampling operations can be summarized in the following equations, in which sequences

_{s+1}*h*(

*i*) and

*g*(

*i*), of 2

*Q*points each (

*i*=0, 1,…2

*Q*-1), are the impulse response functions of the low-pass and high-pass filters

*H*and

*G*respectively.

(where lowercase *i* denotes the index of elements in the filtering sequence, whilst 2*Q* denotes the number of points in the sequence). It is worth noting that most works use the data vector itself as the input to the tree algorithm (instead of **c _{0}**). This results in the

**d**coefficients obtained being actually approximations to the wavelet transform coefficients. This approximation improves as

_{s}*s*increases. However, in most cases, it may not be essential to obtain the actual wavelet transform coefficients, but just to generate coefficients in which the energy of the signal is compressed. In this case, both

**c**and

_{s}**d**coefficients can be regarded as the result of an energy compressing linear transform and can thus be used for filtering and feature extraction. In this filter bank, the low-pass filtering result undergoes successive filtering iterations with the number of iterations

_{s}*S*(

*S*<log

_{2}(

*J*)) chosen by the analyst. The final result of the decomposition of data vector

**L**is a vector

**p**resulting from the concatenation of row vectors

**c**

_{s}and

**d**

_{s}in the following manner:

with coefficients in larger scales (e.g. **d**
_{S}, **d**
_{S-1}, **d**
_{S-2}) associated with broad features in the data vector, and coefficients in smaller scales (e.g. **d**
_{1}, **d**
_{2}
**d**
_{3}…) associated with narrower features such as sharp peaks. In this manner, the double index (*s,r*) of the filter bank result can be replaced with a single index *k* indicating the position of the coefficient in vector **p**.

If the impulse responses {*h*
_{0}, *h*
_{1}, …, *h*
_{2Q-1}} and {*g*
_{0}, *g*
_{1}, …, *g*
_{2Q-1}} satisfy the following conditions:

with *l* an arbitrary point in the length of the sequence, the structure in Fig. 1 is termed a quadrature-mirror filter (QMF) bank [22–24]. A QMF bank is said to enjoy a perfect reconstruction (PR) property, because vector **L** can be reconstructed from vector **p** which means that there is no loss of information in the decomposition process. Moreover if the convolution operation is performed in a circular manner [25], after data are periodized to minimize border effects, the transformation from **L** to **p** can be represented by a matrix operation **p**
_{1×J}=**L**
_{1×J}
**V**
_{J×J} where **V** is an orthogonal matrix, i.e. **V**
^{T}
**V**=**I**. For a single decomposition level, for instance, the **V** matrix is given as **V**=[* H*|

*], where*

**G***and*

**H***are sub-matrices built from the low-pass and high-pass filter impulse responses as*

**G**If more decomposition levels are used, the left hand side of matrix **V** will include products between the low-pass and high-pass filter coefficients due to the successive splitting of the low-pass channel. It is worth noting that the transform is entirely defined by the low-pass filter response, because the low-pass and high-pass filters are related by Eq. (16b). Since there are *Q* restrictions on the 2*Q* coefficients of the filter response, there remain *Q* degrees of freedom that can be exploited to adjust the transform to a particular data-set. Given a suitable objective function, these adjustments of the transform can be formulated as a non-linear constrained optimisation problem. To circumvent the difficulty caused by the restrictions we use the formulation proposed by Vaidyanathan [24], as adopted by Sherlock and Monro [26] to parametrise orthonormal wavelets of compact support, by *Q* angles {*θ*
_{0}, *θ*
_{1}, …, *θ*
_{Q-1}} in a lattice structure as shown in Fig. 2. The problem then becomes one of unconstrained optimization in *R ^{Q}* which can be solved by standard numerical search procedures.

Now, let *p*(*k*) be the *k*
^{th} coefficient resulting from the decomposition of **L** by the filter bank shown in Fig. 1, that is,

where *v _{k}*(

*n*) is the

*n*

^{th}element of the

*k*

^{th}column of the transformation matrix

**V**associated with the filter bank. Since

*L*(

*ω*)=

_{n}*L*(

_{m}*ω*)+

_{n}*M*(

*ω*), it follows that:

_{n}Let ${p}_{m}\left(k\right)=\sum _{n}{L}_{m}\left({\omega}_{n}\right){v}_{k}\left(n\right)\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\mathrm{and}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}f\left(k\right)=\sum _{n}M\left({\omega}_{n}\right){v}_{k}\left(n\right)\phantom{\rule{.2em}{0ex}}\mathrm{so}\phantom{\rule{.2em}{0ex}}\mathrm{that}\phantom{\rule{.2em}{0ex}}p\left(k\right)={p}_{m}\left(k\right)+f\left(k\right)$ .

Assuming that the coefficients of high-pass and low-pass filtering sequences (and thus the elements of **V**) are real-valued, the noise term *f*(*k*) has the following properties:

$$=E\left[\sum _{{n}_{1}}\sum _{{n}_{2}}M\left({\omega}_{{n}_{1}}\right){M}^{*}\left({\omega}_{{n}_{2}}\right){v}_{k}\left({n}_{1}\right){v}_{k}\left({n}_{2}\right)\right]=\sum _{{n}_{1}}\sum _{{n}_{2}}E\left[M\left({\omega}_{{n}_{1}}\right){M}^{*}\left({\omega}_{{n}_{2}}\right)\right]{v}_{k}\left({n}_{1}\right){v}_{k}\left({n}_{2}\right)$$

$$={\sigma}^{2}\underset{\mathit{NAF}\left(k\right)}{\underbrace{\sum _{{n}_{1}}\sum _{{n}_{2}}\frac{Z\left({\omega}_{{n}_{1}}-{\omega}_{{n}_{2}}\right)}{{B}_{a}\left({\omega}_{{n}_{1}}\right){{B}_{a}}^{*}\left({\omega}_{{n}_{2}}\right)}{v}_{k}\left({n}_{1}\right){v}_{k}\left({n}_{2}\right)}}={\sigma}^{2}\mathit{NAF}\left(k\right)$$

where the evaluation of the double sum is performed noting that *Z*(0)=∑_{t}*w*
^{2}(*t*) and *Z*(-*ω*)=*Z**(*ω*). It follows that the noise power in the *k*
^{th} wavelet coefficient is related to the noise power in the sample interferogram by a noise amplification factor denoted by NAF(*k*). Thus, under the hypothesis that *σ*
^{2} is constant, the signal-to-noise ratio of that wavelet coefficient is improved by maximizing the following objective function relating to the energy in each particular wavelet coefficient divided by a factor that measures the importance in the noise propagation:

The objective, therefore, is to maximize the energy of the wavelet coefficient and minimize the effect of the noise. This is achieved after a variable selection procedure which eliminates some of the wavelet coefficients (thresholding). The optimization is carried out separately for each coefficient as the objective function is related to one coefficient at a time. The flexible polyhedron algorithm [27] available in the Matlab Optimization Toolbox was employed to search for the optimal values for the lattice angles *θ*, using the parameters associated with the db4 wavelet as a starting point.

Moreover, if data sequence *x*(*t*) is taken as the difference between two sample interferograms, the wavelet coefficients *p*(*k*), *k*=0, 1, …, *J*-1 can be used as discriminant values. In this case, *Ξ*(*k*) will reflect the signal-to-noise ratio in the *k*
^{th} discriminant coefficient. Wavelet coefficients with low signal-to-noise ratio are discarded. The maximization of the objective function *Ξ*(*k*) circumvents the main difficulty of the wavelet transform which is how to choose the best mother wavelet for a particular application.

Figure 3 illustrates the operations involved in the optimization process described above for the discrimination of interferograms at four particular time instants *t*
_{1}, *t*
_{2}, *t*
_{3}, *t*
_{4}. Each interferogram can be regarded as a point in a four-dimensional space. After the Fourier transformation and upon elimination of the right half of the resulting vectors, the patterns become spectra with real and imaginary parts at two frequencies *ω*
_{1} and *ω*
_{2}. In Figs. 3a and 3b, the two patterns to be discriminated are represented as points in two complex planes, each plane associated with one of the frequencies. The larger circles represent the noise associated with the patterns, which is equal in all directions because the Fourier transform is orthogonal. The difference Δ*X _{a}* is represented in each plane (Figs. 3(c) and (d)) as a vector joining the two points. Ratioing against the background at each frequency is equivalent to rotating and contracting the difference vectors, which then become differences in complex insertion loss Δ

*L*. Since the background intensity is different at each frequency, the noise level becomes different in Δ

*L*(

*ω*

_{1}) and Δ

*L*(

*ω*

_{2}). Since, the noise is no-longer equal in all directions, as shown in Figs. 3(e) and (f), the Euclidean distance classifier is no-longer optimal [16]. The optimization of the objective function

*Ξ*(

*k*) performs a further rotation of the four axes trying to find a single direction for the

*k*

^{th}wavelet coefficient in which the projection OB of the difference is maximal with respect to the projection AC of the noise.

## 4. Classification of spectra in the wavelet domain

To illustrate the above procedure, Fig. 4 show typical spectra of lycra and leather samples of thickness used by the clothing industry as measured in transmission in the 100 GHz to 1.2 THz frequency range using a Fourier transform spectrometer [10]. Noise up to 0.15 THz is due to the low source power and the restricted aperture of the detector, whereas the arrows and the low background spectrum above 1.2 THz are due to water vapour lines as the spectrometer is not evacuated.

Figure 5 presents the real and imaginary parts of the continuous db4 wavelet transform as contour plots for the leather and lycra interferograms, as well as for the difference between them. At small scales, the effect of the noise at frequencies up to 0.15 THz can be easily separated from the other spectral regions, because each wavelet coefficient is only associated with a narrow frequency region. However, as the scale increases, the wavelet coefficients become associated with wider frequency regions, and the low-frequency noise becomes mixed with the information from other parts of the spectrum. This phenomenon can be better visualized in Fig. 5, by noting that the low-frequency noise displays a cone of influence that encompasses frequencies up to 0.6 THz as the scale increases to 256. A similar phenomenon can be seen for the noise at frequencies larger than 1.2 THz. This reasoning shows that in situations where noise is concentrated in parts of the signal, the use of very large scales in the wavelet transform is not recommended. Moreover, in such situations, linear transforms that do not preserve spatial information, such as principal component analysis (PCA) [16], may not provide an adequate separation of signal from the noise. This is the reason why PCA has been used in conjunction with statistical techniques that conveniently select certain spectral regions [28].

Figure 6 presents the results of carrying out the wavelet processing stage with a single-level filter bank using the db4 mother wavelet (attached Matlab file). Figure 6c shows the objective function values for the coefficients resulting from the filter bank decomposition. As can be seen, the detail coefficients (red region) have a very small *Ξ* value, which means that they are basically associated with the noise. Moreover, by comparing Fig. 6a and 6c, it can be seen that large coefficients may not have a large *Ξ* value, since their noise propagation factor may also be large. In this sense, a deterministic approach to the selection of wavelet coefficients, which only takes into account the magnitude of the coefficients (Fig. 6(a)) without considering the noise, might not yield good results. Function *Ξ* (Fig. 6(c)) can be employed to select the wavelet coefficients to include in the classification model. Statistical criteria [29] could be used to weight the discriminative power of the model against its complexity, in order to determine the best number of wavelet coefficients to be used. Without loss of generality, the 5 wavelets with the largest *Ξ*-value will be employed in a simulated study to illustrate the utility of the proposed optimization procedure.

The optimization was carried out separately for each of the 5 coefficients. As a result, 5 different sets of filter weights were obtained. The inset in Fig. 7 compares the *Ξ*-values before and after the optimization. As can be seen, a considerable improvement was obtained in the objective function for all coefficients. For illustration purposes, Fig. 7 depicts the effect of the optimization on the low-pass filter weights *h*(*i*) for the wavelet coefficient that had the largest *Ξ*-value. The original db4 filter weights are represented as circles, whereas the filter weights resulting from the optimization are represented as squares. The thick vertical lines are used to help the visualization of the magnitude of the weights. This graph shows that the filters were significantly modified by the optimization algorithm.

A classification test was carried out by adding an artificial zero-mean gaussian noise to the interferograms, in order to simulate signals acquired with a smaller integration time. Each noisy pattern was apodized, Fourier-transformed, ratioed against the background and then wavelet-transformed. The five wavelet coefficients selected in the optimization process were then used as classification variables. These are taken from regions where the background signal to noise ratio is highest. The classification was accomplished by analyzing the norm of the difference between the wavelet coefficients of the noisy pattern and the corresponding wavelet coefficients of the original lycra and leather patterns. For comparison, the classification was also carried out in the original spectral variable domain, by using a standard Euclidian distance classifier.

Figure 8 presents the result of this simulated test, in terms of classification errors. To generate this graph, the standard deviation of the noise was varied from 0.001 to 0.5. For each noise level, 250 noisy patterns were generated for each class (lycra and leather). As can be seen, the classification is much more robust to noise when carried out in the wavelet domain than in the original domain. Moreover, the robustness to noise is further increased by the optimization of the wavelet transform.

## 5. Conclusion

This paper presented a detailed description of the propagation of noise from time-domain interferograms to the wavelet transform of the complex insertion loss function after taking into account the apodization function used in the time domain interferograms. This description was used to develop a wavelet coefficient selection criterion based on a function that reflected the signal-to-noise ratio of the coefficients (Fig. 6(c)). The use of this function allows the noise propagation to be taken into account in the wavelet coefficient selection process, unlike standard thresholding techniques based only on the magnitude of the wavelet coefficients (Fig. 6(a)). As a result, it is possible to discard large coefficients that are associated with complex insertion loss peaks caused by ratioing against a small background intensity, a typical situation arising in certain regions of the spectrum when the spectrometer is not evacuated. An optimization procedure for increasing the signal-to-noise ratio of selected wavelet coefficients was also proposed.

The utility of the proposed approach was illustrated in a simulation study, in which noisy lycra and leather interferograms were to be discriminated by comparison with standards with a smaller noise level. The results showed that the classification is much more robust to noise when it is carried out on the basis of a selected number of wavelet coefficients instead of the original spectral variables. In addition, the proposed optimization procedure, which increases the signal-to-noise ratio of the selected wavelet coefficients, further improved the classification performance for a fixed noise level. Conversely, to achieve a given classification error rate, the optimized wavelet model tolerates a larger noise level than the non-optimized one. In Fig. 8, for a 15% classification error rate, the non-optimized wavelet model requires the noise in the co-averaged interferograms to have a standard deviation smaller than 0.30, whereas the optimized model achieves that classification performance with a noise standard deviation of 0.35. Since increasing the integration time improves the noise standard deviation by the square-root of the number of observations, the integration time required for our optimized model would be 27% smaller than the time required for the non-optimized model of the same classification capability.

## References and Links

**1. **D.M. Mittleman, R.H. Jacobsen, and M.C. Nuss, “T-Ray Imaging,” IEEE J. Sel. Top. Quantum Electron. **2**, 679–692 (1996). [CrossRef]

**2. **X.-C. Zhang, ‘Next Rays? T. Ray!’, Plenary session, 26^{th} International Conference on Infrared and millimeter waves, Toulouse France, September 2001

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

**4. **Z. Jiang and X.-C. Zhang, “Single-Shot Spatial-Temporal THz Field Imaging,” Opt. Lett. **23**, 1114–1116 (1998). [CrossRef]

**5. **Z. Jiang and X.-C. Zhang, “2D measurement and spatio-temporal coupling of few-cycle THz pulses,” Opt. Express **5**, 243–248 (1999), http://www.opticsexpress.org/oearchive/source/13775.htm. [CrossRef] [PubMed]

**6. **T. Löffler, T. Bauer, K.J. Siebert, H.G. Roskos, A. Fitzgerald, and S. Czasch, “Terahertz dark-field imaging of biomedical tissue,” Opt. Express **9**, 616–621 (2001), http://www.opticsexpress.org/oearchive/source/37294.htm. [CrossRef] [PubMed]

**7. **S.W Smye, J.M. Chamberlain, A.J. Fitzgerald, and E. Berry “The interaction between Terahertz radiation and biological tissue,” Phys. Med. Biol. **46** No 9 R101–R112 (2001). [CrossRef] [PubMed]

**8. **A.J. Fitzgerald, E Berry, N.N. Zinovev, G.C. Walker, M.A. Smith, and J.M. Chamberlain, “An introduction to medical imaging with coherent terahertz frequency radiation,” Phys. Med. Biol. **47** No 7 R67–R84 (2002). [CrossRef] [PubMed]

**9. **P. Haring-Bolivar, M. Brucherseifer, M. Nagel, H. Kurz, A. Bosserhoff, and R. Büttner “Label-free probing of genes by time-domain terahertz sensing,” Phys. Med. Biol. **47**, 3815–3822 (2002). [CrossRef] [PubMed]

**10. **S. Hadjiloucas, L.S. Karatzas, and J.W. Bowen, “Measurements of Leaf Water Content Using Terahertz Radiation,” IEEE Trans. Microwave Theory Tech. MTT. **47**, 142–149 (1999). [CrossRef]

**11. **P.Y. Han, G.C. Cho, and X.-C. Zhang, “Time-domain transillumination of biomedical tissue with terahertz pulses,” Opt. Lett. **25**, 242–244 (2000). [CrossRef]

**12. **D.D. Arnone, C. Ciesla, and M. Pepper, “Terahertz imaging comes into view,” *in Issue April 2000 of Physics World*, (Institute of Physics and IOP Publishing Limited2000), pp. 35–40.

**13. **R.M. Woodward, B. Cole, V.P. Wallace, D.D. Arnone, R. Pye, E.H. Linfield, M. Pepper, and A.G. Davies, “Terahertz pulse imaging of in-vitro basal cell carcinoma samples,” in OSA *Trends in Optics and Photonics (TOPS)*56, *Conference on Lasers and Electro-Optics (CLEO 2001)*, Technical Digest, Postconference Edition (Optical Society of America, Washington D.C., 2001), 329–330.

**14. **P Knobloch, C. Schildknecht, T. Kleine-Ostmann, M. Koch, S. Hoffmann, M. Hofmann, E. Rehberg, M. Sperling, K. Donhuijsen, G. Hein, and K. Pierz, “Medical THz imaging: an investigation of histo-pathological samples,” Phys. Med. Biol. **47**, 3875–3884 (2002). [CrossRef] [PubMed]

**15. **R.M. Woodward, B.E. Cole, V.P Wallace, R.J. Pye, D.D. Arnone, E.H. Linfield, and M. Pepper “Terahertz pulse imaging in reflection geometry of human skin cancer and skin tissue,” Phys. Med. Biol. **47**, 3853–3864 (2002). [CrossRef] [PubMed]

**16. **S. Hadjiloucas, R.K.H. Galvao, and J.W. Bowen, “Analysis of spectroscopic measurements of leaf water content at THz frequencies using linear transforms,” J. Opt. Soc. Am A **19**, 2495–2509 (2002). [CrossRef]

**17. **R.K.H. Galvão, S. Hadjiloucas, and J.W. Bowen, “Use of the statistical properties of the wavelet transform coefficients for the optimization of integration time in Fourier transform spectrometry,” Opt. Lett. **27**, 643–645 (2002). [CrossRef]

**18. **D.M. Mittleman, G. Gupta, R. Neelamani, R.G. Baraniuk, J.V. Rudd, and M. Koch “Recent advances in terahertz imaging” Appl. Phys. B **68**, 1085–1094 (1999). [CrossRef]

**19. **B. Ferguson and D. Abbott, “Wavelet de-noising of optical terahertz imaging data,” Fluctuation and Noise Lett. , **1**, L65–L70, (2001). [CrossRef]

**20. **B. Ferguson and D. Abbott, “De-noising techniques for terahertz responses of biological samples,” Microelectron. J. , **32**, 943–953 (2001). [CrossRef]

**21. **J. W. Handley, A.J. Fitzgerald, E. Berry, and R.D. Boyle “Wavelet compression in medical terahertz pulsed imaging,” Phys. Med. Biol. **47**, 3885–3892 (2002). [CrossRef] [PubMed]

**22. **G. Strang and T. Nguyen, *Wavelets and Filter Banks*, (Wellesley-Cambridge Press, Wellesley, 1996).

**23. **I. Daubechies, *Ten lectures on wavelets*, (Society for Industrial and Applied Mathematics, SIAM, Philadelphia, 1992). [CrossRef]

**24. **P. P. Vaidyanathan, *Multirate Systems and Filter Banks*, (Prentice-Hall, Englewood Cliffs, 1993).

**25. **A.V. Oppenheim and R.W. Schafer, *Discrete-Time Signal Processing*, (Prentice-Hall, Englewood Cliffs, 1989).

**26. **B.G. Sherlock and D. M. Monro, “On the space of orthonormal wavelets,” IEEE Trans. Signal Processing , **46**, 1716–1720, 1998. [CrossRef]

**27. **J. A. Nelder and R. Mead, “Simplex method for function minimization,” Computer , **7**, 308–313 (1965).

**28. **H.C. Goicoechea and A.C. Olivieri, “Wavelength selection by net analyte signals calculated with multivariate factor-based hybrid linear analysis (HLA). A theoretical and experimental comparison with partial least-squares (PLS)”, Analyst , **124**, 725–731 (1999). [CrossRef]

**29. **E. R. Malinowski, *Factor analysis in chemistry*, (Wiley, New York, 1991).