Sweep Distortion Removal from THz Images via Blind Demodulation

Heavy sweep distortion induced by alignments and inter-reflections of layers of a sample is a major burden in recovering 2D and 3D information in time resolved spectral imaging. This problem cannot be addressed by conventional denoising and signal processing techniques as it heavily depends on the physics of the acquisition. Here we propose and implement an algorithmic framework based on low-rank matrix recovery and alternating minimization that exploits the forward model for THz acquisition. The method allows recovering the original signal in spite of the presence of temporal-spatial distortions. We address a blind-demodulation problem, where based on several observations of the sample texture modulated by an undesired sweep pattern, the two classes of signals are separated. The performance of the method is examined in both synthetic and experimental data, and the successful reconstructions are demonstrated. The proposed general scheme can be implemented to advance inspection and imaging applications in THz and other time-resolved sensing modalities.


I. Introduction
Due to fine time resolution and broad spectral coverage, terahertz time domain spectroscopy (THz-TDS) has become a leading method in THz spectroscopy [1,2], imaging [3] and nondestructive testing [4] of dielectric structures.There is an extensive literature on improving the THz imaging capability.In general, a large body of literature is focused on improving the hardware as THz-TDS power levels are usually at sub microwatt levels [5,6,7].This perspective can be extended to improving THz acquisition methodologies such as compressive [8] and wide field acquisitions [9,10,11].On the signal processing side the mainstream research has been focused on improvement of signal to noise ratio (SNR) [12,13].Many of these SNR improvement methods are developed for transmission mode spectroscopy and have appreciable alignment with infrared spectroscopy [14,15,16].However, when it comes to imaging, inspection, and content extraction of surfaces and layered structures, THz-TDS can suffer significantly from phase distortions along the sample surface.These sweeping distortions in THz time domain imaging appear as dominant challenge for in depth imaging and content extraction in densely layered structures, irregular surfaces or any 2D geometry.Unfortunately, these heavy distortions are directly induced by the nature of the sample and THz-TDS measurement scheme, and cannot be addressed by minor hardware improvement, conventional denoising or SNR enhancement techniques.
In this paper we propose and demonstrate a mathematical framework that enables demodulation of the recorded signal from sweep distoritions induced by slight depth variations or the layered structure inter-reflections.We view the problem as a demodulation of the distortion profiles from the sample texture, based on the reflected wave measured at different instances of time.The problem of interest is modeled as a bilinear inverse problem (BIP), which can be addressed using a low-rank matrix recovery framework.We propose an alternating minimization scheme, where a prior structure is considered for each factor in the BIP.More specifically, we assume the distortion profiles belong to a low dimensional subspace (which can be extracted from the raw data) and the layer structure is of binary nature.The model considered for the layer texture is in fact a shape-based approximation (see [17,18,19] for examples), a phase corresponding to the main texture and another phase representing the anomalies and inclusions.Using the given priors for each factor in the BIP, we alternatively solve the demodulation problem to exclude the undesired sweep distortions from the THz images.
Our mathematical presentation mainly relies on multidimensional calculus.We use bold characters to denote vectors and matrices.Considering a matrix A and the index sets Γ 1 , and Γ 2 , we use A Γ 1 ,: to denote the matrix obtained by restricting the rows of A to Γ 1 .Similarly, A :,Γ 2 denotes the restriction of A to the columns specified by Γ 2 , and A Γ 1 ,Γ 2 is the submatrix with the rows and columns restricted to Γ 1 and Γ 2 , respectively.
The remainder of this paper is organized as follows.In Section II we overview the physics of the problem and discuss the main demodulation problem to be addressed.In Section III we present the main algorithmic framework to address the demodulation problem using an alternating minimization scheme.Finally, in section IV the sensitivity and performance of the algorithm are assessed with synthetic data.We further experimentally demonstrate the blind demodulation of the data recorded from both single and multilayered structures and report the algorithm outcomes for experimental data.We conclude this section with some remarks and future directions of research.

II. Physics of the Problem
Consider an electric field that is linearly polarized in the x-direction, propagating along the z-direction in the free space.The waveform is considered to be a finite duration THz pulse χ(t), for which the traveling field along every point (x, y) is and c is the wave speed.Due to confocal nature of the measurement at each (x, y) position, the problem of interest can be considered as extracting the contents of a dielectric slab placed perpendicular to the wave propagation axis, based on analyzing the reflected field For a homogeneous layer of width d, with reflection coefficient ρ and refraction index n ρ , the returned signal can be analytically expressed by convolving the pulse χ(.) with a train of impulse functions with decaying coefficients.More specifically [20], where (2) and τ ρ = n ρ d/c.Each impulse term in (2) corresponds to a reflection.Particularly, the first term corresponds to the reflection from the front surface of the slab, the second term (m = 1) corresponds to the back surface reflection and the remaining terms correspond to the returned waves after a number of inter-reflections within the slab.
The response in (1) can still be reasonably accurate when the dielectric slab is homogeneous along the zaxis, i.e., ρ(x, y, z) = ρ 0 (x, y) for 0 ≤ z ≤ d.It only requires plugging the point-wise values of ρ and τ ρ in the formulation.
In the case of perfect or close to perfect reflection, where ρ ≈ 1, the first impulse term in (1) dominates all the other terms and the returned wave is a simple modulation of ρ with the pulse waveform.When 0 < ρ 1, the most dominant terms in u(τ) correspond to the first and second impulses, and the remaining terms decay exponentially fast.Specifically, in many practical scenarios where the pulse width is sufficiently small and the slab width d is large, the first reflected waveforms does not overlap with the subsequent reflections.In this case, still a modulation of ρ with the pulse waveform is observed in different time intervals.For multilayer dielectric slabs, the impulse response can still be cast as an impulse train, however deriving closed-form expressions for the coefficients is a sophisticated task and beyond the scope of this paper.When the reflection coefficients of the layers are small (hence transmission is the leading phenomenon), the dominant terms associated with different layers remain reasonably distinct from one another.If the pulse width is sufficiently small, for each layer the corresponding dominant reflection is still proportional to the reflection coefficient of that layer.
In all the situations stated above, for a measurement location (x, y, z 0 ) with a fixed z component, and sampling time t = t 0 , the reflected signal is ρ(x, y)u(t 0 + z 0 /c), which is a constant multiple of ρ(x, y).In other words, the reflected images should contain the same pattern as ρ.However, due to non-ideal configurations such as a tilted or uneven sample, the effective measurement points are (x, y, z 0 + ε(x, y)), where ε(x, y) is a function of small magnitude.For instance, in the case of a tilted sample ε(x, y) = α 1 x + α 2 y, where α 1 and α 2 are small constants.
Technically, the demodulation problem of interest in this paper corresponds to extracting the sample struc-ture ρ(x, y), based on observing the reflected signal E − 0 (x, y, z, t) at a fixed z coordinate and time samples t = t 1 , • • • , t M .Due to the non-ideal configurations stated above, the returned signal, sampled at a time t j , is in the form of ρ(x, y)u j (x, y) + n j (x, y), where u j (x, y) depends on the pulse waveform, sampling time and ε(x, y).The additive term n j (x, y) models the noise uncertainty and undesired electromagnetic interactions.We will refer to the u j (x, y) factor as a sweep distortion profile, which is often a slowly varying function when the misconfiguration is small.Specifically, inspired by the applications presented in this paper, we are interested in single or multilayer composite slabs of binary nature, where for each layer and ρ 0 and ρ 1 are constant.An example of a composite slab which roughly follows such model is a card marked with deep ink (see example in Figure 1(b)).
Figures 1(c), (d) show examples of the observed reflection from a binary slab (letter "M" printed on a card) at two different time instances.We are willing to make a binary characterization of the sample content by demodulating and excluding the multiplicative sweep distortion profiles affecting the desired image.In the sequel we will present an inversion scheme to characterize the slab contents based on such multiplicative observations.III.Methods

I. General Setup
For a more concise formulation we use the more compact spatial notation x = (x, y).Based on the discussions in the previous section, our observation is in the form of where y j (x) represents the j-th observed reflected image dependent on the planar coordinate x, ρ(x) is the true dielectric texture profile common across all the layers, u j (x) is the sweep distortion profile corrupting layer j, n j (x) models the noise and uncertainty and D is the domain of imaging for x.As a reasonable model, we consider a normal distribution for the noise, that is, n j (x) ∼ N (0, σ 2 ) are independent and identically distributed (i.i.d) random variables for j = 1, • • • , M and x ∈ D.
The ultimate goal is recovering the image ρ(x) and the distortion profile u j (x) based on the observation of y j (x).Our strategy to address the problem is considering known subspace models for the signals u j (x).More specifically, u j (x) ∈ S j , where S j = span s j 1 (x), • • • , s j N j (x) and N j is the j-th subspace dimension.In Section I.3 we introduce an algorithm to construct a set of low dimensional subspaces S 1 , • • • , S M from the set of observations y 1 (x), • • • y M (x).Also, based on the structure of the problem, we may benefit from prior assumptions about ρ(x) as will be detailed in the sequel.

I.1 Problem Reformulation as a Low-Rank Recovery Scheme
Consider discretizing the domain D into a collection of pixels {x i } P i=1 .Based on the model in ( 4), a natural way of inverting the data for the image and sweep distortion profiles is through the following minimization: Here denotes the Hadamard (pointwise) product and the bold scripts are vectors of length P, corresponding to the variables in (4).Consider S j to be a matrix representation of the subspace S j and Q is a subspace in which the signal ρ lives (in the case of no such assumption, Q can be simply the identity matrix spanning the canonical basis).Under this assumption, the minimization in (5) can be cast as min where the search for the unknown vectors is performed in the corresponding subspaces.Suppose T and P j to be a suitable selection matrix (subset of the rows of the identity matrix) such that α j = P j α.For the point-wise product of the image and the sweep distortion profiles we have In other words considering X = βα T to be the outerproduct of the unknown vectors β and α, the noise-free observations correspond to the inner-products of X with known matrices Q T i,: S j i,: P j .Basically, if we define a linear operator the noise-free observations are linearly dependent on X, while the original version of the problem in ( 6) is nonlinear in terms of β and α.However, since the outerproduct of two nonzero vectors is a rank-one matrix, an equivalent reformulation of ( 6) is the minimization where Y = [y 1 , • • • , y M ] and .F denotes the Frobenius norm.Despite the linear relationship between the unknown variable and the noise-free observations, the rank constraint in ( 9) makes it a combinatorial problem.As a remedy, a well-known approximation strategy is to use the nuclear norm of X, denoted by X * , as a convex proxy to rank(X): It is an interesting fact that under certain conditions on A (e.g., see [21]), the approximate solution of ( 10) can coincide with the solution of ( 9).Once a rank-one solution X * is available, the factors α and β can be determined up to a constant multiple (β = X * :,1 and α T = X * 1,: /X * 1,1 ).Other than the relaxation technique, there are other methods reported in the literature to address instances of (9).We are specifically interested in alternating minimization approach [22], since it allows incorporating prior information into the reconstruction of each bilinear factor.

I.2 Alternating Minimization: A Maximum a Posteriori Framework
In minimizing the non-convex program (5) a reliable technique would be to setup an alternating minimization scheme.Basically, by fixing any of the operands in the bilinear term in (5) (either ρ or u j ) the problem turns into a standard least squares problem.The process would be to initialize one of the operands in the bilinear model, say ρ; with this quantity fixed, the resulting least squares in terms of u j is solved.By plugging in the acquired solutions for u j 's, we can solve another least squares problem in terms of ρ, and proceed with such alternation until convergence.This alternation approach is capable of approximating the solution to (9).Again under certain incoherence conditions on A (slightly stricter than the ones stated for nuclear norm minimization [22,21]), the approximate solution coincides with the solution of (9).The major advantage of using an alternation scheme over the nuclear norm surrogate is the possibility of incorporating prior information beyond the subspace constraint on each factors of the problem in (4) (e.g., see [23]).
Basically, the proposed scheme corresponds to an alternation between the maximum likelihood (ML) estimates of ρ and u j .When some level of prior knowledge about the ρ and/or u j exists, the framework could be generalized to an alternation between the maximum a posteriori (MAP) estimates of ρ and u j , as will be detailed in the sequel.
Consider f (ρ, {u j } M j=1 {y j } M j=1 ) to be the joint probability density function (pdf) of ρ and u j 's given the measurements y j .The MAP estimates of ρ and u j correspond to the maximizing parameters of the posterior distribution: As a natural generalization, the alternating MAP estimation corresponds to the process sketched in Algorithm 1.
In the sequel we elaborate on how to implement the proposed scheme using certain facts about the problem.We would like to note that we make some reasonable assumptions and provide some slight modifications in implementing the steps outlined in Algorithm 1.

I.3 Iterative Reconstruction of the Sweep Distortion Profiles
In this section we mainly focus on addressing the first stage of Algorithm 1 (line 3), which is determining the distortion profiles u j based on the observations {y j } M i=1 and the knowledge of ρ.
To address this problem we only assume that based on the physics of the problem, the signals u j reside in a low dimensional subspace.In other words, u j = S j α j for the discrete representation.We however make no prior assumptions about the coefficient vector α j .Since ρ u j = diag(ρ)u j , we simply consider solving the least squares problems of the form The solution to (12) is where A † denotes the pseudoinverse of A. In a sense, by making no assumptions about α j , instead of performing a MAP estimation (as outlined in Algorithm 1), for a given ρ (k) we are performing an ML estimation of the sweep distortion profiles by solving the least squares problems (12).The choice of the subspace plays a key role for this problem.In the remainder of this subsection we present a one-time process which generates embedding subspaces S j based on the observations y j .
From a technical standpoint, y j (x) ≈ ρ(x)u j (x), where as stated earlier, u j is a rather smooth function and for our application ρ(x) is a function of almost binary nature.Therefore, if we have a multi-scale representation of y j , there should be a good overlap between u j and the course-level approximation of y j .The high frequency components of y j should be mainly in hold of ρ and the noise.
Based on this argument, consider a wavelet representation of the signal y j as where φ 0 ,w are the scaling functions, ψ ,w are the wavelet basis and 0 is a fixed scaling level [24].For a fixed subspace order N j , our proposed subspace selection for u j corresponds to the top N j scaling/wavelet basis functions with the largest coefficients (in magnitude).In other words, we sort the wavelet coefficients in descending magnitude order and select the basis functions associated with the top N j coefficients.The selected basis (in discrete form) correspond to the columns of S j .This process is performed only once during the entire reconstruction and once the subspaces are determined, they remain fixed throughout the alternation scheme.

I.4 Iterative Reconstruction of the Image Profile
Inspired by the second stage proposed in Algorithm 1, the main objective in this section is finding the MAP estimate for ρ(x).In order to maintain a more compact formulation and closed-form expressions we make some independence assumptions about the pixel values as will be detailed in the sequel.
Based on the binary model considered, we assume that each pixel in the image ρ belongs to one of the two classes 0 or 1, and denote the pixel class by C(x).For our estimation problem, we are willing to address the maximization problem In other words, given the observation y j and the distortion profiles u (k) j , we are willing to acquire MAP estimates of the image profile value at each pixel along with the binary class that pixel belongs to.
Generally speaking, for arbitrary random variables ρ, y, u and C (discrete), using the Bayes' rule we have where we used f (C, u) = P(C|u) f (u) to derive the second equality.Based on this result we have where for the sake of convenience, in all the expressions we dropped the prior knowledge of {u Basically, in (13), the vectors {u (k) j } M j=1 are known and deterministic and our intention to derive ( 14) is to make the MAP estimation in terms of probability density functions which all make such prior assumption.Also since the denominator term in ( 14) is independent of ρ and C, in (15) we neglected the corresponding term as it will be a constant factor in the maximization.
To further proceed with simplifying the MAP objective, we assume that the pixel values in ρ are independent of each other.Based on this assumption and the i.i.d nature of the noise we can state that where y j,i is the i-th element of y j .The point-wise distributions appeared in ( 16) can now be modeled based on the problem setup.
Clearly, based on (4), f (y ), which is an immediate result of the normal noise model.We also assume that where p 0 and p 1 = 1 − p 0 are somehow the prior estimates of the portion of pixels belonging to each class (e.g., knowing that roughly 20% of the image pixels correspond to the anomaly).If no such information is available a priori, a reasonable assumption is p 0 = p 1 = 0.5.Finally, with reference to the term f (ρ i |C i ), knowing that a pixel belongs to class 0 (or 1) we know that its value is likely to be close to ρ 0 (or ρ 1 ).In practice f (ρ i |C i ) should model the uncertainty on how the values of each class concentrate around the mean class value.While we can use different distributions to model such concentration, for simplicity and in order to obtain closed form expressions, we assume that f (ρ i |C i = 0) and f (ρ i |C i = 1) are in the form of truncated normal distributions, taking positive arguments and mainly concentrating around ρ 0 and ρ 1 , respectively.By definition, a random variable truncated to values greater than zero takes the following pdf: Here, μ and σ2 are roughly the mean and variance and γ is a normalizing factor to assure the pdf integrates to 1.To model an almost binary nature of the pixel values, we assume f (ρ Figure 2 shows the underlying truncated distributions.The positivity assumption in ( 17) is imposed by the physics of the problem, where the reflection coefficient of the slab is always considered to be a positive quantity.
Having the constituting terms modeled in ( 16), we can proceed with the pixel value and class label MAP estimation.To maximize G(ρ, C) we may minimize the negative log function Based on the proposed distribution models, one can easily verify that ρ 0 0 where and K is a constant in terms of σ, γ i and the constant factor in (15).The minimization objective g(ρ, C) is separable in terms of the cost for each pixel and we can minimize each term individually.Also minimizing g(ρ i , C i ) can be performed in a serial manner by first minimizing with respect to ρ i and then with respect to the class label C i .More concretely, min For a fixed given label C i = C, we can minimize g(ρ i , C) with respect to ρ i by setting ∂g(ρ i , C)/∂ρ i to zero.Based on the formulation in (18) a simple calculation yields arg min The function g(ρ i , C) is convex in terms of ρ i , and the solution to the constrained problem (ρ i ≥ 0) is the expression in (20) if feasible, or 0 otherwise.Basically, We are now only left with minimizing g(ρ * i , C i ) with respect to C i .Since the class labels C i only takes binary values, we can conveniently find the corresponding minimizer through the following binary comparison: .
Equation (20) in some sense generates the optimal solution ρ i , but leaves us with an ambiguity about the class label.The label can be obtained based on the above comparison, to find the ultimate MAP estimator for the image profile.Algorithm 2 summarizes the steps sketched in this section, and presents the overall process to perform the demodulation task.
Algorithm 2 Decoupling algorithm As noted in Section I.3 3: ρ ← 1 4: repeat 5: u j ← S j α j 8: , 0 10: if g(w 0 , 0) ≤ g(w 1 , 1) then g(., .) is defined in (18) 12: ρ i ← w 0 13: else 14: While the parameters σ 0 and σ 1 have statistical interpretations, from an optimization standpoint they can be considered as free parameters controlling the algorithm's performance.When σ 0 and σ 1 are set to be small quantities, the algorithm converges faster at the expense of more sensitivity to the initialization (possibility of recovering a local minimizer).Assigning relatively larger values to these quantities paves the path towards identifying the global minimizer in more number of iterations.This is also a more reliable parameter selection in the case of noisy observations.

IV. Simulation and Experiments
In this section we assess the performance of the proposed technique in demodulating simulated and real data.The simulation results mainly highlight the performance of the method for various number of frames and noise levels in the data.The second set of experiments demonstrate the method's success in removing multiplicative sweep distortion profiles from actual THz measurements.

I. Simulated Data
We simulate the reflected signal from a dielectric slab with the x-y profile depicted in Figure 3(a).The dielectric width is d = 100 µm.For the binary dielectric slab the reflection coefficients are ρ 0 = 0.3 and ρ 1 = 0.1.The sample is probed with a bipolar THz pulse (simply derivative of a Gaussian) as where t 0 = 1 ps and T = t 0 /4.To model a sweep modulated signal, we assume that instead of point-wise observations at a constant z = z 0 , the measurements are performed at z = z 0 + α 1 x + α 2 y, where (x, y) is the planar coordinate, and α 1 = 10 −6 , α 2 = 10 −4 are small constants.
We take M = 10 uniform samples of the reflected wave in a time period of 0.8 ps.Among the recorded frames Figure 3(b) shows three sample images, where an almost vertically moving sweep distortion across the frames is observable.The samples are synthetically corrupted with white noise (SNR = 10 dB).
In recovering the binary profile, for the subspace construction step associated with this example we use N 1 = • • • = N 10 = 100 most dominant wavelet basis.Consistently through this example and the remaining experiments we use the symlet wavelet family.
Figure 3(c.1)reports the reconstructed reflection profile ρ, using our proposed scheme.The algorithm converges in less than 20 iterations, and as may be observed, the result is no more in hold of the distortion traces.A simple rounding of the result to the closest binary value (between ρ 0 and ρ 1 ) yields Figure 3(c.2),which is reasonably close to the original profile depicted in Figure 3(a).
Figure 3(d.1)reports the reconstruction using the nuclear norm framework proposed in Section I.1.While some level of distinction between the two binary phases is achieved, the fact that our proposed scheme efficiently uses the binary prior helps us outperform the nuclear norm reconstruction.A similar binary approximation of the nuclear norm reconstruction is demonstrated in Figure 3(d.2),which clearly fails to characterize some of the main details in the original binary profile.To more extensively analyze the performance of our algorithm, we proceed by the sensitivity analysis of the main parameters affecting the reconstruction.For this purpose, two main scenarios are considered.First, an ideal case where the exact sweep distortion subspace is known a priori.From a theoretical standpoint, one way to characterize a subspace for the sweep distortion profiles is to run exactly the same experiment with a For each MSE report M frames are randomly selected from the 20 frames and passed to the algorithm.The process is performed multiple times for each M and the average MSE is reported homogeneous slab and take the observations as the subspace basis.While experimentally such accurate characterization of the subspace is hard, here we present it as a hypothetical case for comparison purposes and to compare the results against an ideal setup.
For the second class of experiments, we use the proposed wavelet thresholding framework to determine the distortion subspaces.For each proposed scenario we test the algorithm for various values of M (number of available frames) and the SNR in the observations.Both the reconstructed ρ as well as a binary rounded version are compared with the reference image depicted in Figure 3(a) and the mean squared error (MSE) is reported.
Figure 4 shows the MSE values for various number of available frames and a fixed SNR value of 10 dB.The dark solid curve corresponds to the case of exact subspace characterization and the dashed curve standing close to it is the MSE value for the reconstructed ρ rounded to the closest binary value.We can see that for only 7 frames (or more), a perfect or close to perfect recovery of the binary profile is possible.For fewer number of available frames, the MSE still remains at a controlled level.
The lighter solid line (and the dashed line standing close to it) report the MSE values of similar experiments when the sweep distortion subspaces are characterized through the wavelet thresholding.We can see that aside from a slight performance gap, the MSE values follow a similar reduction pattern as the ideal case.Based on the reported MSE values, the recovery is still close to the reference image.The slight performance gap between this case and the ideal case is the result of error in exact characterization of the sweep distortion subspaces.
Figure 5 shows the MSE values for the case of fixed number of frames (M = 10) and varying observation noise.For this experiment, in the case of known distortion subspace, a perfect recovery of the sweep distortion profile is possible for SNR values approximately exceeding 10 dB (by rounding to the closest binary value).In the case of characterizing the sweep distortion subspace through wavelet thresholding, the MSE values follow a similar trend as the case of knowing the distortion subspace a priori.The only difference is a slight performance gap which is again linked to the error in exact characterization of the sweep distortion subspaces.

II. Experimental Data
To experimentally apply the demodulation process we used a standard THz-TDS system in reflection geometry.In the first experiment (Figure 6) we used a binary sample with metallic surface and a cross shaped pattern carved into a 15 mm × 15 mm metal surface with 2 mm depth.As noted in Figure 6(a1-3), the spatial sweep distortions are dominantly present in the time domain observations.For this experiment M = 25 frames were available, three of which are shown here.
Figures 6(b1-3) and 6(c) report the decoupled sweep distortions from the binary cross profile.We can observe a good characterization of the cross profile thanks to the large number of available frames and high SNR observations.A more challenging experiment with less number of available frames and noisier data is presented next.A more challenging experiment uses a three layer paper sample with letters M-I-T printed on each page from front to back respectively.The page dimensions are 15 mm × 30 mm and the thicknesses are 300 µm, stacked flush next to each other (∼ 30 µm gap) similar to the structure of a book.Figure 7(a1-3) shows the reflected signal from the first layer at three different instances of time.Figures 7(b1-3) and 7(c1-3) show reflection samples from the second and third layers respectively.The available number of frames corresponding to the first, second and third layer are M = 25, 6 , 8, respectively.
Followed by the observations in each column of Figure 7, the demodulation results are presented for the three letters at different pages of the sample, as the THz pulse travels deeper into the pages.The rows 4-6 in Figure 7 correspond to the sweep distortion profiles decoupled from the observed data.The last row corresponds to the recovered character profiles.
For this multilayer setup, after the first layer the SNR drops and a larger number of inter-reflections is induced, thus the recovered letters "I" and "T" in panels (b.7) and (c.7) have more noise compared to the recovered letter "M" on the first page (panel (a.7)).
The clean results from the first experiment (Figure 6) on the metallic surface are anticipated due to high SNR and perfectly flat surface of the polished steel, however for the paper the reflection is rather small (<0.15), the distortion contrast is at the same level of the signal contrast It is worth noting that the major burden in character-izing the sample inhomogeneity profiles (characters in this example) is the varying signal level in the observed images.For instance, while the contrast between the character "M" and the background can be visually detected in Figure 7(a.3), the pixel values drastically vary from one portion of the letter to the other.As a result, characterizing the letter based on the pixel values becomes an inconceivable task.However, once a binary profile is reconstructed through the proposed algorithm, other post-processing schemes can be employed to characterize the reconstructions, especially in the case of the more noisy reconstructions such as Figure 7(c.7).The interested reader is referred to [25,26,27,28], where advanced shape composition techniques are employed to accurately identify binary inclusions in noisy images.This class of techniques can accurately characterize the inclusion in an image like Figure 7(c.7),even when the noise is higher, parts of the object are missing due to signal occlusions, and we have clutter or overlapping characters caused by the inter-reflections from neighboring layers.

III. Concluding Remarks
The main outcome of this research is developing a general demodulation scheme to process the reflected signals, conveniently applicable to THz imaging.While a careful processing of the THz data can produce high resolution images thanks to its high frequency, dealing with phenomena such as inter-reflections, noise and modulated distortions is inevitable.We showed that reformulation of the problem as a bilinear inverse problem and making practical assumptions, such as the binary prior, can assist us with a successful inversion of THz measurements in THz-TDS systems.
Using basic statistical assumptions about the image and the observation, we were able to develop an iterative inversion scheme, which is computationally efficient, distributable and scalable to big data sets.In fact, the main computational load at every step of the algorithm is a least squares solve.The remaining computations are in the form of thresholding-type operations.
Since in many applications such as water profilometry, extracting coded structures, extracting sub surface cracks, the main objective is a binary characterization of the inhomogeneity in the samples, we used a binary prior in our modeling.The algorithm can be naturally generalized to the case of polyadic prior, where more than two levels are considered for the image profile.Basically, the proposed algorithm can be further modified or combined with post or pre-processing tools to handle a larger class of imaging applications.

Figure 1 :
Figure 1: (a) A dielectric slab, emitted and reflected electrical fields; (b) setup schematics, blue is the emitted THz field and the red waveform is indicative of a typical reflected signal in one pixel; (c,d) induction of sweep distortion: observation of the returned field at two different time instances Following the categorization in (3), the pixel values in class 0 concentrate around ρ 0 and the pixel values associated with class 1 concentrate around ρ 1 .Consider ρ = [ρ 1 , • • • , ρ P ] T and C = [C 1 , • • • , C P ] T to be vectors containing the pixel values and the class labels for the discrete image.

Figure 4 :
Figure 4: MSE as a function of available number of frames, M,for various scenarios: a total of 20 frames using uniform sampling in time is made available through the simulation.For each MSE report M frames are randomly selected from the 20 frames and passed to the algorithm.The process is performed multiple times for each M and the average MSE is reported

Figure 5 :
Figure 5: MSE as a function of SNR for a fixed number of available frames (M = 10)

Figure 6 :
Figure 6: Experimental demonstration of sweep decoupling for a metallic surface; (a.1-3) three time instances of the recorded image from the sample, vertical sweep distortions are dominantly over the cross pattern; (b.1-3) the recovered distortion profiles corresponding to the observations shown; (c) the recovered binary profile