1. Introduction

Typically, terahertz (THz) pulse is a wave consisting of several oscillations of an electrical field with a corresponding spectrum in the range 0.1-10 THz. Such THz radiation finds its wide application in systems for detecting various substances (drugs, explosives) hidden behind optically opaque obstacles, as well as in methods of spectroscopy [1, 2] of various biological tissues [3–5]. Due to the development of femtosecond laser technology, microelectronics, and digital signal processing methods, it became possible to receive and detect extremely short in the number of oscillations electromagnetic THz waves [6].

There are several approaches to THz generation techniques including optical rectification of femtosecond pulses [7], excitation of THz waves by femtosecond pulses in the near-surface layer of semiconductors [6, 8], THz generation in optical filament [9–12] and in jets of water and ethanol mixtures [13]. Methods of detection are also subdivided according to the type of receiver system: electrooptical detection [8, 14, 15], photoconductive antenna [16, 17], cross-polarized scheme [18]. All these methods can register time-domain THz pulses, in form of spatio-temporal profiles of the electromagnetic field. This allows different imaging applications.

Holographic approach occupies an important place among the various methods of imaging techniques in THz range. Note, that pulsed THz holography completely differs from continuous wave (CW) methods [19]. Phase information, however, due to incoherent registration of the THz signal in this CW layout is available only by the interferometric [20] or holographic [21–25] methods, which requires a finely adjusted and more complex setup. Moreover, monochromatic CW THz radiation do not provide the information about wavelength dependence of wavefield amplitude and phase characteristics. The exception could be the phase sensitive fiber-based THz spectrometers which are now commercially available [26]. In contrast, broadband pulsed THz radiation generated by means of ultra short laser pulses contains this spectral information and allows the coherent time-resolved amplitude/phase detection.

Method of THz pulse time-domain holography (THz PTDH) [27–29] based on classical time-domain spectroscopy (TDS) has, however, several significant peculiarities in comparison with it. Firstly, THz PTDH is performed in the collimated beam, thus the investigated object can be stationary placed and only temporal scan is needed [30]. Secondly, such optics-free configuration in collimated beam allows higher reconstruction quality and is not limited by the beam waist size as in TDS [31]. Moreover, the resolution in PTDH can be increased incorporating multiple frequencies to achieve sub-wavelength resolution [32, 33], thus realizing hyperspectral data model.

Theoretically, THz PTDH method is based on the spectral approach for description of propagation of radiation, the spectrum of which, both temporal and spatial, can be wide [34]. One of the advantages of this approach is that the equations of the spectral evolution do not require approximation of either a slowly varying envelope or a slowly changing profile. But most importantly, the spectral equations, in contrast to the field equations, allow to describe the non-paraxial diffraction of single-cycle THz pulses [35, 36]. Moreover, resultant field propagation equations in the spectral approach have a simple algebraic form. In particular using the Fourier transform [37] it is convenient to model and analyze diffraction phenomena numerically and qualitatively within the framework of the spectral approach. This is especially true for radiation with a wide spatial and temporal spectrum as in the case of a single-cycle THz pulse.

Thus, THz PTDH is a powerful tool for obtaining amplitude-phase THz fields with high resolution, which allows to investigate the spectroscopic information of the object, also accounting the optical parameters such as a complex refractive index dispersion of an object’s material and a propagation media. The method demonstrates the ability to reconstruct smooth and stepped relief objects [27] or samples that are transparent in the THz region, as well as to estimate the complex spatial-temporal evolution of THz field itself [29, 38]. On the other hand, THz PTDH has some valuable restrictions connected with low signal-to-noise ratio which becomes a serious issue in coherent measurements. There are different algorithms for hyperspectral data denoising. However, a comprehensive study on the use of these algorithms for volumetric THz data in the spatio-temporal and related spatial-spectral domain has not yet been investigated.

In this research we investigate the application for THz image denoising of Block-Matching Three Dimensional (BM3D) algorithms adapted for video-filtering (VBM3D) [39] and for complex-domain denoising (CDBM3D) of complex valued data [40–42]. A preliminary version of this work with a short abstract publication has been presented on OSA topical meeting [43]. This paper involves substantially more analytic, simulation and experimental development, specifically for time-domain wavefront representation improvement, spectral amplitude/phase imaging and object’s relief reconstruction.

The proposed denoising algorithms have two purposes: first, to suppress the noise in the time-space measurements and second, to suppress the noise in spectral-space data. The both filtering applied successively are targeted on improved reconstruction of phase/amplitude characteristics of the object (or wavefront) of interest.

In this paper we are focused on the following problems: the object’s relief reconstruction (as in [27, 28, 44, 45]) and the pulsed wavefront sensing aimed on the study of spatio-temporal broadband wave field dynamics (as in [29, 38]). Multiple simulation and physical experiments demonstrate high-efficiency of the developed algorithms especially for phase imaging. For both tasks, the development and utilization of high-efficient noise reduction algorithms is extremely relevant.

It should also be noted that the issues discussed in this paper engage aspects of the currently actively emerging direction of hyperspectral holographic imaging [46, 47] developed in visible frequency range [48]. Except the another frequency range and specific technology of detection of THz waves, there are also ideologically important distinctions of the THz holography we developed, what is necessacy to account during the designing denoising algorithms. Since it operates withextremely ultrashort pulses, the spatio-temporal evolution of waveform itself is of interest for research. As a consequence, in the spectral domain, spatially distiributed wave field has a deterministic spatio-spectral phase. Therefore developed algorithms should gentle work with this data.

2. Principles of THz PTDH

In this section we consider the aspects of proposed mathematical model of the THz PTDH and denoising algorithms, carefully formalizing the initial conditions and equations to describe the dynamics of a single-cycle THz fields propagation. It should be noted that this model provides equal opportunities for numerical simulations as well as for experimental measurements. Thus, direct [28, 29] and inverse problems [27] of wavefront propagation and object reconstruction could be formalized by described approach.

2.1. Formulation of the problem

Consider the evolution of the electric field of the linearly polarized THz radiation with arbitrary spatial configuration. The dielectric medium in which the THz pulse propagates is assumed to be homogeneous and isotropic with an arbitrary dependence on the frequency ν of the linear refractive index $n\left(\nu \right)$. A general scheme of the field propagation is presented in Fig. 1. In most of the considered tasks we suppose that input THz field, irradiated by the source is ideally collimated, so its wavefront has a Gaussian amplitude spatial distribution and plane phase with the transverse coordinates $\left(x,y\right)$. In the temporal domain input signal represents a single-cycle pulse with duration of several picoseconds. The further experimental setup is building in dependence on the problem to be solved. As it was already mentioned in the introduction, there are two general classes of the problems what we have considered previously: wavefront sensing and object properties measurement.

2.1.1. Object reconstruction

For object’s reconstruction the solving tasks depend on priory object information.

(i) If optical characteristics of the sample material in the THz range are known, its thickness profile can be derived from phase information. Broadband spectrum allows to improve the thickness reconstruction quality.

(ii) If a sample profile is measured by any other profilometric technique, material properties in THz range can be reconstructed for every point of the object.

(iii) For object with several materials featuring strongly distinguished differences in their spectrum THz PTDH allows determining the lateral distribution of all matter within such objects. This kind of objects could be found in different areas: explosives [44], drugs [49], different biological tissues [50–52] or any other spectral selective objects [45]. This is possible because the derivation of the thickness from the reconstructed wavefront at different frequencies will show divergences only on characteristic spectral lines where the material is in the object. Thus, it is possible to determine not only the position of the hidden matter, but also the matter itself.

2.1.2. Wavefront sensing

In addition to the ability to assess object’s properties we mainly have the tool for analyzing the wavefront propagation dynamics. The main interest of wavefront sensing is connected with a peculiarities of behavior of various THz beams and fields with complex spatial and temporal nature. This is especially actual for any spatially structured beams such as THz Gauss-Bessel beams [29] or vortex beams [38]. It is because such beams of pulsed broadband radiation contain specific properties due to the presence of a wide temporal spectrum in which the phase versus frequency dependency is a smooth piecewise continuous function. The fact of the determinism of spectral phase significantly distinguishes the problem of pulsed broadband wavefront’s sensing from hyperspectral holographic imaging in visible range. In the last case a broadband source of continuous radiation with low time coherence is used and do not provide any phase determinism.

2.2. THz PTDH Principle Scheme

Let us consider physical layout of the key THz PDTH elements to formulate then a generic mathematical model. For the integrity of the description we assume, that the initial THz spherical divergent wavefront emitted by the point source is ideally collimated. Then, it passes through a set of modulation components (Fig. 1). Under this term, we will understand both the object under study and a system of elements that define the wavefront of a specific structure, which dynamics and evolution is supposed to be investigated. In the task of beam shaping these components represent different objects, such as phase axicons [29, 53], spiral phase plates [54, 55] or specially composed systems, e.g. THz broadband uniformly topologically charged vortex beam generator [38] or even more complicated devices using active THz metamaterials [56].

On a certain output plane z₀ the wavefront acquires a specific structure determined by the properties of the modulator, or the object under study. Then the wavefront propagates through the media and undergo changes of spatio-temporal and spatio-spectral properties. At some specific plane z the electrical field $E\left(x,y,t,z\right)$ is registered and then whole broadband wavefront propagation dynamic could be revealed by means of methods of digital holography.

 

Fig. 1 Principle scheme of pulsed THz wavefront propagation. This scheme is based on the concept of THz PTDH described and experimentally approved in [27]

Download Full Size | PDF

In terms of spatial coordinates the registration of $E\left(x,y,t,z\right)$ can be performed by 2D raster scanning procedure [54] or by using detector’s matrix array [57]. Then the inverse problem could be solved, reconstructing the spatio-temporal and spatio-spectral information of THz field in the object plane z₀ or in any other plane z. The registered THz field could be described as 4D data: 3 spatial coordinates $\left(x,y,z\right)$ and temporal t or spectral ν coordinate. Thus, THz PTDH provides the coherent output data of THz field $E\left(x,y,t,z\right)$ as well as $G\left(x,y,\nu ,z\right)$.

The description which follows further formalizes the mathematical model of THz PTDH.

2.3. Mathematical model

2.3.1. Input THz wavefield

In the scheme at Fig. 1 the input real valued function of THz pulse in temporal domain can be represented by the equation, describing the single-cycle electrical field amplitude $E_{thz}\left(t\right)$ [58, 59]:

This equation is an approximation of THz pulse generated through optical rectification of femtosecond pulses in crystals [7] and by a photoconductive semiconductor emitter irradiated by pulsed femtosecond lasers [60]. This kind of near-perfect single-cycle THz wave, obtained in experiments, can be found in various papers [7, 61, 62]. Eq. (1) describes the shape of the THz pulse $E_{thz}\left(t\right)$ specified on the time axis. This signal is real-valued and does not contain the phase in an explicit form. However, the full phase of any real signal can be defined as the complex phase of its analytic signal [63]. The spectral representation$G_{thz}\left(\nu \right)$ of this temporal form can be obtained by the Fourier integral in the time domain with the analytic solution:

This spectral density is the complex function $G_{thz}\left(\nu \right)=\left|G_{thz}\left(\nu \right)\right|\exp \left(i{\phi }_{thz}\left(\nu \right)\right)$ that contains both amplitude $\left|G_{thz}\left(\nu \right)\right|$ and phase ${\phi }_{thz}\left(\nu \right)$ components. The spatial distribution of this initial THz beam will be considered in the axially symmetric Gaussian form with transverse spatial size ρ:

2.3.2. Object

Consider the formalization of the frequency dependent spatial modulation components in case of amplitude-phase object. Let the object transmittance function be $O\left(x,y,\nu \right)=T\left(x,y,\nu \right)\exp (i\cdot {\phi }_{obj}\left(x,y,\nu \right))$, where T is an amplitude transmittance and ${\phi }_{obj}$ is a phase delay caused by the object.

Knowing the information about the thickness $H\left(x,y\right)$, complex refractive index dispersion $n_{obj}\left(x,y,\nu \right)$ and assuming the refractive index of media around the object $n_{ref}\left(x,y,\nu \right)=1$ (e.g. for the air) we can formalize the complex amplitude-phase transmittance caused by the object:

In the present study, we suppose the approximation where the THz beam is collimated and is not affected by diffraction during small propagation distances from the THz source to the set of modulation components and between them [64].

2.3.3. Wavefront propagation

To perform propagation of the initial wavefront to the plane z, we use the angular spectrum method. Two-dimensional complex wave field $G\left(x,y,\nu ,z_o\right)$ is decomposed to the angular spectrum by the spatial two dimensional Fourier integral:

The non-paraxial propagation of electromagnetic waves in the spectral approach is theoretically described in the works [34, 36]. In these papers authors have demonstrated the applicability of the spectral approach for broadband single-cycle THz radiation. The equation for the linearly polarized transverse spectral component’s propagation in homogeneous isotropic dielectric media with an arbitrary complex refractive index has been derived (see eq. (30) in the work [34]).

In numerical calculations and in THz studies it is convenient to use the linear frequency rather than angular:

Thus, the equation for transverse spectral component’s propagation can be formalized as:

Further, passing back from the spatial frequencies ($f_x,f_y$) to the original coordinates $\left(x,y\right)$ by the inverse 2D Fourier integral, we obtain the spectrum of the field propagated over the distance z:

Then, the 1D inverse Fourier transform for the propagated spectrum $G\left(x,y,\nu ,z\right)$ provides the field $E\left(x,y,t,z\right)$, which corresponds to the pulse temporal form:

We assume that the observation of E is corrupted by the additive zero-mean Gaussian noise and can be presented as follows:

This noisy data represents the model of measurements which can be experimentally obtained. In that case, the registered spatial time-domain measurements at the plane z = z_r can be treated as 3D arrays $\tilde{E}\left(x,y,t,z\right)$ with the coordinates $\left(x,y,t\right)$.

A reconstruction of the spatio-spectral wavefronts $G\left(x,y,\nu ,z\right)$ from the measurements $\tilde{E}\left(x,y,t,z\right)$ is the inverse problem usually ill-posed. The latter fact makes this reconstruction difficult as leads to amplification of measurement noise in the variable of interest $G\left(x,y,\nu ,z\right)$ which is the most informative for characterization of the object. Nevertheless, $G\left(x,y,\nu ,z\right)$ is the most informative variable.

A reconstruction of this wavefront allows also to study the dynamics of wavefront propagation as well as investigate the spectroscopic properties of the object [27], accounting a complex refractive index dispersion of the object’s material [28].

2.3.4. Hyperspectral data denoising

One of the crucial points of the considered problem is that the object to be reconstructed $O\left(x,y,\nu \right)$ is complex valued and each of the 2D object images obtained from O for the fixed z, in general different for each frequency, is complex valued having 2D phase- and amplitude-images. Another important point is that these images are obtained from indirect observations as solutions of the ill-posed inverse problem. The later conventionally leads to a serious amplification of noisy components in the object reconstructions, in particular, as speckle noise in phase images.

In this paper the denoising is implemented as a two stage procedure applied to two different groups of variables and, respectively, two different types of the filters are applied. These 2D variables are stacked together and represented as 3D cubes with spatial coordinates and the third coordinate for the time t in the first group of variables and for the frequencies ν in the second group of variables. Thus, we deal with space-time and space-frequency 3D cubes of variables for filtering.

We define the 3D space-time cube as a set of real-valued noisy observations (eq. 11):

Further, we define the 3D space-frequency cube as a set of complex-valued variables in the time Fourier transform domain as

The coordinate z may take values from z₀ till z_r. Then, this cube can be used as data for reconstruction of the wavefront for all z in this interval. For the object plane this $z=z_0$.

Thus, in what follows, we deal with two types of 3D cubes of variables: real-valued space-time $Q_T\left(x,y,t,z_r\right)$ and complex-valued space-frequency $Q_{\mathcal{V}}\left(x,y,\nu ,z\right)$.

One of the most powerful and effective modern filtering approaches is based on the grouping and collaborative filtering paradigm embodied in the BM3D image denoising algorithm [65]. This algorithm is based on sparse transform domain representations of filtered variables. The enhancement of the sparsity is achieved by grouping similar 2-D fragments (patches) of the image into 3-D data arrays which are called “groups”. Due to the similarity between the groupedpatches, the noise can be well separated in the transform domain. In this way, the collaborative filtering of the patches reveals even the finest details of the images destroyed in the original noisy data.

As it is already mentioned, in this work we use the extension of the BM3D algorithm to video-filtering VBM3D approach [39]. The 3D grouping in this algorithm is performed by a specially developed predictive block-matching spatio-temporal search similar to the technique used for motion estimation. This grouping significantly reduces the computational cost of the search for similar blocks and makes no difference between the patches in same and different video-frames. A two-step video-denoising algorithm is developed, where the predictive-search block-matching is combined with collaborative hard-thresholding in the first step and with collaborative Wiener filtering in the second step. At a reasonable computational cost, this algorithm achieves high efficiency denoising results in terms of both PSNR and visual quality. Contrary to BM3D designed for denoising 2D data the VBM3D algorithm is directly applicable for filtering real-valued space-time 3D cubes.

Applying VBM3D to the data $Q_T\left(x,y,t,z_r\right)$ we can represent the filtered estimate of $\tilde{E}\left(x,y,t,z_r\right)$ in the form

Hyperspectral imaging (HSI) is originated from the earth remote sensing [66, 67] and widely explored for numerous applications, including control of water and vegetation resources [68], control of food quality [69], biomedicine [70], etc. Acquiring two-dimensional images across a wide range of electromagnetic spectrum with hundreds and thousands spectral channels HSI provide valuable information on material, size, shape and other characteristics of objects of interest.

Conventionally, 2D measured hyperspectral images are stacked together and represented as 3D cubes with the third coordinate for the frequencies. Recently, HSI becomes a hot topic in spectrally resolved digital holography for hyperspectral phase imaging (e.g. [29, 47, 71]).

The averaging by the sample mean over the frequency is used routinely for denoising [71] in this sort of the problems. Naturally, this approach is able to filter out noise but may result in serious oversmoothing the true signal provided averaging over a large wavelength interval. Application of the more sophisticated BM3D algorithms for separate filtering of phase and amplitude appeared as a much more efficient tool [72].

The further development of the BM3D concept is targeted to deal with complex valued noisy data directly. Let $G\left(x,y,\nu ,z\right)\subset {ℂ}^{N\times M}$ be a slice of the hyperspectral noisy cube (eq. 9) of the size $N\times M$ on $\left(x,y\right)$ provided a fixed frequency ν and a given z. The Complex Domain BM3D algorithms (CDBM3D) [41, 42, 73] filters the complex-valued images of the hyperspectral cube for the each frequency separately:

Note, that the MATLAB codes of the both VBM3D and CDBM3D algorithms are publicly available. (http://www.cs.tut.fi/~foi/GCF-BM3D/)

When the estimates $\widehat{G}(x,y,\nu ,z_0$) are found from eq. (15) for the object plane $z=z_0$, the corresponding estimates for the object $O\left(x,y,\nu ,z\right)$ are calculated due to eq. (4) by separating them from the THz spectrum $G_{thz}\left(x,y,\nu \right)$.

3. Numerical implementation of the algorithms

Temporal and spectral forms of the initial THz pulse defined by the eq. (1) are presented in Fig. 2.

 

Fig. 2 The left curve is the normalized electrical field amplitude E_thz/E₀ and the right curve is the amplitude spectrum of the corresponding input THz pulse.

Download Full Size | PDF

In our simulation and experimental study, we used the following parameters: THz pulse duration τ was 0.5 ps with the corresponding spectrum interval approximately from 0.1 to 1.2 THz. The number of points N in the temporal profile was 2048, the time window size T was equal to 16,6 ps. The investigated phase object has a form of the Greek latter λ with the width $H\left(x,y\right)$ of the maximum value equal to 1.5 mm. The object is a square in $\left(x,y\right)$ of the size 24×24 $m{m}^2$. The refractive index of the object’s material $n\left(\nu \right)$ is approximately equal to 1.46, which corresponds to the Teflon usually exploited in experiments because of the high degree of transparency in the THz frequency range [27, 74]. This phase object was studied in [27] and its above mentioned parameters were obtained in this work. Usage of the same parameters allows to compare the experimental results with the corresponding simulation tests.

In the numerical computations the Fourier integrals are carried out by the Fast Fourier Transform (FFT) algorithm: 1D FFT and 2D FFT for temporal and spatial frequencies transforms correspondingly. The temporal FFT in eq. (2) is calculated for the time window T eight times larger than 16.6 ps in order to obtain a higher frequency spectral resolution. It also often happens, that the information on the signal spectrum is available only for a positive branch of frequencies. A reconstruct of the temporal signal from its positive frequency spectrum without loss of a number of points in the signal and without resource-intensive interpolation procedures is possible using the analytic signal approach [63].

The structure of data transformation and application of the denoising algorithms is demonstrated in Fig. 3. We initialize the complex spectrum $G\left(x,y,\nu ,z_0\right)$ in the object plane according to the eq. (4). This spectrum is numerically propagated to the registration plane z = z_r, where the signal is measured. The additive Gaussian noise is added to this modeled signal according to eq. (11). The value of the standard deviation σ in simulation corresponds to σ beforehand estimated from experimental data using the algorithm $function\_stdEst2D.m$ from MATLAB tool-box LASIP (http://www.cs.tut.fi/~lasip/). The 2D wavelet domain approximation of an input is used in this standard deviation measurement algorithm.

As described above the two different filter types are applied for two different groups of data. The first data type represents the spatio-temporal noisy 3D data cube $\tilde{E}\left(x,y,t,z_r\right)$ depicted by the deep yellow color in Fig. 3. In physical experiments, it corresponds to the measured data and in simulation it is a wavefront propagated to the registration plane z_r plus noise.

We apply the VBM3D filter to this data in order to suppress this noise in time domain. It provides preliminary filtered data $\widehat{E}\left(x,y,t,z_r\right)$ according to eq. (14): the cube is depicted by the yellow color in Fig. 3. The corresponding spectral cubes are presented below these temporal cubes and calculated by FFT.

Since the main goal is connected with object’s property reconstruction, the main variables of interest are given by the complex spectrum at the object plane. Thus, we propagate the wavefield back to the object plane getting $G\left(x,y,\nu ,z_0\right)$: the green cube. Then we apply the second CDBM3D filter to this complex domain spectral cube providing the final filtered data $\widehat{G}\left(x,y,\nu ,z_0\right)$ in accordance with eq. (15). This final denoised spectral data is depicted by the blue color cube.

 

Fig. 3 Algorithm of hyperspectral data denoising in temporal and spectral domain.

Download Full Size | PDF

The both algorithms VBM3D and CDBM3D have the threshold (Th) parameters which defines the level of noise suppression. These thresholds are calculated according to the following generic formula:

The lower index t indicates that this threshold is calculated for VBM3D filtering of time domain data. For spectral complex domain data and CDBM3D we use the threshold calculated according to the formula:

The standard deviations σ_t and σ_s are calculated for the corresponding variables using the program function_stdEst2D.m mentioned above.

4. Experimental and simulation results

4.1. Experimental setup

In the experimental measurements (see the optical setup in Fig. 4) we obtain the complete spatio-temporally characterized pulsed THz electrical field $E\left(x,y,t,z_r\right)$. In this scheme the laser beam of a femtosecond laser (FL) is divided into two beams: the pump beam and the probe beam. A more powerful pump beam is used to generate a THz pulse. The probe beam is firstly passed the delay line and then focused onto InAs crystal to generate THz beam. The radiation from InAs generator is collimated and passes through the object. The coherence of optical pulses (pump and probe) and the THz pulse allows the coherent field’s measurement [75]. In this case, coherence is understood as the connection of the phase of aTHz pulse with a certain characteristic (in this case intensity) of an optical pulse. This connection with high accuracy is constant in time.

The fluctuations of the laser pulse repetition frequency or mechanical vibrations are the same for the THz pulse and the pump pulse and therefore do not violate this coherence. The pump pulse and probe pulse are replicas of the same pulse and therefore also retain coherence. Thus, THz and probe pulses are tied to each other in phase. Due to this binding the probe pulse interacts in the detector (photoconducting antenna - PCA) with the same portion of THz pulse.

The duration of the probe pulse is 20 femtoseconds, which is much less than the period of the THz pulse (300-1000 fs), so we can assume that the probe radiation interacts with some constant field. By introducing a time delay for the probe pulse with respect to the pump pulse (and the associated THz pulse), it is possible to select the relative arrival time of the probe and THz pulses in the detector and detect a different portions of the pulse with a time resolution corresponding to the probe pulse duration. By scanning the time delay interval using a mechanical delay line (DL), we get a waveform of a THz pulse. We recorded the 2-D wavefront by moving a raster scanning aperture (XY with scanning pinhole PH) placed at $z=z_r=5$ mm behind the object inside the collimated beam, thus obtaining full temporal profile in the registration plane $\tilde{E}\left(x,y,t,z_r\right)$ according to the eq. (11).

 

Fig. 4 Experimental setup: FL - Ti:Sa femtosecond laser system (20 fs, 790 nm, repetition rate 70 MHz, 300 mW), BS - beam splitter, C - mechanical chopper, M - mirror, DL - delay line, O - object, InAs - THz generator, XY - 2-D scan stage, PH - scanning pinhole, PCA - photoconductive antenna, L - lens, Lthz-lens for THz beam, LIA - lock-in-amplifier, ADC - digitizer, PC - computer.

Download Full Size | PDF

4.2. THz signal measurement

In our experiment we used PCA as a receiver for THz radiation. The photoconductive antenna is one of the most commonly used as a generator and as a receiver of THz radiation. PCA consists of two metal electrodes located at some distance from each other on a semiconductor substrate. When the gap between the electrodes is illuminated by an ultrashort laser pulse, the concentration of charge carriers in the semiconductor increases sharply for a short time (about units or tens of picoseconds).

The resulting free carriers are accelerated by the field applied to the gap, resulting in a short-term current pulse $J\left(t\right)$. Thus, the ultrashort laser pulse acts as an ultra-fast switch for the antenna, transferring it from the insulating to the conducting state. As the resultant photocurrent $J\left(t\right)$ changes in time, it emits an electromagnetic pulse whose shape is determined by the time derivative of the photocurrent function.

Similarly, the PCA can be used as a THz detector. The current is proportional to the electric field of the THz pulse at the time of arrival of the probe pulse (THz field changes slowly in comparison with the duration of the laser pulse):

Here N denotes the average electron concentration, e is an elementary charge, μ is electron mobility, T is the time delay between THz and probe pulse.

Scanning the delay T, we measure the waveform of THz pulse. One of the measured THz pulse in the beam center is illustrated in Fig. 5(a). We calculate also the standard deviations σ of the random noise by averaging pixel-wise observations over the time interval. This distribution of σ is random without any traces of the THz signal, Fig. 5(b). The histograms of these σ are shown in Fig. 5(c-d). We may conclude from these histograms, that the distribution of σ is well concentrated around the peak value in a quite narrow interval. Due to this point, we assume in our simulation tests that the noise is additive i.i.d with spatially invariant standard deviation σ.

 

Fig. 5 (a) is the profile $\tilde{E}\left(t\right)$ in THz beam center; (b) is a spatial distribution of estimated standard deviation σ; (c) and (d) are histograms of σ and normalized σ, respectively.

Download Full Size | PDF

4.3. Phase images denoising: simulation results

In simulation results, we compare the reconstructed phase images with the given true phase. The Root Mean Square Error (RMSE) is used for characterization of the accuracy of this reconstruction. Under ‘true spectral phase’ or ‘true phase’ we assume the phase of the wavefront $G\left(x,y,\nu ,z_o\right)$ just behind the object plane as it is defined by eq. (4). To clearly compare phase images in these results, we applied phase unwrapping: 1D unwrapping is applied for $\phi \left(\nu \right)$ in each $\left(x,y\right)$ point and then for each frequency ν we applied 2D phase unwrapping [76].

We calculate RMSE values for the differences of the true spectral phase with the phase after VBM3D but before CDBM3D filtering (data shown as the green cube in Fig. 3), and the phase of FFT of the observed noisy data back-propagated to the object plane (data shown as as the white cube with black frame in Fig. 3). The curves of these RMSEs versus the frequency ν are shown in Fig. 6 by the green and dotted black curves, respectively. The phase after CDBM3D filtering is shown as the blue curve.

A comparison of these curves allows to evaluate the efficiency of the applied filters. The maximum values of RMSEs are for the black curves corresponding to the noisy observed data. Naturally, this data are the most noisy.

The further results with smaller RMSE, green curves, correspond to the filtered data and then back propagated to the object plane.

The best results with smallest RMSEs, blue curves, correspond to the reconstructions after second CDBM3D filtering.

We introduce also a ‘reference’ phase image, which is a free-noise reconstruction of the true phase for each spectral component (blue dotted-line in Fig. 6 which corresponds to white cube with blue frame). Formally, it can be interpreted as a data processing from the initial white cube (see Fig. 3) and going to the final blue cube with parameters σ = 0, $t{h}_t=0$ and $t{h}_s=0$. Thus, it is the reconstruction of the true phase provided noiseless observations.

Compare these dotted blue curves with the solid blue curves (results after CDBM3D) we may conclude that the proposed filtering provides nearly perfect results close to those which can be achieved, provided that the observations are ideal noiseless.

The calculations for the results in Fig. 6 are produced for various values of $t{h}_t$ in eq. (16), what demonstrates importance of the proper selection of this parameter.

 

Fig. 6 RMSE frequency dependencies for the stages of data filtration in spectral domain. Black dotted curve is a reconstruction from the initial noisy data without any filtering. Blue dotted curve is a ‘reference’ free-noise image reconstruction. Green and blue curves correspond to the reconstructed spectral data after VBM3D and CDBM3D. a) corresponds to th_t = 0, b) th_t = 0.1; c) th_t = 0.2; d) th_t = 0.3; e) th_t = 0.4; f) th_t = 0.5; g) th_t = 0.8; h) th_t = 1.

Download Full Size | PDF

The main impact in filtering is associated with the complex domain algorithm CDBM3D due to filtration of amplitude and phase images together. This correlation between phase and amplitude improves the image quality even if VBM3D in time domain is not applied (case (a) in Fig. 6).

In order to demonstrate the importance of the joint processing of phase and amplitude, we show in Fig. 7 the results which are obtained when instead of CDBM3D we apply BM3D separately for phase and amplitude for individual frequencies. The RMSE decreases only for a much higher value of $t{h}_t>2$ and with a significant sawtooth like oscillations. For the frequencies upper 0.4 THz the RMSE became even higher than initial. Thus, the advantage of CDBM3D is clear.

 

Fig. 7 RMSE frequency dependencies for the stages of data filtration in spectral domain by BM3D applied separately for amplitude and phase. Green and blue curves corresponds to the spectral data in Fig. 3. Blue dotted line is a ‘reference’ free-noise image reconstruction. (a) corresponds to th_t = 0.4, (b) th_t = 1; (c) th_t = 2; (d) th_t = 3.

Download Full Size | PDF

To observe the impact of time-domain filtering we plotted 2D image of RMSE via frequency ν and time domain threshold factor $t{h}_t$ (see Fig. 8). Indeed, picture (b) is a result of applying CDBM3D filter to the data presented in (a). In terms of color cubes (a) corresponds to the green cube and (b) to the blue. From this spectral dependency of RMSE in case (b) the optimal filtration region could be identified: $t{h}_t$ in range $\left[0;0.4\right]$ provides the lowest RMSE for the frequencies 0.1-0.6 THz. The further increasing of $t{h}_t$ leads to the growth of error in the final phase images. The valuable fact is that only VBM3D filter does not allow to resolve this optimal region, but complex spectral domain filter does.

However, then this gap become more narrow via threshold factor increase. Thus, the optimal frequency-threshold region exists, concentrating approximately in range $\nu \in \left[0.1;0.6\right]$ THz and $t{h}_t\in \left[0;0.4\right]$. Note, that increasing of RMSE after 0.6 THz is associated with the phase retardation on the object greater than 2π which leads to difficulties with 2D phase unwrapping and phase image quality decreasing.

It is illustrated on cross-section graphs in Fig. 8(c,d) for fixed frequency components 0.25 THz and 0.35 THz. This is explained by the over blurring of necessary information in time domain, where all the frequency components exist simultaneously due to the broadband pulsed nature of this investigated THz signal.

 

Fig. 8 Dependency of RMSE versus frequency ν and time domain threshold factor th_t. (a) is a result before CDBM3D and (b) is after it; (c,d) depicts the cross-sections indicated by line 1,2 from case (a,b) for ν = 0.25 THz and ν = 0.35 THz correspondingly. Green and blue curves correspond to the same cube colors.

Download Full Size | PDF

Smaller values of RMSE result in denoising and visualization of phase images are illustrated in Fig. 9. These phase images are presented for the fixed frequency ν equal to 0.25 THz and 0.35 THz with $t{h}_t=0.4$.

 

Fig. 9 Spectral domain phase images $\widehat{\phi }\left(x,y,\nu \right)$ for 0.25 THz (1st row) and 0.35 THz (2nd row) for simulation with fixed threshold value th_t = 0.4. The columns correspond to the filtration stage, pointed by the colors of the spectral cubes.

Download Full Size | PDF

4.4. Phase images denoising: experimental results

Now, let us go to discussion of denoising in experimental results. Only quantitative visual effects can be evaluated since the true phase is unknown and RMSE cannot be calculated. Fig. 10 demonstrates the VBM3D filtering with $t{h}_t=0.4$ and Fig. 11 corresponds to the much higher threshold value $t{h}_t=1.5$. In the both cases we can talk about improving the quality of the reconstructed phase images and noise suppression, although over smoothing effects can be noted in Fig. 11$.$

 

Fig. 10 Spectral domain phase images $\widehat{\phi }\left(x,y,\nu \right)$ for 0.25 THz (1st row) and 0.35 THz (2nd row) for experiment with fixed threshold value th_t = 0.4. The columns correspond to the filtration stage, pointed by the colors of the spectral cubes.

Download Full Size | PDF

 

Fig. 11 Spectral domain phase images $\widehat{\phi }\left(x,y,\nu \right)$ for 0.25 THz (1st row) and 0.35 THz (2nd row) for experiment with fixed threshold value th_t = 1.5. The columns correspond to the filtration stage, pointed by the colors of the spectral cubes.

Download Full Size | PDF

The applicability of various denoising strategies depends, of course, on the task at hand: resolution enhancement of the reconstructed object’s relief or general image quality for several monochromatic frequency components. The relief analysis which follows further illustrates the denoising efficiency.

4.5. Object’s relief denoising: experimental results

We calculate the object phase as follows: ${\widehat{\phi }}_{obj}=\widehat{\phi }-{\phi }_{thz}$, where as above $\widehat{\phi }$ is the phase after the object plane, i.e. the phase of the wavefront $\widehat{G}\left(x,y,\nu ,z_0\right)$, and ${\phi }_{thz}$ is the phase of the initial THz beam. Note, that this difference of the phases for ${\widehat{\phi }}_{obj}$ is calculated for absolute (unwrapped) phases. Then, variation of the object relief $H\left(x,y\right)$ is given by the equation

‘No filtering’ relief has been already compared with the original relief measured by 3D scanner presented in paper [27]. That is why we took this ’no filtered’ data as a reference to be compared. The results for relief reconstructions are presented below (Fig. 12) without filtering and after two types of the BM3D algorithms. Fig. 12(I) corresponds to VBM3D with $t{h}_t=0.4$ and Fig. 12-(II) corresponds to the threshold value $t{h}_t=1.5$. The frequency range for summarizing in (eq.19) is chosen 0.25-0.35 THz.

In both cases image quality improvement occurs, however for higher $t{h}_t$ image becomes over smoothed: noise is deeper suppressed but the depth resolution is decreased. Thus, the balance of image denoising and resolution depth is a point to be investigated and depends on soft temporal pre-filtration by VBM3D.

 

Fig. 12 Reconstructed relief for two values of temporal VBM3D filter: I) - th_t = 0.4 and II) - th_t = 1.5. Inset (a) corresponds to the reconstruction of initial data without any denoising; (b) is the case of only VBM3D filtering; (c) is the case of both VBM3D and CD-BM3D. (d) and (e) are the cross-sections depicted by dotted lines 1 and 2 correspondingly.

Download Full Size | PDF

4.6. Time-domain wavefront denoising

The amplitude of the THz pulse in experiments may be so low in relation to the noise, low signal-to-noise (SNR) ratio, that it has a significant impact on the analysis of temporal evolution of the pulses. This problem of field time-domain analysis and visualization could also be significant in itself, especially for different spatially structured beams. For those purposes VBM3D algorithm with higher filtration parameters allows significant denoising of temporal data profiles.

Suppressing of the noise in temporal data cube by VBM3D filter is presented in Fig. 13 with comparison of simulation and experimental results. Cross-sections of time-domain filtered cube are presented for higher parameter $t{h}_t$ contrary to the values of the best RMSE in the spectral domain, because these VBM3D settings do not provide enough good denoising in temporal domain. To demonstrate the ability of the developed technique we also present the time-domain data denoising for THz vortex beam experiment partially published in [77] and for simulation of THz wavefront propagation through the random mask.

 

Fig. 13 Time domain wavefront central cross-sections for three different modulation objects: phase object in form letter “λ”, spiral phase plate with input data from [77], and phase random mask. Columns 1,2 and 3 correspond to the different stages in time domain data processing. 1 is noisy data, 2 is denoised data with VBM3D and 3 is this denoised data propagated to the object plane. This is marked by yellow colors in cubes and corresponding captions $\tilde{E}\left(x,y,t,z_r\right)$, $\widehat{E}\left(x,y,t,z_r\right)$ and $\widehat{E}\left(x,y,t,z_o\right)$.

Download Full Size | PDF

Middle row in Fig. 13 shows the example of the time variation of the signal corresponding to the THz wavefront passed through the spiral phase plate (at the central cross-section). Resultant vortex wavefront can contain dislocations or singularity points, firstly discovered by Nye and Berry [78]. Spiral phase plate exactly produces this edge relief dislocation which causes the wavefront dislocation in vortex beam. The huge advantage of coherent THz measurements is a possibility to detect directly the electrical field amplitude $E\left(t\right)$, providing the tools for a phase sensitive measurements. Thus, in such experiments the dislocations of whole broadband wave train may be observed purely in time-domain.

Compared to continuous wave monochromatic wave fields the case of ultrafast vortices is associated with the formation of a wavepacket, the space-time analysis of which with the help of THz PTDH allows to detect the dislocations in the wave field. These dislocations behave as ’markers’ in the wave field because such vortex beams are recognizable in some degree even in the presence of noise. Our denoising makes these dislocations clearly distinguishable bringing this approach to a practical application. Thus, denoising in this work is applicable not only to the study of object itself but also for analysis of wavefront evolution in time-domain.

The bottom row in Fig. 13 demonstrates a more complex case in terms of noise reduction, when THz wavefront passes through the scattering medium. Here the noise is more difficult to distinguish from the signal, which also has a stochastic spatio-temporal nature. By analogy with work [79], to evaluate the performance of our time-domain denoising technique for simulation we selected the phase object, containing three regions with various scattering features. As can be seen from the results, the filter we use selectively distinguishes the scattered field even for the higher spatial-frequency components (x in range 10-24 mm), thus allowing the study of the wavefields dynamics that have passed through the scattering media.

5. Conclusion

THz PTDH is a very promising technique for measurement and representation of THz pulse parameters in space, time, and spectral domain. It is suitable for investigation of diffractive evolution of broadband wave fields as well as for object’s properties reconstruction. This technique could be interpreted as THz analogue of the state-of-the-art spatio-temporal metrology techniques in the visible wavelength range [80, 81]. The variety of the hyperspectral 3D techniques provides a possibility to analyze 4D and various cross-sections of complex-valued arrays of measured data.

However, the low SNR typical for experiments with small scanning pinhole imposes strong restrictions on imaging quality. In this work we investigated what kind of improvement can be achieved by the modern denoising algorithm applied to the hyperspectral THz pulse time domain holography. Using the block-matching algorithms adapted for spatio-temporal (VBM3D) and spatio-spectral (CDBM3D) data, we studied performance and optimized the parameters of these algorithms.

It is proved by simulation tests that the development techniques can significantly improve both the visual quality of phase imaging as well as the numerical accuracy of the phase reconstruction. Processing of the experimental data confirms this conclusion.

We provided also the analysis of the temporal wavefront evolution. This time-domain field visualization is significant in itself, especially for different spatially structured beams (e.g. vortex, Gauss-Bessel beams, randomly structured speckle fields). For those purposes, we demonstrated the applicability of the VBM3D algorithm with higher values of the threshold parameters for time-domain data denoising. VBM3D works rather effectively for temporal data which makes it valuable for wave field evolution observation.

However, application of VBM3D to temporal arrays provides denoising in a spectral images: starting from a certain threshold $t{h}_t$ the algorithm decrease RMSE in phase images. More efficient CDBM3D algorithm applied directly for a hyperspectral data with amplitude and phase correlation can provide even lower RMSE. Combination of these two filters leads to more precise phase images quality: the soft temporal pre-filtering by VBM3D can increase further results of CDBM3D applicability in terms of object’s relief reconstruction.