Open Access
Add to BibSonomy
Abstract
Terahertz sparse deconvolution based on an iterative shrinkage and thresholding algorithm (ISTA) has been used to characterize multilayered structures with resolution equivalent to or finer than the sampling period of the measurement. However, this method was only studied on thin samples to separate the overlapped echos that can’t be distinguished by other deconvolution algorithms. Besides, ISTA heavily depends on the convolution matrix consisting of delayed incident pulse, which is difficult to precisely extricate from the reference signal, and thereby fluctuations caused by noise are occasionally treated as echos. In this work, a terahertz sparse deconvolution based on a learned iterative shrinkage and thresholding algorithm (LISTA) is proposed. The method enclosed the matrix multiplication and soft thresholding in a block and cascaded multiple blocks together to form a deep network. The convolution matrices of the network were updated by stochastic gradient descent to minimize the distance between the output sparse vector and the optimal sparse representation of the signal, and subsequently the trained network made more precise estimation of the echos than ISTA. Additionally, LISTA is notably faster than ISTA, which is important for real-time tomographic-image processing. The algorithm was evaluated on terahertz tomographic imaging of a high-density poly ethylene (HDPE) sample, revealing obvious improvements in detecting defects of different sizes and depths. This technique has potential usage in nondestructive testings of thick samples, where echos reflected by minor defects are not discernible by existed deconvolution algorithms.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The pulsed terahertz (THz) radiation generated by photo-conductive antenna is versatile in characterizing substances due to its wide frequency band (0.1-10 THz) and unique wavelength range ($30 \mathrm{\mu} {\textrm{m}}$- $3 {\textrm{mm}}$) between those of microwave and infrared, which allows it to penetrate numerous nonmetallic materials that are opaque for infrared and visible light [1]. Besides, the low photon energy makes THz time-domain measurement a safe substitute or supplement to traditional testing approaches [2]. A THz transient consists of one or two wave cycles and is reflected due to the discontinuity of the Fresnel coefficients between interfaces and thus carries the depth information of the inner boundaries. This property makes reflected THz pulses ideal for imaging multi-layer structures and bulks with inside voids and defects [3,4]. For example, reflected THz imaging has been used to identify the stratigraphy and hidden features of classic oil paintings [3,5,6], the quality of automobile paints [7,8], and the coating uniformity of pharmaceutical tablets [9,10]. The heat insulation foam of the space shuttle with inner voids and defects that can’t be resolved by other techniques such as X-ray was successfully characterized [4,11]. For each point of the tomographic scan, the reflected time-domain signal comprises multiple time-shifted and amplitude-attenuated pulses, which are formed by the temporal waveform of the incident pulse convoluted with the impulse responses of the interfaces [3], therefore the impulse-response function of the system can be estimated by deconvoluting the time signal by the temporal waveform of the incident pulse. In practice, the echoes reflected by the interfaces overlap if the layer or defect is too thin, and the echoes may submerge by noise in poor detection, therefore another effect of deconvolution is to discern individual impulse responses [3].
Numerical deconvolution method for terahertz time-domain spectroscopy (THz-TDS) was first investigated to map the distribution of layers for single-layer and multi-layer paint films [12], by which the minimum resolvable thickness of a THz paint meter was determined by the width of the echo. A method based on multiple regression analysis was later proposed to improve the resolution [13], however, this method is restricted by the signal-to-noise ratio (SNR), which is greater than 100 in the literature, and the proximity of the regression model to the real time signal, and hence is not practical in mitigating the environmental artifacts, such as air turbulence or mechanical vibration, in real-life scenarios. Other researches attempted to realize higher resolution of thickness measurement by compressing the pulse width of the pump laser to generate narrower THz pulses for enhanced resolution [14], but this technique is still confined by the physical limit of the pump source. Subsequently, deconvolution focusing on filtering out the noise spikes in the higher frequency region of the impulse-response function by thresholding was investigated and successfully reduce the minimum resolvable thickness to the size of the sampling interval, which is substantially smaller than the pulse width [15]. Further study was proposed by Chen et al. to signify the impulse responses by combining Wiener deconvolution and wavelet shrinkage to eliminate the low-SNR region of the impulse-response function to better resolve the overlapped pulses and low-SNR echoes [8]. This approach was named hybrid frequency-wavelet domain deconvolution (FWDD) and sufficiently improved the estimation of impulse-response functions in processing the reflected THz signals from objects such as human palms [8], oil paintings [5], and polymer coatings [16]. However, the impulse responses retrieved by FWDD are broadened by the cutoff of the low-SNR region, which limits its resolvability for overlapped echoes [16].
The reflected signals from structures with discontinued Fresnel coefficients typically comprise a limited number of echoes, and the corresponding impulse-response functions have sparse representations, which means only a limited number of data points have nonzero values. This feature enables us to retrieve impulse-response function by sparse deconvolution. Compared with conventional deconvolution, sparse deconvolution estimates the impulse-response function by adjusting the sparse vector in time domain and therefore doesn’t introduce spikes to the high-frequency region, and thus avoids the problems of the methods depending on frequency-domain analysis and enables higher resolvability for overlapped echoes. The sparse deconvoluton also simplifies the identification of impulse responses because the background values are mostly zeros. Sparse deconvolution in THz has been investigated by previous researches to resolve thin stratigraphic structures [3,17], but to the best of our knowledge hasn’t been used to discern defects in thick materials. Besides, the conventional sparse deconvolution was addressed with iterative shrinkage and thresholding algorithm (ISTA), where a $l_1$ regularized optimization problem was solved by adjusting the components of the sparse vector with a fixed step size. Normally, hundreds of iterations were required to make the penalised least-squared criteria converge [18], which was excessively time-consuming for tomographic imaging [19].
In this paper, the application of learned iterative shrinkage and thresholding algorithm (LISTA) in THz reflection imaging to detect echoes that are compromised by system noise is discussed. Compared with other established deconvolution methods, LISTA demonstrates better performance in resolving low-SNR echoes and suppressing shot noise, which is highly important in processing the signal obtained by reflective scans of thick samples. To the best of our knowledge, this study is the first report of LISTA in respect to THz applications. The principle of LISTA is to estimate the sparse representation of the impulse-response function by a deep network, in which each layer of the network is analogous to an ISTA iteration [20]. According to previous studies and verified by our work, the number of layers in LISTA is an order of magnitude smaller than the number of iterations required for ISTA, which gives a notable improvement in computational time [20,21].
2. Theory and equipment
2.1 Experimental setup
As shown in Fig. 1(a), the experiment apparatus comprises a 2 mm thick silica wafer as the beam splitter to direct the reflected THz beam from the sample to the receiver, three plano-convex polymethylpentene (TPX) lenses to collimate and focus the beam, emitter and receiver photo-conductive antennas (PCA, Menlo Systems, Martinsried, Germany), and a all-fiber THz spectrometer (Quenda Terahertz Technologies, Qingdao, China). The THz spectrometer consists of a femtosecond laser with a 1560 nm central wavelength and a 100 MHz repetition frequency, an optical delay line, and a lock-in amplifier connected to a computer to process the received signal. The source laser pumps the emitter PCA to generate THz radiation, and the THz pulses accepted by the receiver PCA are probed with laser pulses that are splitted from the source and delayed by the optical delay line. In this study, the sampled time signal has a step size of 0.01 ps and a time span of 90 ps. The PCAs are connected to the spectrometer by 1560 nm single-mode optical fiber, and the width of the THz pulse is approximately 1 ps. As illustrated in Fig. 1(b), the sample is a high-density polyethylene wafer with prebuilt circular tunnels of diameters ranging from 0.8 mm to 5 mm, and 5 tunnels of the same size and increasing depths are aligned as a group. The tunnel closest to the top surface (highlighted in Fig. 1(b)) was scanned and the reflected signal was analyzed to elaborate the improved resolvablility of LISTA, which will be discussed in the next section. The sample was supported by a mechanical stage, which moved by 0.1 mm steps in both x and y axes.
Fig. 1. (a)The experimental setup.(b)The HDPE sample with prebuilt circular tunnels of different sizes and depths, as highlighted by the red dashed circles, and tunnels of diameters 1.0 mm and 0.8 mm are not shown in the image.
2.2 Iterative shrinkage and thresholding
In the time domain, the reflected THz signal $y\left (t\right )$ is the convolution of the incident THz pulse $h\left (t\right )$ with the impulse-response function $f\left (t\right )$, which represents the structural properties of the point of scan.
In this study, the reference signal is the average of signals received from 10 s repetitive acquisitions of the beam reflected by a gold coated mirror. The main pulse of the reference signal represents the incident pulse and is separated by a time window to describe $h\left (t\right )$. According to Fresnel’s law, the reflections by the gold coated mirror and the silica wafer give phase shifts of $\pi$ respectively, therefore $h\left (t\right )$ has the same phase as the incident pulse. Because the refractive index of HDPE is mostly constant from 0.1 THz to 3THz, which contains most of the frequency components of the pulse, the change of the temporal profile caused by dispersion is ignored during the propagation [22]. Therefore, the main echo of the reference signal is suitable to describe the reflected pulsed from the interfaces inside the HDPE sample. For materials with non-ignorable variation of refractive index, a correction factor accounting for the temporal spreading should be convolved with $h\left (t\right )$ [3]. Regrading the noise originating from the measurement, the discrete format of Eq. (1) follows,
where y is the discrete form of $y\left (t\right )$, H is the convolution matrix whose row are delayed versions of the reversed vector of $\textbf {h}^\intercal$, which is the transpose of the discrete form of $h\left (t\right )$ [18], f is the discrete form of $f\left (t\right )$, e accounts for the measurement noise.The sparse deconvolution aims to approximate the received THz signal y with Hf, and due to the aforementioned sparsity of f, Eq. (2) can be solved by $l_1$ regularized optimization [3], which is defined as,
where $||\textbf {H}\textbf {f} - \textbf {y}||_2^2$ denotes $\sum _{i}^{N}(\textbf {H}f_i - y_i)^2$, $||\textbf {f}||_1$ denotes $\sum _{i}^{N}|f_i|$, $N$ is the length of the vector, and $\lambda$ is the regularization parameter. According to the properties of the penalty function, where $\sigma$ is the standard deviation of the system noise, which is the standard deviation of the first 60 sampling points of the time signal in this study, and $||\textbf {h}||_2 = \sqrt {\sum _{n}|h(n)|^2}$. The detailed derivation of $\lambda$ is discussed in Ref. [22,23].ISTA is a very effective method to solve the regularized optimization problems, which leads to faster converge of the optimizer compared with the classical solvers [24]. In ISTA, each iteration involves matrix multiplication regarding H and $\textbf {H}^\intercal$ followed by a soft-thresholding step. The general iterative procedure is given by,
For the time signal of 9000 sampling points, Eq. (5) is iterated for 200 times to estimate f. As shown in Table 1, the average computational time of ISTA is 23.87 s. For a raster scan over a 1 cm by 1 cm area with a 0.1 mm step size, the number of pixels is 10,000, which takes ISTA 4 hours to analyze by parallelizing the computation into the 16 threads of a 8-core Intel i7 processor. Therefore, an improvement of ISTA towards faster speed is required by practical applications.
Table 1. Evaluation of ISTA and LISTA network of different depths.
2.3 Learned iterative shrinkage and thresholding
Gregor and LeCun developed a technique called algorithm unrolling that generalized iterative algorithm to neural network architecture and proved that the number of layers in a deep network was much smaller than number of iterations required in an iterative algorithm [21]. The unrolled network inherits the interpretability of the iterative algorithm and requires less training data than conventional neural networks. By determining the parameters through end-to-end training, the unrolled network achieves better performance than the iterative algorithm, where the parameters are defined by cross-validation and analytical derivations [20]. The unrolled counterpart of ISTA is named learned iterative shrinkage and thresholding algorithm (LISTA), where each layer of the network represents an iteration of ISTA, and numerous layers are cascaded into a network, as shown in Fig. 2(a).
Fig. 2. (a)An overview of the LISTA network.(b)The training history of networks consisting of 1-4 layers. (c)The outputs of the trained networks consisting of 1-4 layers, the input is the time signal produced by a measurement of the 1 mm tunnel.
As described in Fig. 2(a), the time signal y is given to each layer of the network, besides the initial value of f is fed to the first layer, then the output vector is passed to the subsequent layer, and the procedure is repeated until the end of the network. The weight parameters $\mathbf {H_{t}}$ and $\mathbf {H_{e}}$ are derived from Eq. (5) to generalize the original parametrization [20], which follows $\mathbf {H_{t}} = \mathbf {I} - \tau \mathbf {H}^\intercal \mathbf {H}$ and $\mathbf {H_{e}} = \tau \mathbf {W}^\intercal$. The shrinkage modality is the same as defined in Eq. (7). During the training, parameters $\mathbf {H_{t}}$, $\mathbf {H_{e}}$, and $\lambda$ are optimized, and the trained network achieves higher efficiency compared with ISTA even when the waveform of the incident beam is not precisely defined, i.e., the extrication of the incident pulse is affected by the noise level of the signal and the symmetry of the pulse, where equal size windows extending from the maximum and minimum of the pulse may not enclose the entire temporal waveform. f is initialized as 0, and the initial values of $\mathbf {H_{t}}$ and $\mathbf {H_{e}}$ are calculated using H defined in Eq. (2). The training uses stochastic gradient descent (SGD) scheme with an initial learning rate of $10^{-3}$, a decay rate of $10^{-6}$, a momentum of 0.9, and Nesterov accelerated gradient for faster convergence [25]. About 500 measurements are made for each tunnel, and $80\%$ of data is assigned for training and the rest is left for testing. The training runs for a maximum of 1000 epochs with a batch size of 16 and stops early if the loss doesn’t decrease for consecutive 10 epochs. Mean squared error (MSE) is used to estimate the loss, which follows
where $f_{i}$ denotes the target value of f, which will be discussed below, and $\hat {f}_{i}$ is the estimated value of f. As shown in Fig. 2(b), with the increase of layers, the network takes less epochs to converge, as deeper networks that are more capable of handling complex patterns converge by less epochs to configure the signal that can be tackled by shallower networks [26]. To compare the goodness of fit of the predictive models, R-squared criteria is used in Table 1. R-squared follows $\sum _{i}^{N}(\hat {f}_i - \overline {\mathbf {f}})^2 / \sum _{i}^{N}(f_i - \overline {f})^2$, where $\overline {\mathbf {f}}$ represents the average of all the values in f. R-squared ranges from 0 to 1, where 1 represents the perfect fit to the target, and 0 denotes that the fit complies with the mean of the target. As shown in Table 1 and Fig. 2(c), the networks consisting of 1 to 4 layers generate similar R-squared values and impulse-response functions, however, the time consumption grows linearly from 1 to 4 layers, regrading the efficiency, 1-layer network is used for the following discussion. Because all the parameters in H are calibrated by back propagation during the training to minimize the distance between f and the target, the parameters are adjusted to oppose the anomalous fluctuations in y that confuse ISTA to generate false impulse responses, subsequently the 1-layer LISTA uses the trained convolution matrix to retrieve the actual impulse responses and realize the sparsity of f. Compared with LISTA, H in ISTA is constructed by delayed versions of h and doesn’t adapt to the anomalous fluctuations in the signal, therefore the false impulse responses are manifested through the iterations. The higher efficiency of LISTA is shown by Table 1, where the time consumption of LISTA is more than one magnitude smaller than the time cost by ISTA to converge.To make the target value of f, the time signals acquired from each tunnel are averaged, and ISTA is applied to the averaged signal to make rough estimation of the impulse-response function. Comparing the averaged signal with the result of ISTA, the spikes that represent the actual impulse responses are selected manually, and the rest of the impulse-response function is set to zero, as shown in Fig. 3. According to the tunnel diameter, each item of the training data is assigned to the corresponding target, then a moving average of 300 points is applied to estimate the low-frequency baseline of the signal, which is later subtracted from the original data.
Fig. 3. (a)The time signal averaged from 546 measurements of the 1 mm tunnel. Labels 1-7 are: 1. pulse transmitted by the fast axis of the optical fiber and then reflected by the top surface of the sample and the front surface of the beam splitter, 2. pulse transmitted by the slow axis of the optical fiber and then reflected the top surface of the sample and the front surface of the beam splitter, 3. pulse reflected by the top surface of the tunnel and the front surface of the beam splitter, 4. pulse reflected by the bottom surface of the tunnel and the front surface of the beam splitter, 5. the secondary reflection of 1 (reflected by the rear surface of the beam splitter), 6. the secondary reflection of 2, 7. the secondary reflection of 3, 8. the secondary reflection of 4. (b)The result of ISTA. To make target data, the labeled spikes are kept and the rest of the vector is set to zero.
3. Results and discussion
Signal reflected by the 3 mm tunnel is selected from the testing data to illustrate the improvements made by LISTA. In the unprocessed data, as shown in Fig. 4(a), echoes labeled by 5, 7, and 8 are barely resolvable, and the reflections from the top and the bottom surfaces of the tunnel, 7 and 8 are just discernible from the maximums of the pulses.
Fig. 4. (a) Time signal obtained from a single measurement of the 3 mm tunnel and the impulse-response functions obtained by (b) FWDD, (c) ISTA, and (d) LISTA. The labeled echoes are the same as defined in Fig. 3(a). The deconvolved echoes from the top surface (e) and the bottom surface (f) of the tunnel, where $f_0$ represents the system impulse response. The results of FWDD, ISTA and LISTA are respectively depicted by the black, green and red lines. e. The frequency spectra obtained by making Fourier transforms to the original signal, $H$, and the impulse-response functions achieved by FWDD, $H_{FWDD}$, ISTA, $H_{ISTA}$, and LISTA, $H_{LISTA}$.
As shown in Fig. 4(b)-d, the signal is then deconvolved with FWDD, ISTA, and LISTA. The FWDD algorithm uses the same parameters as defined in Ref. [5], where the cutoff frequency $f_c = 4 {\textrm{THz}}$ to ensure the deconvolved signal has both a high resolution and a satisfactory SNR, and stationary wavelet denoising is implemented with a maximum level of 7 to eliminate the residual noise [3]. The deconvolution result of FWDD in Fig. 4(b) sufficiently improves the SNR of the signal, despite the high-frequency fluctuation at the beginning and the end, all the echoes are clearly resolved. Theoretically, the impulse-response function follows
An algorithm was then developed to locate the top and bottom surfaces of the tunnels. By finding the peaks in the time signal and regarding that a phase shift of $\pi$ exists between the echoes from the top and bottom surfaces due to Fresnel’s law [27], the echoes from the top and bottom surfaces of the tunnel can be precisely determined, and the diameters of the tunnels can be estimated accordingly, as shown in Fig. 5.
Fig. 5. Reflection time signals of (a) 0.8 mm, (b) 1 mm, (c) 2 mm, (d) 3 mm, (e) 4 mm, (f) 5 mm tunnels. The red vertical lines mark the positions of the echoes from the top and bottom surfaces of the tunnels. The signals are selected from the testing data.
Finally, deconvolution is used to deblur the tomographic scans of the sample. Because localizing a defect becomes more difficult as the defect gets smaller, we choose the cross-section images of the group of 0.8 mm tunnels to examine the effects of the deconvolution algorithms, as shown in Fig. 6. The 3D tomographic data is obtained by raster scanning the sample with the mechanical stage. Due to the vibration and the acceleration during the translation, the acquired time signals are noisier than those used in the discussion above, which are obtained by acquisitions from stationary positions.
Fig. 6. The tomographic images of the 0.8 mm tunnel that is closet to the top surface, as described in Fig. 1. (a) The original data, (b) deconvolved by FWDD, (c) deconvolved by ISTA, (d) deconvolved by LISTA. The bright line is the secondary reflection of the sample’s top surface, and the tunnels are labeled by the red squares.
The image generated from the unprocessed data displays obvious patterns caused by noise, and the tunnel further from the top surface is not fairly identifiable. The result of FWDD shows improvements in suppressing the noise patterns, but the contrast between the tunnels and the background is low. The ISTA notably eliminates the noise patterns, but the tunnels and the secondary reflection from the top surface of the sample are blurry due to the surrounding noise spikes. The LISTA completely removes the noise patterns and the visibility of the structures are greatly improved. To conclude, LISTA has obvious advantage over other deconvolution algorithms in analyzing the THz tomographic images, even when the input data has noise patterns that are not included in the training data.
4. Conclusions
In this paper, we described the application of LISTA in THz time signal processing and compared the deconvolution result with those of FWDD and ISTA. It is shown that LISTA is superior than others in resolving the impulse responses from low-SNR signals. Compared with ISTA, LISTA reduces computational time by over 100 folds, which is especially important in tomographic data processing. Because the HDPE material used in this study has ignorable dispersion effect that causes the broadening of the pulses, the discussions are made under the assumption that the waveform doesn’t change during the propagation. Further studies are required to account for the pulse spreading and the uneven distribution of refractive index through the medium.
Funding
Qingdao Quenda Terahertz 276 Technologies, Ltd; key research project of Qingdao.
Acknowledgments
Thank Wenping Li of Qingdao Quenda Terahertz Technologies for the experiment and data collection.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. D. M. Mittleman, S. Hunsche, L. Boivin, and M. C. Nuss, “T-ray tomography,” Opt. Lett. 22(12), 904–906 (1997). [CrossRef]
2. K. Su, Y. Shen, and J. A. Zeitler, “Terahertz sensor for non-contact thickness and quality measurement of automobile paints of varying complexity,” IEEE Trans. Terahertz Sci. Technol. 4(4), 432–439 (2014). [CrossRef]
3. J. Dong, X. Wu, A. Locquet, and D. S. Citrin, “Terahertz superresolution stratigraphic characterization of multilayered structures using sparse deconvolution,” IEEE Trans. Terahertz Sci. Technol. 7(3), 260–267 (2017). [CrossRef]
4. D. Zimdars, J. A. Valdmanis, J. S. White, G. Stuk, W. P. Winfree, and E. I. Madaras, “Time domain terahertz detection of flaws within space shuttle sprayed on foam insulation,” in Conference on Lasers and Electro-Optics/International Quantum Electronics Conference and Photonic Applications Systems Technologies (Optical Society of America, 2004), p. CThN4.
5. J. Dong, J. Bianca Jackson, M. Melis, D. Giovanacci, G. C. Walker, A. Locquet, J. W. Bowen, and D. S. Citrin, “Terahertz frequency-wavelet domain deconvolution for stratigraphic and subsurface investigation of art painting,” Opt. Express 24(23), 26972–26985 (2016). [CrossRef]
6. J. B. Jackson, J. Bowen, G. Walker, J. Labaune, G. Mourou, M. Menu, and K. Fukunaga, “A survey of terahertz applications in cultural heritage conservation science,” IEEE Trans. Terahertz Sci. Technol. 1(1), 220–231 (2011). [CrossRef]
7. S. Krimi, J. Klier, J. Jonuscheit, G. v. Freymann, R. Urbansky, and R. Beigang, “Highly accurate thickness measurement of multi-layered automotive paints using terahertz technology,” Appl. Phys. Lett. 109(2), 021105 (2016). [CrossRef]
8. Y. Chen, S. Huang, and E. Pickwell-MacPherson, “Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy,” Opt. Express 18(2), 1177–1190 (2010). [CrossRef]
9. H. Lin, Y. Dong, Y. Shen, and J. Axel Zeitler, “Quantifying pharmaceutical film coating with optical coherence tomography and terahertz pulsed imaging: An evaluation,” J. Pharm. Sci. 104(10), 3377–3385 (2015). [CrossRef]
10. H. Lin, Y. Dong, D. Markl, B. M. Williams, Y. Zheng, Y. Shen, and J. A. Zeitler, “Measurement of the intertablet coating uniformity of a pharmaceutical pan coating process with combined terahertz and optical coherence tomography in-line sensing,” J. Pharm. Sci. 106(4), 1075–1084 (2017). [CrossRef]
11. Z. Hua, X. Jingzhou, X. Xu, Y. Tao, R. Reightler, E. Madaras, and Z. Xi-Cheng, “Nondestructive defect identification with terahertz time-of-flight tomography,” IEEE Sens. J. 5(2), 203–208 (2005). [CrossRef]
12. T. Yasui, T. Yasuda, K.-i. Sawanaka, and T. Araki, “Terahertz paintmeter for noncontact monitoring of thickness and drying progress in paint film,” Appl. Opt. 44(32), 6849–6856 (2005). [CrossRef]
13. T. Yasuda, T. Iwata, T. Araki, and T. Yasui, “Improvement of minimum paint film thickness for thz paint meters by multiple-regression analysis,” Appl. Opt. 46(30), 7518–7526 (2007). [CrossRef]
14. J. Takayanagi, H. Jinno, S. Ichino, K. Suizu, M. Yamashita, T. Ouchi, S. Kasai, H. Ohtake, H. Uchida, N. Nishizawa, and K. Kawase, “High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser,” Opt. Express 17(9), 7533–7539 (2009). [CrossRef]
15. G. C. Walker, J. W. Bowen, J. Labaune, J. B. Jackson, S. Hadjiloucas, J. Roberts, G. Mourou, and M. Menu, “Terahertz deconvolution,” Opt. Express 20(25), 27230–27241 (2012). [CrossRef]
16. J. Dong, A. Locquet, and D. S. Citrin, “Terahertz quantitative nondestructive evaluation of failure modes in polymer-coated steel,” IEEE J. Sel. Top. Quantum Electron. 23(4), 1–7 (2017). [CrossRef]
17. P. P. Edward, M. Y. S. Stanley, B. Thierry, P. W. Vincent, and P.-M. Emma, “Terahertz pulsed imaging in vivo: measurements and processing methods,” J. Biomed. Opt. 16(10), 1–9 (2011). [CrossRef]
18. E. Carcreff, S. Bourguignon, J. Idier, and L. Simon, “High-resolution deconvolution applied to non destructive testing,” (2012).
19. D. M. Mittleman, “Twenty years of terahertz imaging [invited],” Opt. Express 26(8), 9417–9431 (2018). [CrossRef]
20. V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” IEEE Signal Processing Magazine 38(2), 18–44 (2021). [CrossRef]
21. K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” in ICML, (2010).
22. M. Sypek, J.-L. Coutaz, A. Kolodziejczyk, M. Makowski, and J. Suszek, “Aberrations of the large aperture attenuating thz lenses,” Proc. SPIE 8261, 826110 (2012). [CrossRef]
23. I. Selesnick, “Penalty and shrinkage functions for sparse signal processing,” (2013).
24. M. Elad, B. Matalon, J. Shtok, and M. Zibulevsky, “A wide-angle view at iterated shrinkage algorithms,” in Wavelets XII, vol. 6701D. V. D. Ville, V. K. Goyal, and M. Papadakis, eds., International Society for Optics and Photonics (SPIE, 2007), pp. 15–33.
25. C. Liu and M. Belkin, “Accelerating sgd with momentum for over-parameterized learning,” arXiv: Learning (2020).
26. X. Chen, J. Liu, Z. Wang, and W. Yin, “Theoretical linear convergence of unfolded ista and its practical weights and thresholds,” ArXiv abs/1808.10038 (2018).
27. K. Yamamoto, A. Masui, and H. Ishida, “Kramers–kronig analysis of infrared reflection spectra with perpendicular polarization,” Appl. Opt. 33(27), 6285–6293 (1994). [CrossRef]
