## Abstract

We investigated data denoising in hyperspectral terahertz pulse time-domain holography. Using the block-matching algorithms adapted for spatio-temporal and spatio-spectral volumetric data we studied and optimized parameters of these algorithms to improve phase image reconstruction quality. We propose a sequential application of the two algorithms oriented on work in temporal and spectral domains. Experimental data demonstrate the improvement in the quality of the resultant time-domain images as well as phase images and object’s relief. The simulation results are proved by comparison with the experimental ones.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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*_{0} 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.

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*_{0} 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]:

*E*

_{0}is the amplitude of electrical field, and

*τ*is the pulse duration.

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|\mathrm{exp}\text{}\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)\mathrm{exp}\text{}(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:

*c*is the speed of light.

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:

*σ*is the noise standard deviation and

*ε*denotes a random variable having normalized Gaussian distribution $N\left(0,1\right)$.

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 2*D* object images obtained from *O* for the fixed *z*, in general different for each frequency, is complex valued having 2*D* 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 2*D* variables are stacked together and represented as 3*D* 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 3*D* cubes of variables for filtering.

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

*X*and

*Y*spaces for $\left(x,y\right)$ coordinates of the sensor pixels and

*T*is the time interval of observations, the registration plane coordinate

*z*is fixed as

*z*=

*z*.

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

*G*is defined from eq. (10) by direct Fourier transform from temporal form $E\left(x,y,t,z\right)$.

The coordinate *z* may take values from *z*_{0} 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 3*D* 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 3*D* 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, 2*D* measured hyperspectral images are stacked together and represented as 3*D* 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 {\u2102}^{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.

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 3*D* 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.

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:

where $t{h}_{t}$ is a coefficient (factor) of the threshold parameter.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

*σ*are calculated for the corresponding variables using the program

_{s}*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).

#### 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 *σ*.

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

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.

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.

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

#### 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$.$

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.

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

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.

## Funding

Russian Foundation for Basic Research (18-32-20215$\backslash 18$) and Academy of Finland, project no. 287150, 2015-2019.

## Acknowledgments

The authors express their honor and respect to Prof. V.G. Bespalov, the founder of the THz PTDH.

## References

**1. **I. Azuri, A. Hirsch, A. M Reilly, A. Tkatchenko, S. Kendler, O. Hod, and L. Kronik, “THz Spectroscopy of 2, 4, 6-trinitrotoluene Molecular Solids from First Principles,” Bull. Am. Phys. Soc. (2018).

**2. **K. Kaltenecker, B. Zhou, K.-H. Tybussek, S. Engelbrecht, R. Lehmann, S. Walker, P. U. Jepsen, and B. Fischer, “Ultra-broadband THz spectroscopy for sensing and identification for security applications,” in 2018 43rd International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz), (IEEE, 2018), pp. 1–2.

**3. **O. Smolyanskaya, I. Schelkanova, M. Kulya, E. Odlyanitskiy, I. Goryachev, A. Tcypkin, Y. Grachev, Y. Toropova, and V. Tuchin, “Glycerol dehydration of native and diabetic animal tissues studied by THz-TDS and NMR methods,” Biomed. Opt. Express9 (2018). [CrossRef] [PubMed]

**4. **O. Smolyanskaya, E. Odlyanitskiy, K. Zaytsev, and M. Kulya, “Propagation Dynamics of the THz Radiation Through a Dehydrated Tissue by the Pulse Time Domain Holography Method,” in 2018 43rd International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz), (IEEE, 2018), pp. 1–2.

**5. **O. Smolyanskaya, Q. Cassar, M. Kulya, N. Petrov, K. Zaytsev, A. Lepeshkin, J.-P. Guillet, P. Mounaix, and V. Tuchin, “Interaction of terahertz radiation with tissue phantoms: numerical and experimental studies,” in *EPJ Web of Conferences*, vol. 195 (EDP Sciences, 2018), p. 10012. [CrossRef]

**6. **V. Bespalov, A. Gorodetskii, I. Y. Denisyuk, S. Kozlov, V. Krylov, G. Lukomskii, N. Petrov, and S. Putilin, “Methods of generating superbroadband terahertz pulses with femtosecond lasers,” J. Opt. Technol. **75**, 636–642 (2008). [CrossRef]

**7. **K.-L. Yeh, J. Hebling, M. C. Hoffmann, and K. A. Nelson, “Generation of high average power 1 kHz shaped THz pulses via optical rectification,” Opt. Commun. **281**, 3567–3570 (2008). [CrossRef]

**8. **K. Ahi, “Review of GaN-based devices for terahertz operation,” Opt. Eng. **56**, 090901 (2017). [CrossRef]

**9. **A. A. Andreev, V. G. Bespalov, A. A. Gorodetskii, S. A. Kozlov, V. N. Krylov, G. V. Lukomskii, E. V. Novoselov, N. V. Petrov, S. E. Putilin, and S. A. Stumpf, “Generation of ultrabroadband terahertz radiation under optical breakdown of air by twofemtosecond pulses of different frequencies,” Opt. Spectrosc. **107**, 538–544 (2009). [CrossRef]

**10. **V. A. Andreeva, O. G. Kosareva, N. A. Panov, D. E. Shipilo, P. M. Solyankin, M. N. Esaulkov, P. González, Alaiza de Martínez, A. P. Shkurinov, V. A. Makarov, L. Bergé, and S. L. Chin, “Ultrabroad Terahertz Spectrum Generation from an Air-Based Filament Plasma,” Phys. Rev. Lett. **116**, 063902 (2016). [CrossRef] [PubMed]

**11. **S. Smirnov, Y. V. Grachev, A. Tsypkin, M. Kulya, S. Putilin, and V. Bespalov, “Spatial-temporal dynamics of the terahertz field generated by femtosecond filament,” in *Journal of Physics: Conference Series*, vol. 735 (IOP Publishing, 2016), p. 012065.

**12. **S. Smirnov, M. Kulya, A. Tcypkin, S. Putilin, and V. Bespalov, “Detection of the polarization spatial distribution of THz radiation generated by two-color laser filamentation,” Nanosystems: physics, chemistry, mathematics8 (2017).

**13. **A. Tcypkin, E. Ponomareva, S. Putilin, S. Smirnov, S. Shtumpf, M. Melnik, S. Kozlov, and X.-C. Zhang *et al.*, “Concentration dependence of terahertz generation in jets of water and ethanol mixtures,” in *Infrared, Millimeter-Wave, and Terahertz Technologies V*, vol. 10826 (International Society for Optics and Photonics, 2018), p. 1082603.

**14. **G. Gallot and D. Grischkowsky, “Electro-optic detection of terahertz radiation,” JOSA B **16**, 1204–1212 (1999). [CrossRef]

**15. **Q. Wu, T. Hewitt, and X.-C. Zhang, “Two-dimensional electro-optic imaging of THz beams,” Appl. Phys. Lett. **69**, 1026–1028 (1996). [CrossRef]

**16. **S. Matsuura, M. Tani, and K. Sakai, “Generation of coherent terahertz radiation by photomixing in dipole photoconductive antennas,” Appl. Phys. Lett. **70**, 559–561 (1997). [CrossRef]

**17. **Z. Piao, M. Tani, and K. Sakai, “Carrier dynamics and terahertz radiation in photoconductive antennas,” Jpn. J. Appl. Phys. **39**, 96 (2000). [CrossRef]

**18. **J. Van Rudd, J. L. Johnson, and D. M. Mittleman, “Cross-polarized angular emission patterns from lens-coupled terahertz antennas,” JOSA B **18**, 1524–1533 (2001). [CrossRef]

**19. **N. Karpowicz, H. Zhong, J. Xu, K.-I. Lin, J.-S. Hwang, and X. Zhang, “Comparison between pulsed terahertz time-domain imaging and continuous wave terahertz imaging,” Semicond. Sci. Technol. **20**, S293 (2005). [CrossRef]

**20. **J. L. Johnson, T. D. Dorney, and D. M. Mittleman, “Interferometric imaging with terahertz pulses,” IEEE Journal of selected topics in quantum electronics **7**, 592–599 (2001). [CrossRef]

**21. **K. Xue, Q. Li, Y.-D. Li, and Q. Wang, “Continuous-wave terahertz in-line digital holography,” Opt. letters **37**, 3228–3230 (2012). [CrossRef]

**22. **L. Rong, T. Latychevskaia, D. Wang, X. Zhou, H. Huang, Z. Li, and Y. Wang, “Terahertz in-line digital holography of dragonfly hindwing: amplitude and phase reconstruction at enhanced resolution by extrapolation,” Opt. express **22**, 17236–17245 (2014). [CrossRef] [PubMed]

**23. **E. Hack and P. Zolliker, “Terahertz holography for imaging amplitude and phase objects,” Opt. express **22**, 16079–16086 (2014). [CrossRef] [PubMed]

**24. **M. Locatelli, M. Ravaro, S. Bartalini, L. Consolino, M. S. Vitiello, R. Cicchi, F. Pavone, and P. De Natale, “Real-time terahertz digital holography with a quantum cascade laser,” Sci. reports **5**, 13566 (2015). [CrossRef]

**25. **P. Zolliker and E. Hack, “THz holography in reflection using a high resolution microbolometer array,” Opt. express **23**, 10957–10967 (2015). [CrossRef] [PubMed]

**26. **A. Roggenbuck, H. Schmitz, A. Deninger, I. C. Mayorga, J. Hemberger, R. Güsten, and M. Grüninger, “Coherent broadband continuous-wave terahertz spectroscopy on solid-state samples,” New J. Phys. **12**, 043017 (2010). [CrossRef]

**27. **N. V. Petrov, M. S. Kulya, A. N. Tsypkin, V. G. Bespalov, and A. Gorodetsky, “Application of terahertz pulse time-domain holography for phase imaging,” IEEE Transactions on Terahertz Sci. Technol. **6**, 464–472 (2016). [CrossRef]

**28. **M. Kulya, N. Balbekin, I. Gredyuhina, M. Uspenskaya, A. Nechiporenko, and N. Petrov, “Computational terahertz imaging with dispersive objects,” J. Mod. Opt. **64**, 1283–1288 (2017). [CrossRef]

**29. **M. S. Kulya, V. A. Semenova, V. G. Bespalov, and N. V. Petrov, “On terahertz pulsed broadband Gauss-Bessel beam free-space propagation,” Sci. reports **8**, 1390 (2018). [CrossRef]

**30. **Y. Zhang, W. Zhou, X. Wang, Y. Cui, and W. Sun, “Terahertz digital holography,” Strain **44**, 380–385 (2008). [CrossRef]

**31. **J. F. Federici, B. Schulkin, F. Huang, D. Gary, R. Barat, F. Oliveira, and D. Zimdars, “THz imaging and sensing for security applications–explosives, weapons and drugs,” Semicond. Sci. Tech. **20**, S266–S280 (2005). [CrossRef]

**32. **L. Zhang, Y. Zhang, C. Zhang, Y. Zhao, and X. Liu, “Terahertz multiwavelength phase imaging without 2π ambiguity,” Opt. Lett. **31**, 3668–3670 (2006). [CrossRef] [PubMed]

**33. **L. Zhang, H. Zhong, Y. Zhang, N. Karpowicz, C. Zhang, Y. Zhao, and X.-C. Zhang, “Terahertz wave focal-plane multiwavelength phase imaging,” JOSA A **26**, 1187–1190 (2009). [CrossRef] [PubMed]

**34. **A. A. Ezerskaya, D. V. Ivanov, S. A. Kozlov, and Y. S. Kivshar, “Spectral approach in the analysis of pulsed terahertz radiation,” J. Infrared, Millimeter, Terahertz Waves **33**, 926–942 (2012). [CrossRef]

**35. **S. Kozlov and D. Ivanov, “Diffraction of a terahertz single-period train of the field at a slit,” J. Opt. Technol. **77**, 734–736 (2010). [CrossRef]

**36. **A. Ezerskaya, D. Ivanov, V. Bespalov, and S. Kozlov, “Diffraction of single-period terahertz electromagnetic waves,” J. Opt. Technol. **78**, 551–557 (2011). [CrossRef]

**37. **R. N. Bracewell and R. N. Bracewell, *The Fourier transform and its applications*, vol. 31999 (McGraw-HillNew York, 1986).

**38. **M. Kulya, V. Semenova, A. Gorodetsky, V. G. Bespalov, and N. V. Petrov, “Spatio-temporal and spatiospectral metrology of terahertz broadband uniformly topologically charged vortex beams,” Appl. Opt. **58**, A90–A100 (2019). [CrossRef] [PubMed]

**39. **K. Dabov, A. Foi, and K. Egiazarian, “Video denoising by sparse 3D transform-domain collaborative filtering,” in 2007 15th European Signal Processing Conference, (IEEE, 2007), pp. 145–149.

**40. **V. Katkovnik, I. Shevkunov, N. V. Petrov, and K. Egiazarian, “Computational super-resolution phase retrieval from multiple phase-coded diffraction patterns: simulation study and experiments,” Optica (2017).

**41. **V. Katkovnik, M. Ponomarenko, and K. Egiazarian, “Sparse approximations in complex domain based on BM3D modeling,” Signal Process. **141**, 96–108 (2017). [CrossRef]

**42. **V. Katkovnik and K. Egiazarian, “Sparse phase imaging based on complex domain nonlocal BM3D techniques,” Digit. Signal Process. **63**, 72–85 (2017). [CrossRef]

**43. **M. Kulya, N. Petrov, K. Egiazarian, and V. Katkovnik, “Hyperspectral terahertz pulse time-domain holography: noise filtering,” in *Digital Holography and Three-Dimensional Imaging*, (Optical Society of America, 2019), p. accepted.

**44. **L. Zhang, H. Zhong, C. Deng, C. Zhang, and Y. Zhao, “Terahertz wave reference-free phase imaging for identification of explosives,” Appl. Phys. Lett. **92**, 091117 (2008). [CrossRef]

**45. **N. S. Balbekin, M. S. Kulya, A. V. Belashov, A. Gorodetsky, and N. V. Petrov, “Increasing the resolution of the reconstructed image in terahertz pulse time-domain holography,” Sci. reports **9**, 180 (2019). [CrossRef]

**46. **“Spectrally resolved incoherent holography: 3D spatial and spectral imaging using a Mach–Zehnder radial-shearing interferometer,” Opt. Lett.39, 1857 (2014). [CrossRef]

**47. **D. Claus, G. Pedrini, D. Buchta, and W. Osten, “Accuracy enhanced and synthetic wavelength adjustable optical metrology via spectrally resolved digital holography,” J. Opt. Soc. Am. A **35**, 546 (2018). [CrossRef]

**48. **K. Itoh, T. Inoue, T. Yoshida, and Y. Ichioka, “Interferometric supermultispectral imaging,” Appl. Opt. **29**, 1625 (1990). [CrossRef] [PubMed]

**49. **M. Massaouti, C. Daskalaki, A. Gorodetsky, A. D. Koulouklidis, and S. Tzortzakis, “Detection of harmful residues in honey using terahertz time-domain spectroscopy,” Appl. spectroscopy **67**, 1264–1269 (2013). [CrossRef]

**50. **M. M. Nazarov, A. P. Shkurinov, E. Kuleshov, and V. V. Tuchin, “Terahertz time-domain spectroscopy of biological tissues,” Quantum Electron. **38**, 647 (2008). [CrossRef]

**51. **K. I. Zaytsev, K. G. Kudrin, V. E. Karasik, I. V. Reshetov, and S. O. Yurchenko, “In vivo terahertz spectroscopy of pigmentary skin nevi: Pilot study of non-invasive early diagnosis of dysplasia,” Appl. Phys. Lett. **106**,053702 (2015). [CrossRef]

**52. **K. I. Zaytsev, K. G. Kudrin, S. A. Koroleva, I. N. Fokina, S. I. Volodarskaya, E. V. Novitskaya, A. N. Perov, V. E. Karasik, and S. O. Yurchenko, “Medical diagnostics using terahertz pulsed spectroscopy,” in *Journal of Physics: Conference Series*, vol. 486 (IOP Publishing, 2014), p. 012014.

**53. **N. V. Petrov, V. G. Bespalov, and M. S. Kulya, “Terahertz pulse time-domain holography for studying of broadband beams propagation dynamics,” in *Digital Holography and Three-Dimensional Imaging*, (Optical Society of America, 2018), pp. DTu2F-7. [CrossRef]

**54. **M. S. Kulya, N. V. Petrov, A. N. Tcypkin, and V. G. Bespalov, “Influence of raster scan parameters on the image quality for the THz phase imaging in collimated beam with a wide aperture,” J. Physics: Conf. Ser. **536**, 012010 (2014).

**55. **V. Semenova, M. Kulya, N. Petrov, Y. V. Grachev, A. Tsypkin, S. Putilin, and V. Bespalov, “Amplitude-phase imaging of pulsed broadband terahertz vortex beams generated by spiral phase plate,” in Infrared, Millimeter, and Terahertz waves (IRMMW-THz), 2016 41st International Conference on, (IEEE, 2016), pp. 1–2.

**56. **W. L. Chan, H.-T. Chen, A. J. Taylor, I. Brener, M. J. Cich, and D. M. Mittleman, “A spatial light modulator for terahertz beams,” Appl. Phys. Lett. **94**, 213511 (2009). [CrossRef]

**57. **S. Withington, C. Tham, and G. Yassin, “Theoretical analysis of planar bolometric arrays for THz imaging systems,” in *Millimeter and Submillimeter Detectors for Astronomy*, vol. 4855 (International Society for Optics and Photonics, 2003), pp. 49–63. [CrossRef]

**58. **A. Koulouklidis, V. Y. Fedorov, and S. Tzortzakis, “Spectral bandwidth scaling laws and reconstruction of THz wave packets generated from two-color laser plasma filaments,” Phys. Rev. A **93**, 033844 (2016). [CrossRef]

**59. **Y. A. Kapoyko, A. A. Drozdov, S. A. Kozlov, and X.-C. Zhang, “Evolution of few-cycle pulses in nonlinear dispersive media: Velocity of the center of mass and root-mean-square duration,” Phys. Rev. A **94**, 033803 (2016). [CrossRef]

**60. **Y.-S. Lee, *Principles of terahertz science and technology*, vol. 170 (Springer Science & Business Media, 2009).

**61. **Y. Shen, T. Watanabe, D. Arena, C.-C. Kao, J. Murphy, T. Tsang, X. Wang, and G. Carr, “Nonlinear cross-phase modulation with intense single-cycle terahertz pulses,” Phys. Rev. Lett. **99**, 043901 (2007). [CrossRef] [PubMed]

**62. **A. Gürtler, C. Winnewisser, H. Helm, and P. U. Jepsen, “Terahertz pulse propagation in the near field and the far field,” JOSA A **17**, 74–83 (2000). [CrossRef] [PubMed]

**63. **Y. A. Shpolyanskiy, “Envelope, phase, and frequency of ultrabroadband signal in a transparent medium,” J. Exp. Theor. Phys. **111**, 557–566 (2010). [CrossRef]

**64. **A. Belashov, A. Gorodetsky, M. Kulya, and N. V. Petrov, “Stepwise approach to numerical simulation of broadband femtosecond pulses propagation through amplitude and phase objects,” in *Digital Holography and Three-Dimensional Imaging, OSA Technical Digest (online)*, (Optical Society of America, 2019), p. W1A.6.

**65. **K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-d transform-domain collaborative filtering,” IEEE Transactions on Image Process. **16**, 2080–2095 (2007). [CrossRef]

**66. **W. Wolfe, “Introduction to imaging spectrometers, tutorial text vol,” TT25 SPIEOpt. Eng. Press. Bellingham, Washington, USA (1997).

**67. **A. F. Goetz, “Three decades of hyperspectral remote sensing of the earth: A personal view,” Remote. Sens. Environ. **113**, S5–S16 (2009). [CrossRef]

**68. **M. Govender, K. Chetty, and H. Bulcock, “A review of hyperspectral remote sensing and its application in vegetation and water resource studies,” Water Sa33 (2007).

**69. **Y.-Z. Feng and D.-W. Sun, “Application of hyperspectral imaging in food safety inspection and control: a review,” Critical reviews in food science and nutrition **52**, 1039–1058 (2012). [CrossRef] [PubMed]

**70. **G. Lu and B. Fei, “Medical hyperspectral imaging: a review,” J. biomedical optics **19**, 010901 (2014). [CrossRef]

**71. **S. Kalenkov, G. Kalenkov, and A. Shtanko, “Hyperspectral holography: an alternative application of the Fourier transform spectrometer,” JOSA B **34**, B49–B55 (2017). [CrossRef]

**72. **G. S. Kalenkov, S. G. Kalenkov, I. G. Meerovich, A. E. Shtanko, and N. Y. Zaalishvili, “Hyperspectral holographic microscopy of bio-objects based on a modied Linnik interferometer,” Laser Phys. **29**, 016201 (2018). [CrossRef]

**73. **V. Katkovnik, M. Ponomarenko, and K. O. Egiazarian, “Complex-valued image denosing based on group-wise complex-domain sparsity,” CoRR abs/1711.00362 (2017).

**74. **K. Lee and J. Ahn, “Single-pixel coherent diffraction imaging,” Appl. Phys. Lett. **97**, 241101 (2010). [CrossRef]

**75. **J. X. X.-C. Zhang, *Introduction to Thz Wave Photonics* (Springer, New York, 2010). [CrossRef]

**76. **J. M. Bioucas-Dias and G. Valadao, “Phase unwrapping via graph cuts,” IEEE Transactions on Image process. **16**, 698–709 (2007). [CrossRef]

**77. **V. A. Semenova, M. S. Kulya, N. V. Petrov, Y. V. Grachev, A. N. Tsypkin, S. E. Putilin, and V. G. Bespalov, “Amplitude-phase imaging of pulsed broadband terahertz vortex beams generated by spiral phase plate,” in 2016 41st International Conference on Infrared, Millimeter, and Terahertz waves (IRMMW-THz), (2016), pp. 1–2.

**78. **J. F. Nye and M. V. Berry, “Dislocations in wave trains,” Proc. R. Soc. Lond. A **336**, 165–190 (1974). [CrossRef]

**79. **K. Falaggis, T. Kozacki, and M. Kujawinska, “Hybrid single-beam reconstruction technique for slow and fast varying wave fields,” Opt. Lett. **40**, 2509 (2015). [CrossRef] [PubMed]

**80. **A. Borot and F. Quéré, “Spatio-spectral metrology at focus of ultrashort lasers: a phase-retrieval approach,” Opt. express **26**, 26444–26461 (2018). [CrossRef] [PubMed]

**81. **G. Pariente, V. Gallet, A. Borot, O. Gobert, and F. Quéré, “Space–time characterization of ultra-intense femtosecond laser beams,” Nat. Photonics **10**, 547 (2016). [CrossRef]