Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy

Yang Chen; Shengyang Huang; Emma Pickwell-MacPherson

doi:10.1364/OE.18.001177

1. Introduction

Terahertz (10¹² Hz, THz) technology is continuing to advance and an increasing number of applications are being investigated in areas ranging from terahertz characterization of breast cancer [1] to terahertz spectroscopy of explosives [2]. Different applications impose different requirements on the technology for instance in terms of the characteristics of terahertz devices, data acquisition rates and geometry restrictions for sample access. TeraView Ltd, Cambridge UK, has developed a prototype hand-held terahertz probe which is designed to be suitable for in vivo imaging of breast cancer (in reflection geometry) and we have been working with them to improve its capabilities. Addressing the requirements for fast data acquisition and freedom of movement of the probe (the user can wand the probe over the patient) has compromised the quality of the terahertz signal reaching the sample: both the bandwidth and signal to noise have been reduced. Thus, when we performed our usual deconvolution (inverse filtering coupled with a band pass double Gaussian filter) on data from the palm of the hand acquired using the terahertz probe, we were no longer able to resolve time domain features in the impulse response function that we had seen in data from our flat-bed imaging system (the TPI Imaga 1000, TeraView Ltd) processed in this way [3,4]. Adjusting the band bass filter to remove sufficient noise resulted in loss of too much signal. Therefore we have been investigating alternative processing methods with a view to being able to improve the resolution of the impulse response function extracted from the probe data.

The time–frequency localization and regularity of wavelet basis functions mean that wavelet transforms are able to efficiently represent terahertz pulses [4–7]. Indeed Mittleman et al observed back in 1996 that terahertz pulses are similar in form to wavelet basis functions and proposed that wavelet filtering would be very beneficial to analysis of terahertz pulses [5]. Since then good results have been demonstrated by applying wavelet techniques to terahertz data. For example, Ferguson and Abbott performed wavelet denoising of terahertz measurements of biological samples [6] using a soft thresholding method [7] to successfully remove Gaussian white noise. They then coupled the wavelet denoising with Wiener deconvolution to further improve the image quality. Wiener filtering is a classical method and is a popular choice for the inverse filtering step. It requires the estimation of the noise variance at each wavelet level, and the estimation of wavelet coefficients from a “pilot” estimate of the unknown signal [8].

The combination of wavelet denoising and frequency domain inverse filtering exploits both the wavelet domain’s efficient representation of the signal and the frequency domain’s economical representation of the noise inherent in deconvolution. Such an approach is often referred to as a “hybrid” process. For example Wiener filtering [9] has been employed in a hybrid Fourier-Wavelet Regularized Deconvolution (ForWaRD) algorithm [10] applied to terahertz wide aperture reflection tomography to successfully retrieve an accurate cross-sectional image of a complex metal object [11]. As with most previous approaches, the ForWaRD algorithm computes the discrete wavelet transforms (DWT) during the wavelet shrinkage step. The shrinkage operation is arguably the most challenging problem in wavelet denoising, for which many approaches have been explored [7,12]

In time series analysis, DWT often suffer from a lack of translation invariance. This problem can be solved by using the un-decimated stationary wavelet transform (SWT). The SWT is similar to the DWT in that the high-pass and low-pass filters are applied to the input signal at each level, but the output signal is never decimated. Instead, the filters are upsampled at each level [12]. The translation-invariant property of the stationary wavelet will ensure that at each level in the SWT the lengths and relative positions of the decomposed approximation and detail coefficients are kept the same as in the time domain, such that for each level the threshold can be estimated from the preset noise intervals along the temporal axis. Thus we employ SWT in our algorithm to further improve on other existing approaches.

In this paper we propose a hybrid Frequency-Wavelet Domain Deconvolution (FWDD) for terahertz reflection imaging and spectroscopy. The proposed FWDD approach works by performing stationary wavelet shrinkage on the Wiener deconvolution of the terahertz pulses. The translation-invariant property of SWT and the localization of the temporal pulses also enable automatic estimation of the wavelet shrinkage threshold to be achieved. In our previous work we have devised a new method to calculate the baseline offset due to reflections from the lower surface of the quartz window on which the sample is measured [13]. We use this baseline method in conjunction with our aforementioned FWDD algorithm to achieve further significant improvement in both the frequency domain and time domain analysis. This combined approach simultaneously removes unwanted noise and increases the bandwidth of sample characterization such that reflection measurements of layered structures are much clearer. Experimental data validate the better performance of the proposed method compared to double Gaussian coupled inverse filtering for extracting the time domain terahertz impulse response function. We also illustrate through simulations that the resolution of the time domain response resulting from our FWDD algorithm is higher than from the ForWaRD algorithm for cases when the signal to noise ratio is low. This is particularly relevant as our terahertz probe has increased noise. The Fourier transform of the FWDD impulse response function was then used to determine the frequency domain terahertz properties. This resulted in closer agreement with transmission spectroscopy data than if the Fourier transforms of the sample and reference raw data were used (as in previous analysis methods [14]) and thus further demonstrated the success of our new approach.

2. Theory

2.1 Baseline calculation

The full details of our baseline calculation approach are given in reference [13]. For brevity we give the key points here. For each sample measurement, two reference measurements are made, one with nothing on the sample window i.e. the quartz/air interface and one with water on the sample window i.e. the quartz/water interface. The complex refractive index data of water measured in transmission is then used in Fresnel’s equations to deduce the expected reflection off water given the measured quartz/air reflection. This is then compared with the measured quartz/water reflection to determine the baseline offset due to the lower reflection. This baseline calculation method has replaced our previous measured baseline whereby an additional identical piece of quartz was placed on the quartz window such that the quartz/quartz interface was measured. The new method overcomes the experimental difficulty of obtaining perfect contact between two windows necessary for the measured baseline method to be accurate.

2.2 Double Gaussian filtered Inverse Filtering (DGIF)

Here, we outline two deconvolution approaches: direct inverse filtering and inverse filtering coupled with a double Gaussian filter, which are referred to as IF deconvolution and DGIF deconvolution for brevity. The IF deconvolution approach is based on the equation:

(g - r) = f * (h - r)

where f is the target impulse response function, and g, h and r are the measured sample pulse, reference pulse and baseline pulse, respectively. Using the IF deconvolution, we obtain f by:

FFT (f) = \frac{FFT (g - r)}{FFT (h - r)} \Rightarrow f = {FFT}^{-1} (\frac{FFT (g - r)}{FFT (h - r)})

To suppress the amplified noise effects caused by the division operation in Eq. (2), a band pass filter can be coupled into the IF deconvolution. Hitherto a double Gaussian filter has been our filter of choice:

f i l t e r = \frac{1}{HF} e^{- \frac{t^{2}}{{HF}^{2}}} - \frac{1}{LF} e^{- \frac{t^{2}}{{LF}^{2}}}

f = {FFT}^{-1} (FFT (f i l t e r) \frac{FFT (g - r)}{FFT (h - r)})

where t represents the time axis with zero in the middle. The parameters HF (relating to the high frequency cut off) and LF (relating to the low frequency cut off) in the filter function need to be set manually. Though suppressing noise to some extent, this DGIF deconvolution tends to over-smooth the resulting deconvolved time-domain pulses.

2.3 Hybrid Fourier-Wavelet Domain Deconvolution (FWDD)

The proposed Frequency-Wavelet Domain Deconvolution (FWDD) approach performs deconvolutions of terahertz pulses by first employing frequency domain Wiener filtering and then attenuating the leaked noise through stationary wavelet shrinkage. We use SWT with Daubechies wavelets because they have similar shape and length to our typical impulse function and produce good results in our experiment. If a function is s-times continuously differentiable, then it is said to have regularity s. The regularity or smoothness of our FWDD result is dependent on the regularity of the wavelet transform. Regularity of the wavelet transform is determined by the length of the filter used in the wavelet transform; the higher the regularity, the longer the filter length. Generally, the higher the differentiable order of the target signal, the higher the regularity of the wavelet representing the signal needs to be. In this work we choose the Daubechies4 (db4) wavelet as it has an appropriate order of differentiability to represent our terahertz signal.

It is observed that our detected terahertz pulses have a non-zero signal before the reflection off the quartz/sample interface. Since the signal before the sample reflection should be zero we deduce that this must be due to system noise or fluctuations. Similarly, when measuring a semi-infinite sample layer, we have observed that there is also a non-zero signal after the reflection due to the sample, but if there were no noise then this should be zero too. We therefore treat this region as noise. Thus, by analyzing the two intervals before and after the main sample reflection, the magnitude and distribution properties of the system noise can be obtained. This is the basis of the threshold estimation for the wavelet shrinkage in our proposed FWDD approach. The new approach consists of the following steps:

2.3.1 Wiener filtering

Frequency domain Wiener filtering is performed through following Eqs. (5-7) to achieve a crude deconvolution but ensuring that there is no blurring of the main impulse function:

f_{W i e n e r} = {FFT}^{- 1} (\frac{FFT (g - r)}{FFT (h - r)} (\frac{{| FFT (h - r) |}^{2}}{{| FFT (h - r) |}^{2} + β \frac{N σ^{2}}{S}}))

S = \frac{{‖ (g - r) - mean (g - r) ‖}_{2}^{2} - N σ^{2}}{{‖ (h - r) ‖}_{2}^{2}}

mean (g - r) = \sum_{t} \frac{(g (t) - r (t))}{N}

where S is the mean power spectral density of subtracted measurements

(g - r)

. N is the total number of points in t time domain axis. The L_p norm operator is denoted by

‖ \cdot ‖

and is defined:

{‖ X ‖}_{p} : = {(\sum_{i = 1}^{n} {| x_{i} |}^{p})}^{\frac{1}{p}}

.

N σ^{2}

is the estimated noise power and the noise variance

σ^{2}

is estimated from the subtracted measurements

(g - r)

using the median estimator on the finest level wavelet coefficients of

(g - r)

[7]. The regularization parameter β can be used to modulate the effect of Wiener filtering in the tradeoff between noise-suppressing and pulse-preserving. Manually setting β to be a small value between 0.001 and 0.05 avoids over smoothing and any noise remaining after the Wiener filtering can be well addressed by applying the following stationary wavelet shrinkage step.

2.3.2 Stationary wavelet shrinkage

The translation-invariant SWT is used to preserve the relative temporal positions in wavelet decomposition and reconstruction. Along the temporal axis ℓ at each level k, it transforms a 1-D signal $f_{F} (n)$ into the approximation coefficients vector $c A_{k, l}$ and detail coefficients vector $c D_{k, l}$ by convolving $f_{F} (n)$ with a low-pass filterΨand a high-pass filterΦ respectively. In contrast to the well-known general discrete wavelet transform, SWT does not down sample the signal and instead upsamples the filters at each decomposition level. Figure 1 shows the schematic decomposition of SWT:

Fig. 1 The schematic decomposition of SWT

Download Full Size | PDF

Wavelet shrinkage is a signal denoising technique based on the idea of thresholding the wavelet coefficients. Wavelet coefficients having small absolute value are considered to encode mostly noise and very fine details of the signal. In contrast, the important information is encoded by the coefficients having large absolute values. Removing the small absolute value coefficients and then reconstructing the signal should produce a signal with less noise. After applying SWT on the result from Wiener filtering, we obtain the approximation and detail coefficients $c A_{k, ℓ}$ and $c D_{k, ℓ}$ , at each level k [15]. We used a maximum level of 5 for the wavelet decomposition as no significant improvement was observed for higher levels to justify the extra computational expense.

[c A, c D] = SWT(f_{W i e n e r})

For the wavelet shrinkage at each level k, we need the threshold parameters $T_c D_{k}$ for the detail coefficients. We define $T_c D_{k}$ to be the maximum magnitudes of the detail coefficients within the indicated noise intervals at each level k. We note that at each level k in the SWT, the lengths and relative positions of the decomposed approximation coefficients $c A_{k, ℓ}$ and detail coefficients $c D_{k, ℓ}$ are kept the same as in the temporal axis of the original data for $f_{W i e n e r}$ . The positions of the noise intervals (indicated by the red rectangles in Fig. 2 ) are thus also kept unchanged for each level k. As a default, we choose the length of each red box to be one quarter of the length of the measured signal and position them to avoid the main pulse. For example, for a measurement with 512 data points, the box preceding the main pulse starts 10 points from the beginning of the signal and finishes at point 10 + 512/4 = 138, and the box following the main pulse starts at 512-138 = 374. Their separation in this case is therefore 374-138 = 236 data points. The separation of the boxes is not too crucial, but we sometimes need to change the position if they overlap with the main pulse, this occasionally happens for measurements with fewer data points (eg 256). By considering the noise after the main pulse as well as before it, we are able to improve the resulting impulse function than if only the noise interval before the pulse was used. This is because in many cases the magnitude of the noise after the main pulse is greater than that preceding the main pulse. With these acquired $T_c D_{k}$ at each level k we then perform shrinkage to obtain the estimated $c {\hat{D}}_{k, ℓ}$ by thresholding the coefficients $c D_{k, ℓ}$ along the whole temporal axis:

Fig. 2 A flow chart to highlight the extra processes involved in our proposed approach compared to our original method.

Download Full Size | PDF

c {\hat{D}}_{k, ℓ} = s i g n (c D_{k, ℓ}) {(| c D_{k, ℓ} | - T_c D_{k})}_{+}

Here, soft thresholding is used with $s i g n (x) = x / | x |$ . ${(x)}_{+} = 0$ if $x < 0$ , and ${(x)}_{+} = x$ if $x \geq 0$ . The maximum detail coefficients’ magnitudes from the main pulse are generally higher than those of the noise within the indicated noise intervals. Coefficients with larger magnitudes tend to represent useful information and smaller coefficients tend to represent noise, therefore in the shrinkage process, we only keep the coefficients which are higher than the highest coefficient from within the noise. Then we can obtain the deconvolved pulse $f_{F W D D}$ by performing an inverse stationary wavelet transform (ISWT, the inverse process of SWT) from the approximation coefficients $c A_{k, ℓ}$ and the thresholded detail coefficients $c {\hat{D}}_{k, ℓ}$ at each level k.

f_{F W D D} = ISWT ([c A, c \hat{D}])

2.4 Calculation of the complex refractive index

We use Fresnel theory to determine the complex refractive index of the sample from sample and reference measurements. The resulting equation (as explained in [9]) is:

{\tilde{k}}_{s a m p l e} \cos θ_{s a m} = \frac{(1 - M) {\tilde{k}}_{q u a r t z}^{2} θ_{q u a r t z}^{2} + (1 + M) {\tilde{k}}_{q u a r t z} θ_{q u a r t z} {\tilde{k}}_{a i r} \cos θ_{a i r}}{(1 - M) {\tilde{k}}_{a i r} \cos θ_{a i r} + (1 + M) {\tilde{k}}_{q u a r t z} θ_{q u a r t z}}

where,

{\tilde{k}}_{s a m p l e} \sin θ_{s a m} = \frac{1}{2}

and

{\tilde{k}}_{s a m p l e} = n - j \frac{c}{2 ω} α

with n and α representing the refractive index and absorption coefficient respectively. In our original approach, M was the frequency ratio of the sample and reference, given as:

M = \frac{FFT (g - r)}{FFT (h - r)}

In this new FWDD approach, we use the Fourier transform of the impulse function resulting from the inverse stationary wavelet transform to determine M such that:

M = FFT (f_{F W D D})

The steps of the FWDD process are outlined and compared to our previous methods in Fig. 2.

3. Experimental methods and equipment

We initially apply our theory to the same data previously used to illustrate the improvement due to our baseline calculation method [13]. In this way we can highlight the additional improvement due to our FWDD approach. These data were acquired from the reflection geometry flat-bed system, TPI Imaga 1000, TeraView Ltd. Details about the system and experimental procedure are given in reference [3] and [13] respectively.

To further test our theory we apply it to in vivo measurements of the palm of the hand acquired using the terahertz probe by TeraView Ltd. In contrast to the TPI Imaga 1000, which employs free space optics, the probe system utilizes optical fibers to guide the laser beam to the terahertz devices. Short 90 fs pulses (centered at 800 nm) from a Ti:Sapphire laser are guided along an optical fiber shielded with an electrical cable to the probe head where a photoconductive emitter and receiver are embedded for generating and detecting THz pulses. The bandwidth achieved by this system is 0.1-2 THz. More details are given in reference [16].

4. Results and discussion

4.1 Frequency domain validation

Reflection data of isopropanol from TPI Imaga 1000 were processed using our new approach with parameter β in the Wiener filtering step set to be 0.0015 and the resulting refractive index and absorption coefficient are compared to our results from previous methods in Fig. 3 .

Fig. 3 The refractive index (a) and absorption coefficient (b) of isopropanol calculated using our newly proposed FWDD approach and compared to our previous methods and transmission spectroscopy data.

Download Full Size | PDF

From the traces in Fig. 3 we can clearly see the improvements as our processing method becomes more refined. The most significant improvements are seen in the absorption coefficient. This is most likely to be because the calculation of the absorption coefficient in reflection geometry is predominantly dependent on the phase information, and this is very sensitive to the experimental set up. When the measured baseline was used, the absorption coefficient barely matched the transmission spectroscopy data up to 0.5 THz. The calculated baseline method extended this up to about 1 THz; and using the FWDD approach combined with the calculated baseline has extended this further to about 1.8 THz. Thus the combined approach has nearly quadrupled the useable bandwidth compared to the standard measured baseline and processing approach.

The reason for this improvement is that the FWDD deconvolution produces a more correct impulse function of isopropanol in the time domain. This is then used in the calculation of the complex refractive index and thus results in a very close match with the standard transmission data up to about 1.8 THz. To illustrate the improvements in the time domain impulse response function we use in vivo measurements of skin as this is a multi-layered structure and better shows the strength of our approach.

4.3 Time domain validation and in vivo skin measurements

To illustrate the effectiveness of the FWDD deconvolution in the time domain processing we took in vivo reflection data of the palm of a hand using the terahertz probe system. As depicted in Fig. 4 , the skin of the palm of the hand has a layered structure: the epidermis lies beneath the stratum corneum and has a different refractive index. There are two troughs in the reflected waveform. This is because the incident terahertz light is reflected when it enters a medium of different refractive index. The first reflection is due to the interface between the quartz window and the stratum corneum, and the second reflection is due to the interface between the stratum corneum and the epidermis. The reflections have the opposite phase to the incident pulse because it has entered a medium of higher refractive index. The separation (optical delay) between the troughs is dependent on the refractive index and the thickness of the stratum corneum. Since the skin has a high water content which attenuates the signal, we are not able to see reflections off the epidermis/dermis interface. In some cases there is a narrow peak before the first trough; this is attributed to surface dryness of the stratum corneum [3].

Fig. 4 (a) A photograph to illustrate the measurement of the palm using the terahertz probe; (b) A schematic diagram of the skin layers contributing to the reflections observed in our terahertz measurements.

Download Full Size | PDF

One of our previous studies measured the palms of several volunteers using the TPI Imaga 1000 [3]. We tried to perform a similar study using the terahertz probe, however due to the reduced signal to noise of the probe (a side effect of its improved acquisition times and fibre coupled devices), our previous processing methods using DGIF did not produce satisfactory results even when the calculated baseline was used. Figure 5(a) is a palm measurement taken using the flat-bed system – it is very similar to the palm response function illustrated by Cole et al [4]. The resulting response functions from the reflection measurements of the first author’s palm measured by the probe are plotted in Fig. 5(b)-5(d). In Fig. 5(b), the data are processed using the same double Gaussian filter that was used to obtain the response function in Fig. 5(a) however the two reflections from the skin are not resolved. The double Gaussian filter parameters are adjusted in Fig. 5(c) to try to resolve the two reflections, but increasing the resolution also increased the noise. However as illustrated in Fig. 5(d), the response function from our FWDD approach (with β now set to 0.01) is both able to resolve the reflections and is smooth – the useful high frequency components have been preserved whilst the noise and fluctuations have been removed. Using this approach we are able to obtain similar response functions to those in previous studies using the flatbed system by Cole et al [4] and Pickwell et al [3].

Fig. 5 The resulting response functions calculated using palm data from a) the flat-bed imaging system, processed using a double Gaussian filter; b) the probe, processed using the same double Gaussian as in (a); c) the probe, processed with various double Gaussian filters; and d) the probe, processed using the new (FWDD) approach.

Download Full Size | PDF

The probe was moved across the palm to take measurements at different positions on the palm [see the inset of Fig. 6(a) ]. By processing these data using our FWDD approach, we see in Fig. 6(a) that the optical delay between the troughs increases for measurements closer to the center of the hand. Assuming that there is minimal change in the refractive index of the skin, this means that the stratum corneum becomes thicker as we move towards the centre of the palm. By plotting the normalized amplitude of the waveform (E(n)/E_max) against the optical delay for each of the positions measured, we have formed a b-scan image in Fig. 6(b). Plotting E(n)/E_max enhances the image by making the reflections from the surface of the palm (at E(n) = E_max) a uniform color (white). This nicely illustrates the depth variation with position and is consistent with previous results from the edge of the hand by Cole et al. (Fig. 7(b) , reference [4]).

Fig. 6 (a) The impulse functions along the palm using the FWDD approach with the inset photograph indicating the position of the measurements on the palm. (b) A b-scan image showing the variation in stratum corneum thickness with palm position. The color bar represents E(n)/E_max. The white dashed line indicates the interface between the quartz window and the stratum corneum; and the black dashed line indicates the interface between the stratum corneum and the epidermis.

Download Full Size | PDF

Fig. 7 Pulse separation when SNR is 32 dB by (a) DGIF, LF = 2048 and HF = 20; (b) FWDD, β = 0.001; and (c) ForWaRD approaches.

Download Full Size | PDF

To highlight the improvement in resolution that our FWDD algorithm has over DGIF and ForWaRD approaches, we simulate the reflections at normal incidence due to a material of refractive index 1.5 with thicknesses ranging from 200 µm down to 20 µm at 20 µm intervals and process the data using the three approaches. In Fig. 7 we perform the simulations for a system with a high signal to noise ratio (SNR). The same DGIF parameters used in the flat-bed system analysis are used for the higher SNR simulation (LF = 2048 and HF = 20). Gaussian white noise is added to the incident terahertz pulse resulting in a SNR of 32 dB. We see that both the FWDD and ForWaRD approaches are able to resolve down to 40 µm (0.4 ps of optical delay) with similar sharpness, but that the DGIF can only resolve down to 80 µm (0.8 ps of optical delay). When the SNR of the system is reduced to 22 dB the DGIF parameters are adjusted to LF = 2048 and HF = 30 so that the resulting reflections are still smooth. As illustrated in Fig. 8 , the DGIF can now only resolve down to 120 µm (1.2 ps of optical delay) and the FWDD and ForWaRD down to 80 µm (0.8 ps). However the FWDD approach results in clearer separation of the reflections than ForWaRD – this is most noticeable for 100 µm and 80 µm. Since our terahertz probe has a lower SNR than the flat-bed system this sharper resolution of the FWDD approach is very important. One reason for this is that ForWaRD’s automatic parameter settings limit its performance in high noise conditions whereas by modulating β the FWDD approach is still able to perform well. Thus our FWDD algorithm copes well with the probe data as β can be adjusted to deal with the noisier signal.

Fig. 8 Pulse separation when SNR is 22 dB by (a) DGIF, LF = 2048 and HF = 30; (b) FWDD, β = 0.005; and (c) ForWaRD approaches.

Download Full Size | PDF

The computation costs for the IFDG, ForWaRD and proposed FWDD method for the 1-D signals calculated in terms of CPU times are 0.0018 seconds, 0.87 seconds, and 0.80 respectively. The simulation was performed in Matlab 7.1 using an Intel Core^TM 2 Quad CPU, 4096 Mb RAM. It is noted that although the FWDD method uses SWT and not DWT, it is slightly computationally less expensive than the ForWaRD algorithm; this is likely to be due to the automatic parameter finding in ForWaRD. The upsampling involved in the SWT process is far more (450 times more in this case) computationally expensive than the DGIF as it creates a much larger sample set and this makes real-time processing less feasible, but still tolerable for 1-D cases.

5. Conclusions

We have shown that by using the FWDD approach combined with our baseline calculation method we are able to nearly quadruple the bandwidth of reflection spectroscopy measurements as compared to standard methods. The FWDD approach also results in vast improvements in the time domain analysis as compared to DGIF such that impulse functions of samples can be retrieved with much less noise and without loss of useful high frequency information. This enables the reflections off layered structures to be much more clearly resolved. By using the FWDD approach on reflection data acquired by the terahertz probe we are able to improve the impulse response function so much so that we are able to compensate for the reduced signal to noise performance of the probe due to its fast acquisition times and fiber coupled devices.

In conclusion, we have devised a robust data processing algorithm that is able to obtain sharper and clearer deconvolved impulse functions and more accurate complex refractive indices.

Acknowledgments

The authors would like to thank Dr David Seery, Department of Applied Mathematics  and Theoretical Physics, Cambridge University, UK, for useful discussions. We gratefully acknowledge partial financial support for this work from the Research Grants Council of the Hong Kong Government and the Shun Hing Institute of Advanced Engineering, Hong Kong.

References and links

1. P. C. Ashworth, E. Pickwell-MacPherson, E. Provenzano, S. E. Pinder, A. D. Purushotham, M. Pepper, and V. P. Wallace, “Terahertz pulsed spectroscopy of freshly excised human breast cancer,” Opt. Express 17(15), 12444–12454 (2009). [CrossRef] [PubMed]

2. M. C. Kemp, P. F. Taday, B. E. Cole, J. A. Cluff, A. J. Fitzgerald, and W. R. Tribe, “Security applications of terahertz technology,” Terahertz for Military and Security Applications 5070, 44–52 (2003).

3. E. Pickwell, B. E. Cole, A. J. Fitzgerald, M. Pepper, and V. P. Wallace, “In vivo study of human skin using pulsed terahertz radiation,” Phys. Med. Biol. 49(9), 1595–1607 (2004). [CrossRef] [PubMed]

4. B. E. Cole, R. Woodward, D. Crawley, V. P. Wallace, D. D. Arnone, and M. Pepper, “Terahertz imaging and spectroscopy of human skin, in-vivo,” Proc. Soc. Photo Opt. Instrum. Eng. 4276, 1–10 (2001).

5. D. M. Mittleman, R. H. Jacobsen, and M. C. Nuss, “T-ray imaging,” IEEE Sel. Top. Quantum Electron. 2(3), 679–692 (1996). [CrossRef]

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

7. D. L. Donoho, “De-noising by soft-thresholding,” IEE Trans. Inf. Theory 41(3), 613–627 (1995). [CrossRef]

8. K. R. Castleman, “Digital Image Processing,” Prentice-Hall, Englewood Cliffs, NJ (1996).

9. S. P. Ghael, A. M. Sayeed, and R. G. Baraniuk, “Improved wavelet denoising via empirical Wiener filtering,” Wavelet Applications in Signal and Image Processing V 3169, 389–399 (1997).

10. R. Neelamani, H. Choi, and R. Baraniuk, “Forward: Fourier-Wavelet regularized deconvolution for Ill-conditioned systems,” IEEE Trans. Signal Process. 52(2), 418–433 (2004). [CrossRef]

11. J. Pearce, H. H. Choi, D. M. Mittleman, J. White, and D. Zimdars, “Terahertz wide aperture reflection tomography,” Opt. Lett. 30(13), 1653–1655 (2005). [CrossRef] [PubMed]

12. J. C. Pesquet, H. Krim, and H. Carfantan, “Time-invariant orthonormal wavelet representations,” IEEE Trans. Acoust., Speech, Signal Process. 44, 1964–1970 (1996).

13. S. Huang, P. C. Ashworth, K. W. C. Kan, Y. Chen, V. P. Wallace, Y. T. Zhang, and E. Pickwell-MacPherson, “Improved sample characterization in terahertz reflection imaging and spectroscopy,” Opt. Express 17(5), 3848–3854 (2009). [CrossRef] [PubMed]

14. S. Y. Huang, Y. X. J. Wang, D. K. Yeung, A. T. Ahuja, Y. T. Zhang, and E. Pickwell-Macpherson, “Tissue characterization using terahertz pulsed imaging in reflection geometry,” Phys. Med. Biol. 54(1), 149–160 (2009). [CrossRef]

15. G. Nason and B. Silverman, “The Stationary Wavelet Transform and some Statistical Applications,” Wavelets and Statistics, Springer Lecture Notes in Statistics 13, 281–300 (1995).

16. P. C. Ashworth, P. O'Kelly, A. D. Purushotham, S. E. Pinder, M. Kontos, M. Pepper, and V. P. Wallace, “An intra-operative THz probe for use during the surgical removal of breast tumors,” 2008 33rd International Conference on Infrared, Millimeter and Terahertz Waves, Vols 1 and 2, 767–769 (2008).

Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy

Abstract

1. Introduction

2. Theory

2.1 Baseline calculation

2.2 Double Gaussian filtered Inverse Filtering (DGIF)

2.3 Hybrid Fourier-Wavelet Domain Deconvolution (FWDD)

2.3.1 Wiener filtering

2.3.2 Stationary wavelet shrinkage

2.4 Calculation of the complex refractive index

3. Experimental methods and equipment

4. Results and discussion

4.1 Frequency domain validation

4.3 Time domain validation and in vivo skin measurements

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (8)

Equations (13)

Optics Express