High-Resolution Long-Range THz Imaging for Tunable Continuous-Wave Systems

Imaging in the terahertz frequency range has attracted growing interests since the first image of a leaf more than 20 years ago, due to its countless applications in basic and applied research, medical imaging, and nondestructive testing. However, most terahertz imaging approaches rely on focusing optics which require knowledge about the imaging scene before the actual imaging takes place. Further, imaging is mostly restricted to short distances and high resolution is only achieved for systems with a high bandwidth. Here, we present a method that enables high-resolution imaging of small metallic and dielectric objects at distances up to 2 m based on a synthetic aperture. We derive a simple approximation for the resolution of partial circular synthetic apertures with limited bandwidth. The bandwidth limitation is encountered by replacing the measured signals with replica signals of high bandwidth and equal round-trip time so that the resolution is only limited by the carrier frequency and signal-to-noise ratio of the measurement system.


I. INTRODUCTION
Due to its unique properties, i.e., nonionizing radiation and high transmission through many dielectric materials, terahertz waves have been subject to intense research for different imaging applications during the last 20 years [1]. Applications demonstrated so far include security screening [2], [3], medical research [4] including cancer screening in particular [5]. Other areas of interest are non-destructive testing [6], hydration monitoring [7] and spatially resolved thickness measurements, e.g. of tablet coatings [8], [9], paint layer thickness [10] and in the field of art and heritage conservation [11], [12].
For these applications, the object under test (OUT) is typically located at the focus of a THz time-domain spectroscopy (TDS) system either in reflection [13] or in transmission [14] mode. These imaging setups provide a high lateral resolution if proper focusing optics are used. However, sharp focusing is achieved only for a very limited depth of The associate editor coordinating the review of this manuscript and approving it for publication was Zihuai Lin . field and the resolution is limited by the frequency-dependent spot size. Furthermore, the orientation of the OUT relative to the beam remains challenging for arbitrarily shaped objects. Finally, these techniques are typically used for short-range imaging where the distance to the OUT is fixed.
A possible solution to overcome these drawbacks are migration based methods for THz TDS which rely on lensless imaging [15] and [16]. In [17], we demonstrated an extended method for lensless imaging for broadband THz TDS systems which enabled the reconstruction of high-resolution images of samples with mm and sub-mm features without a priori information about the imaging scene. The sub-mm feature extraction was enabled by the high bandwidth (>1 THz) of the THz TDS system. However, the inherently low signal-tonoise ratio of a typical THz TDS system with a divergent beam, restricts the imaging to short-range (around 0.3 m range). Further, the imaging range, i.e. the spatial window where a reflection from an OUT can be detected, is limited for THz TDS systems to several hundreds of millimeters due to the commonly used optical sampling. Hence, a different VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ approach is required for long range imaging of unknown objects.
In this article, we present a promising imaging approach for long-range terahertz imaging with resolution in the submm range based on a vector network analyzer (VNA) with two extenders for the THz range. Due to the high signalto-noise ratio (SNR), a VNA system can detect reflections from objects in the mm range even at large distances. The high frequency resolution yields to a large unambiguous measurement distance. However, VNA systems operate in contrast to THz TDS systems with a limited relative fractional bandwidth which limits the resolution. Further, the migration based imaging methods for THz TDS are so far developed for the time-domain while VNA systems operate in the frequency domain. Therefore, we present a frequency domain formulation for lensless imaging suitable for VNA systems. Further, we derive an approximation for the imaging resolution for restricted circular synthetic aperture with limited bandwidth. Based from this approximation and the limited fractional bandwidth of VNA systems, we propose an algorithm to overcome this limitation. The algorithm is based on the image reconstruction with a broadband replica signal with equal round-trip time replacing the narrowband measured signal. In this case, the resolution is only limited by the carrier frequency and the SNR of the system. We demonstrate that a clear image reconstruction of a 5 mm hex wrench and a wooden lead pencil is possible at a distance of 2 m with a bandwidth of just 20 GHz.

II. LENSLESS IMAGE GENERATION
For the image generation, we employ the migration-based radar imaging method adapted for lensless terahertz systems as presented in [17] extended for tunable continuous-wave measurement systems.
As a test bed for the imaging method, the scenario shown in Fig. 1 is considered. Here, we have two single-point targets at position P 1 and P 2 . An emitter-detector antennas pair is used for reflection measurements of the single-point target and after each measurement the antennas are moved along a circular track resulting in a circular synthetic antenna aperture. Using a continuous wave system, we can transmit a single frequency wave with an emitter positioned at E 0 . The target reflects the signal, which after a round-trip time t i,n is captured by the detector at position D 0 . Thus, the reflection coefficient of the target for the given frequency f 0 is measured. If we now perform a frequency sweep starting from f 0 to f 0 + B, we can measure the response of the target to a signal with a rectangular-shaped spectrum with a frequency bandwidth B and a center frequency f c = f 0 + B/2. The detected signal at the receiver is then given as where D n and E n denote the positions of the detector and emitter antennas for the n-th measurement. P i is the position of the i-th point target and r(P i , D n , E n ) is the reflection factor describing the terahertz wave reflected from the target P i onto the detector. r(P i , D n , E n ) is a function of the sample shape, and angle of incidence of the terahertz wave on the sample, assuming a frequency-independent reflection factor for the single-point target.
To generate the image of the target, the frequency-domain signal S R (f , D n , E n ) is transformed into the time-domain signal s R (t, D n , E n ) using the Inverse Fourier Transformation. The detected signal in the time-domain after s R (t, D n , E n ) can be expressed as the where s t (t) denotes a time-domain signal with a rectangularshaped spectrum. Assuming a constant wave propagation velocity c, the round-trip time resulting from the i-th point target and n-th measurement is where c is the speed of light. Subsequent traces are then considered for different antenna positions by moving the emitter-detector antenna pair generating a synthetic antenna aperture. Here, we generate a circular synthetic antenna aperture by rotating the antennas with an angle ϕ around the targets as shown with the dotted line in Fig. 1. The relative position of the antennas to each other is unchanged. Therefore multiple reflections are collected from the same sample. Using the information from the different traces and the knowledge about the respective positions of emitter and detector antennas for each measurement, we can generate a reflectivity map I (x m , y m ) of the sample area, essentially a sum of the back-projected terahertz pulses which is where t m,n = |P m −E n |+|P m −D n | c is the round-trip time of the pulses to and from an assumed single-point target at P m and e (j2πf 0 t) is a phase correction term accounting for the phase rotation during the wave propagation.
Eq. (4) yields a 2-dimensional matrix where every matrix element corresponds to the reflectivity coefficient of a possible single point scatterer at (x m , y m ). Pixels with high intensities show the contour of the object.

III. RESOLUTION EVALUATION
To evaluate the resolution of the THz lensless circular synthetic antenna aperture, we offer a mathematical model for the image resolution of a single point target. Based on the model, numerical solutions for the −3 dB resolution in x-direction (azimuth) and y-direction for partial synthetic aperture and different frequency bandwidths are presented. Furthermore, a simplified set of approximation equations for calculation of the −3 dB resolution for circular synthetic aperture are presented based on the numerical solution. Finally, a simulation is presented where the effects of the interference of two single-point targets along the y-axis (cf. Fig. 1) are analyzed.
First, we present a mathematical model for the image resolution based on the point spread function (PSF) of a single-point target. We consider the effects of the carrier frequency, the system frequency bandwidth and the integration angle θ of the circular synthetic aperture. The angle θ defines the range of the circular track for which a measurement is performed and a reflection is detected. Assuming a single-point target at the center of the coordinate system with frequency-independent reflection factor (r(P i , D n , E n ) = 1), the reconstructed THz image of the given target is I SP (x m , y m ). The Fourier Transform of the image transformed in polar coordinates is with ρ min = 2πf 0 /c, ρ max = 2π(f 0 + B)/c, ϑ min = −θ/2, and ϑ max = θ/2, where θ is the integration angle of the circular synthetic aperture. The point spread function in polar coordinates of the single-point target is then To evaluate the azimuth x and range y −3 dB resolution, the values for r for which I SP (r, ϕ) = −3 dB for ϕ = 0 • and ϕ = 90 • , respectively, have to be calculated. For a complete circular scan (θ = 360 • ) and a narrow frequency bandwidth system (f c B), a solution for (6) is presented in [22] and [23]. As a complete 360 • scan is performed, the function for the PSF in not dependent on ϕ anymore. Hence, the azimuth and range resolution are equal ( x = y) and the −3 dB resolution of the imaging system is where λ c is the wavelength of the carrier frequency. Hence, the imaging resolution of the system is only dependent on the center frequency. Moreover, using a circular synthetic aperture, imaging with sub-wavelength resolution is possible.
In [24] and [25], a solution for (6) is presented for broadband circular synthetic aperture systems with θ = 360 • as the difference of two weighted Bessel functions of 1 st order. Similar to the solution proposed in [22] and [23], here the azimuth and range resolution are equal and the image resolution is only a function of the carrier frequency. The larger bandwidth affects only the amplitude of the side lobes of the reconstructed single-point target.
For practical applications θ < 360 • . Hence, we consider a partial circular track for the reconstruction of the point target, i.e. θ < 360 • . In [26], a solution for (6) is presented for a partial circular track. However, only a narrow-band system is considered. An approach for broadband systems is demonstrated in [27] and a solution for (6) is presented for the calculation of the PSF according to where and I s (r, ϕ, θ ) The azimuth resolution (in x-direction) and the range resolution (in y-direction) can be calculated for r with ϕ = 90 • and ϕ = 0 • . Using equations (8) -(10), we evaluated the PSF for a THz tunable continuous-wave system with the following parameters: 5 • < θ/2 < 90 • , 10 GHz < B < 110 GHz, and f c = 275 GHz. The −3 dB values corresponding to the azimuth resolution x and range resolution y for each evaluated PSF are extracted and shown in Fig. 2 and Fig. 3, respectively. As shown in Fig. 2, the azimuth resolution is mainly dependent on the integration angle θ and the center frequency of the system. However, as demonstrated in Fig.3 the range resolution is dependent on θ, the center frequency and the frequency bandwidth of the system.
Based on the results from Fig. 2 and Fig. 3, we propose a set of equations for the calculation of the range and azimuth VOLUME 8, 2020  −3 dB resolution for broadband systems. As the azimuth resolution ( x) is mainly dependent on the integration angle θ and the center frequency, a simplified equation can be derived as The proposed equation (11) is similar to equation (6) for a complete circular scan. In the first part of the equation for 0 • < θ < 180 • , the azimuth resolution is degraded by the partial circular synthetic aperture. However, for larger synthetic aperture θ > 180 • , the azimuth resolution is approximately the resolution of the imaging system when the complete circular synthetic aperture is considered. Similarly, we propose an equation for the calculation of the −3 dB range resolution for circular synthetic aperture. The new proposed equation is split into two parts depending on the integration angle θ. The proposed approximation function for the −3 dB range resolution is According to the proposed range resolution equation, for 0 • < θ < 120 • the range resolution is dependent on the bandwidth and improves exponentially with the integration angle. Then, for large antenna aperture (θ > 120 • ) the range resolution is dependent only on the carrier frequency. Next, we evaluate the error of the new proposed set of equations. The results from new proposed equations for the azimuth and range −3 dB resolution, (11) and (12), are compared with the results from Fig. 2 and Fig. 3 derived from equation (8). Furthermore, using the same parameters as for Fig. 2 and Fig. 3, the imaging method from the previous section is employed to generate an image using the simulated data for a single-point target at the center of the coordinate system. For the simulation and further evaluation of the imaging algorithm, the antennas remain at a fixed position and the targets are rotated with respect to the center of rotation (0,0) (cf. Fig. 1) to generate a 2-dimensional image of the target. The −3 dB beamwidth of the imaged single-point target in azimuth and range is then extracted and compared with the results from the new proposed equations. A mean value for the error of the new proposed equation for the −3 dB azimuth resolution is calculated to be µ e x = 0.135 · λ c and the mean value for the error of the proposed equation for the −3 dB range resolution is µ e y = 0.085 · λ c . Hence, the two new proposed equation provide a good approximation for the resolution limits of the imaging setup.
To further analyze the resolution limits of the imaging method, we perform simulations using now two single-point targets. Thus, we can evaluate the correlation between the measurement parameters such as frequency bandwidth and synthetic aperture, and the resulting image resolution in the case of interference of two single-point targets. Table 1 lists the parameters for the different simulations. θ defines indirectly the range of angles for which a reflection from the targets was obtained. For the single-point target at P 1 , a reflection is obtained for ϕ 0 − θ/2 < ϕ < ϕ 0 + θ/2 and for the single-point target at With the given parameters from Table 1, we simulate signals from a rectangular shaped spectrum with a frequency bandwidth B and a center frequency of f c . In Fig. 4, examples of three simulated time-domain signals with 20, 60, and 100 GHz bandwidth and a carrier frequency of 275 GHz are shown. The envelope of the signals are depicted in red and the simulated electric fields are depicted in blue. As the fractional bandwidth B fr = B/f c is smaller than 1, multiple oscillations of the electric field are observed under the envelope of the signal.  In the following, the imaging results for the simulation are discussed. Fig. 5 shows the results from the simulation when we consider the worst case for the system bandwidth B = 20 GHz. Here, the frequency bandwidth is constant and the different values for θ are considered. From Fig. 5, we can state that for constant frequency bandwidth the image resolution increases with increasing θ, which corresponds to the synthetic antenna aperture. In the worst-case scenario (θ = 20 • ) the two targets can not be resolved. They merge into a single object centered at the axis origin. By increasing θ to 40 • the resolution in y does not change and only a slight improvement of the resolution in x is observed. The results for θ = 60 • show an improvement for the resolution in y and the two single-point targets are to a certain degree resolved. However, strong interference of the two targets leads to image artifacts and the targets are stretched in the y dimension. Similar results are observed for θ = 90 • , but the point-targets are no longer stretched in the y-direction and the image artifacts are attenuated. Finally, with a frequency bandwidth of just 20 GHz the two targets can be correctly resolved for θ > 180 • . However, for a real-world imaging scenario, the synthetic aperture may be limited to smaller values.
Next, we consider the worst-case scenario for the synthetic antenna aperture of θ = 20 • for different frequency bandwidths. In this manner we can analyze the correlation between image resolution and system bandwidth for a constant synthetic antenna aperture. The results are shown in Fig. 6. The first results for B = 20 GHz is the same as from Fig. 5.
Here, no differentiation between the two targets is possible. By increasing the frequency bandwidth to B = 40 GHz the two targets are still not correctly resolved. With a system bandwidth of 60 GHz or fractional bandwidth B fr ≈ 0.22 the targets are correctly resolved. However, the position of the targets is not correctly reconstructed as the target images in the figure are at (0, −2.95) mm and (0, 2.95) mm for P 1 and P 2 , respectively. By further increasing the system bandwidth, the targets are now correctly resolved and are located at the correct coordinates. Finally, the results from the simulation that have the strongest relation to real-world imaging applications, θ = 20 • , 40 • , 60 • and all different system bandwidths are summarized in Fig. 7. Here, the imaging resolution in y is represented by the profiles of the main lobe of the two reconstructed targets. For this purpose slices of the images along the y-axis for x = 0 are shown to analyze the beam profile of the reconstructed targets. For θ = 20 • the profiles of the beams of the reconstructed targets from Fig. 6 are shown. It can be seen that no differentiation between the targets is possible for B < 60 GHz. Furthermore, we can see the shift in the y-direction of the two reconstructed targets for the case of B = 60 GHz. A reason for this error is the interference between the main lobes of the targets and the side lobes. In the case of B = 80 GHz and B = 100 GHz we can see the correct reconstruction of the two point targets for θ = 20 • and the constructive interference of the side lobes resulting in a strong peak at (0, 0) mm. Similar results are observed in the case of θ = 40 • and θ = 60 • . However, an attenuation of the side lobes and a narrower main lobe of the targets can be observed for the increasing synthetic antenna aperture. To further analyze the beam profiles of the reconstructed targets, the −3 dB beamwidths of the interfering main lobes are given in Table 2. From the table, we can see that for B = 20 GHz the −3 dB beamwidth changes just slightly and no real change is present for the cases of B = 60 GHz, B = 80 GHz and B = 100 GHz. The simulation demonstrates that for a circular synthetic aperture, the image resolution is a function not only of the antenna aperture, but also of the frequency bandwidth B. Hence, in the case of restricted constant antenna aperture, the image resolution is a function of the system bandwidth. In the next section, we propose an imaging algorithm for a continuous-wave system for a restricted antenna aperture and frequency bandwidth.

IV. IMPROVED IMAGING METHOD FOR CONTINUOUS-WAVE SYSTEMS
THz frequency systems, such as frequency modulated continuous wave (FMCW) radars or vector network analyzers (VNA), provide high-resolution frequency-domain signals with a high dynamic range at large distances, but a small fractional frequency bandwidth. The VNA systems provide typically a frequency bandwidth from 60 GHz to up to 250 GHz at the lower THz frequency spectrum (from 100 GHz up to 750 GHz). The FMCW systems have even smaller frequency bandwidth from 60 GHz [19] to 136 GHz [20], which as we show in the previous section limits the imaging resolution of continuous-wave systems for a restricted synthetic antenna aperture.
For the generation of the THz images we employ a Rohde&Schwarz VNA ZVA67 system with two millimeterwave extenders (R&S R ZC330). The system can sweep the frequency from 220 GHz to 330 GHz. The emitter extender is connected to Port 1 and the detector extender is connected to Port 2 in a bi-static configuration. The signal is then transmitted over broadband horn antennas and the reflection from a target is measured as the S 2,1 parameter. To increase the system imaging resolution, we propose a new imaging algorithm for continuous-wave systems with a fractional bandwidth B fr < 1. The algorithm generates a replica signal of the reflected signal using a waveform with a fractional bandwidth B fr > 1.
In Fig. 8, a block diagram of the algorithm is shown. The measured signal S R (f , D n , E n ) is the input of the algorithm. In the first step of the algorithm, a calibration procedure is proposed to compensate for the antenna pattern, antenna crosstalk, errors at the system calibration, and radar clutter in the received THz signal. The signal is calibrated using two reference measurements. To compensate for the antenna patterns and the system transfer function a reference measurement with a large metal plate as a reflector is performed. The second calibration measurement is performed to calculate the frequency-dependent signal-to-noise ratio (SNR) of the system. To calculate the SNR of the system, a measurement, where no target is present is performed for the complete frequency range. This measurement is then considered as the noise component for the calculation of the SNR. The calibration of the measured signal is performed according to where H Ref (f ) is a reference measurement with a large metal plate as a reflector and SNR(f ) is the measured SNR of the system. In Fig. 9 an example for a measurement with 110 GHz bandwidth at 275 GHz center frequency is shown for the raw data and the calibrated signal. As seen from the raw signal (blue curve), there are spikes for specific frequencies In the next two steps of the algorithm, the time-domain signal is calculated and the resolution of the signal is virtually increased by a zero-padding in the frequency domain. We employ the Hermitian signal processing technique for VNA systems [21]. Using this technique the signal is zero-padded from the lowest frequency, here 275 GHz, to direct current (DC). Next, we take the conjugate complex of the signal and reflect it to the negative frequencies. By using the Inverse Fast Fourier Transform algorithm, the symmetric to DC double-sided spectrum signal is transformed in the time-domain resulting in a real time-domain signal and improving the accuracy of the VNA system.
In the following step, the algorithm generates the envelope of the time-domain signal. Using the envelope signal, the reflections from the objects are detected to decrease the noise and clutter.
As a detector of the reflected signal, we employ a cell-averaging constant false alarm rate (CA-CFAR) detector and a reflection is detected according to the following decision rule where s E (t k , D n , E n ) is the value of the envelope of s R,F (t, D n , E n ) at time t k , H 1 is the hypothesis that a target is located after time t k , H 0 is the noise-only hypothesis and δ CA (t k ) is the cell-averaged local threshold adaptively obtained for each s E (t, D n , E n ). Assuming white Gaussian noise, the threshold is derived according to where P FA is the value for the constant false alarm rate, T defines a time-interval around the considered t k , f a is the sampling frequency, and µ t k is the expectation value of the signal s E (t, D n , E n ) for a time period T around t k , or Here, a constant false alarm rate of P FA = 0.1 is chosen. Once the reflections are detected, a new replica signal is generated based on the detected round-trip times extracted from the signal envelope. The replicated signal is reconstructed using the following equation where t m is the m-th detected round-trip time for the m-th detected reflection, M n is the number of detected reflections for the n-th measurement, A m is the amplitude of the signal s R,F (t, D n , E n ) for t = t m and s W (t) is a new waveform we use to generate the replica signal s RS (t, D n , E n ). For this work, we investigate three different broadband waveform signals s W (t): 1) a half Gaussian shape waveform symmetrical to DC with B = 2f c resulting in a full width at half maximum (FWHM) of the time-domain signal of t FWHM = 1/f c , 2) a Dirac impulse response, 3) and an ideal single-cycle Gaussian pulse with a carrier frequency f c , frequency bandwidth B = 2f c resulting in a FWHM of the envelope of the time-domain signal of t FWHM = 1/f c . In Fig. 10, the three different waveforms are shown compared to a simulated signal with f c = 275 GHz and 20 GHz frequency bandwidth. As seen from the figure compared with the raw signal, the waveforms of the replicated signals are much shorter due to the larger bandwidth. To analyze the resolution capabilities of the proposed pre-processing algorithm, first the simulation from the previous section is repeated for the three different waveforms. However, only the parameters for the worst conditions are considered: θ = 20 • and B = 20 GHz. In Fig. 11, the reconstructed images using the proposed pre-processing method are presented. Here, similarly to Fig. 7, the profiles of the main lobes of the two reconstructed targets using the pre-processing algorithm are shown. For all three waveforms, the two targets are reconstructed correctly. VOLUME 8, 2020 FIGURE 11. A single slice of the imaging results using the reconstructed replica signals for x = 0.
The −3 dB beamwidths of the different waveforms are 370 µm, 100 µm, 175 µm, for the half Gaussian-shaped waveform, for the Dirac waveform and the single-cycle Gaussian waveform, respectively. All considered waveforms show smaller −3 dB beamwidth values compared to the results from Table 2 corresponding into an improvement of the system resolution by a factor of 22, 82 and 47, respectively. In Fig. 11 the difference of the dynamic range (the difference between the smallest and largest image values) of the reconstructed images with the three waveforms is demonstrated. The value for the single-cycle Gaussian pulse shows the best results with a peak dynamic range of 70 dB compared to the 39 dB of the Dirac shaped pulse and the 24 dB of the half Gaussian-shaped pulse.
The simulation demonstrates that for a restricted circular synthetic aperture, the image resolution can be improved by employing the image reconstruction algorithm proposed here.

V. EXPERIMENTAL RESULTS
So far, we have only considered simulated signals. Now, we apply the proposed algorithm using real measured data. An object under test is mounted on a rotation unit where the center of rotation of the rotation unit is 1.35 m away from the two VNA extenders as depicted in Fig. 1. The distance between the emitter and detector extenders is 0.25 m. The test object is rotated from 0 • to 360 • with 1 • resolution. After each rotation step a measurement is performed allowing for a complete 360 • 2D circular scan of the test object. The VNA system performs the measurements using a complete frequency sweep from 220 GHz to 330 GHz with 55 MHz resolution resulting in 2001 sampled frequency points of the channel. This raw data from the VNA system is then used for the further generation of the 2D image of the test object.
As an imaging target, a 5 mm hex-wrench is employed resulting in a 2-dimensional object with 6 weak reflectors (the 6 edges of the hex-wrench) and 6 strong reflectors (the 6 sides of the hex-wrench).
To evaluate the performance of the three proposed waveforms as a replica signal, the data is first evaluated using the proposed algorithm without the steps ''Simulated replica'' and ''Reflection extraction (CFAR)''. The algorithms from the previous section is then used to generate images based on the different waveforms as a basis of the replica signal. Three different frequency bandwidths are considered, B = 20 GHz, B = 40 GHz, B = 80 GHz with a center frequency f c = 275 GHz. In Fig. 12 the results from the imaging algorithm are presented.
First let us consider the images a), b) and c) from Fig. 12 which are reconstructed from the original measurement data. Using just 20 GHz of the frequency bandwidth makes the reconstruction of the 5 mm hex-wrench impossible. By backprojecting the reflections from the different details of the object, the interference between these reflections result in a distorted image without any recognizable shape. Similar is the result where we considered 40 GHz frequency bandwidth. By increasing the frequency bandwidth of the system to 80 GHz the general shape of the object becomes clearer, but with strong image artifacts. The shape of the object is now visible. However, the sidelobes of the object's contour is just ≈ −1.5 dB smaller than the main lobe and strong artifacts arise during the reconstruction. Now, for the image reconstruction the three new replica waveforms for the generation of a replica signal are considered. First, in Fig. 12.d) -f) the results using a half Gaussian shape pulse for the replica signal with t FWHM = 1/f c are presented. Here, all edges of the hex-wrench are reconstructed for all different bandwidths. However, due to the high pulse width the target is blurred and no clear contour can be recognized. Although the image artifacts are lower compared to the raw data, the intensity of the artifacts is still high leading to a distorted image.
Next, we consider the Dirac shaped pulse as the replica signal. The results from the imaging algorithms are shown in Fig. 12.g) -i). The object is reconstructed correctly for all different bandwidths considered. The images differ only in the intensity of the image artifacts, where the artifacts intensities decrease with higher bandwidth. The reason for this is the pulse shape used for the replica signal. Using a Dirac shaped pulse makes the image reconstruction more susceptible to errors in the reconstruction (artifacts) in the case of jitter or false reflection detection. If the position of the reflection is not correctly detected due to the broader pulse width of the original data, the position of the time-shifted Dirac pulse is also incorrect leading to errors in the image reconstruction.
Now we present the results of the single-cycle Gaussian pulse for the replica signal in Fig. 12.j) -l). Similarly to the results with the Dirac shaped pulse, the object is correctly reconstructed for every considered frequency bandwidth. However, in contrast to the images obtained with the Dirac shaped pulse, the image artifacts are reduced due to the higher dynamic range as demonstrated by the simulated data. Hence, we can conclude from these measurement results that the single-cycle Gaussian pulse outperforms the Gaussian-shaped and Dirac shaped pulses for image reconstruction. To evaluate the resolution of the imaging method using a single-cycle Gaussian pulse as a replica signal we measured the −3 dB width of the side of the hex-wrench for the three frequency bandwidths and calculated the mean values. The calculated mean values are 153 µm, 148 µm, and 145 µm using 20 GHz, 40 GHz and 80 GHz frequency bandwidth, respectively. Hence, applying replica waveforms for the image reconstruction greatly enhances the quality and the resolution of the reconstructed image.
Finally, to demonstrate the high-resolution imaging capabilities at larger distances, the same object is measured at 1.5 m and 2 m. Furthermore, to show that the algorithm can reconstruct the image of not just metal, but also of dielectric materials a simple wooden hexagon pencil with a diameter of ≈7 mm is imaged at 2 m distance. Only the single-cycle Gaussian pulse is considered for the generation of the images. The results from the image reconstruction are shown in Fig. 13. The images of the 5 mm hex-wrench at 1.5 m and 2 m show a correct reconstruction of the object shape and size. The images differ just slightly from the image at 1.35 m, thus, proving that the proposed algorithm can correctly image small metal targets at larger distances. The shape and size of the wooden pencil are also correctly reconstructed at 2 m distance showing that the algorithm is not only limited to high-reflective materials. Furthermore, as reflections from the core of the pencil are also detected, high-intensity image components inside of the pencil are visible in the image, but the shape and size of the core can not be correctly reconstructed as the refractive index of the wooden pencil is not considered for the image generation. For the correct reconstruction of layered material first the material of the different layers must be correctly characterized.

VI. CONCLUSION
In this article, we presented a new imaging approach for long-range high-resolution THz imaging for a tunable continuous wave systems with limited fractional frequency bandwidth. An approximation for the image resolution for a given bandwidth and for a restricted circular synthetic aperture is retrieved. We demonstrated by a simulation that the image resolution is not only a function of the synthetic antenna aperture and center frequency but also of the fractional bandwidth. This limits the imaging resolution for systems with inherent low bandwidth. To overcome this challenge, we proposed an algorithm which is based on a reconstruction of a replica signal of the reflected signal using a broadband single-cycle Gaussian waveform. By employing the proposed method, we showed that the resolution is in this case only limited by the carrier frequency of the system as long as the position of the peak can be correctly identified. This leads to high-resolution imaging with a −3 dB beamwidth of 170 µm for a bandwidth of 20 GHz with a center frequency of 275 GHz. The imaging technique relies on a divergent THz beam and experimental results showed the capability of the method to resolve objects in the sub-mmregion for long-range imaging without any a priori knowledge about the imaging scenario. The results showed the correct image reconstruction of metallic objects with a diameter of just 5 mm and an edge length of 2.9 mm. Furthermore, a measurement was performed using a wooden pencil with a 7 mm diameter at a range of 2 m. Using the proposed algorithm, the shape and size of the object was correctly reconstructed, while also revealing the internal structure. Further investigation of the imaging method will be performed for object exhibiting concave and convex surfaces before extending the method for 3-dimensional image reconstruction of multilayered objects.
AMAN BATRA (Student Member, IEEE) received the B.Tech. degree in electronics and communication engineering from Maharshi Dayanand University, India, in 2014, and the M.Sc. degree in embedded systems engineering from the University of Duisburg-Essen (UDE), Germany, in 2017, where he is currently pursuing the Ph.D. degree with the Institute of Digital Signal Processing.
In 2012, he underwent an eight weeks internship at the Centre for Railway Information Systems, India, and from January 2014 to May 2014, he did an industrial internship with the Department of Signaling and Telecommunication, Indian Railways, India. During his master's studies tenure, he has worked as a Student Research Assistant at UDE on implementation of hardware and software design on System on Chip and building embedded Linux drivers for software defined radio. He is also working for the project Collaborative Research Center / Transregio 196 Mobile Material Characterization and Localization by Electromagnetic Sensing (MARIE). In this job, he is also working on Terahertz and Millimeter Wave radar imaging with real-time operations. His research domain includes radar and massive MIMO signal processing, SAR imaging, high-performance computing and embedded programming. Moreover, he is also a part of organizing committee of the IEEE IWMTS series and has also coordinated various seminars and workshop under the agenda of MAIRE. In December 2019, he has also visited the Blekinge Institute of Technology, Sweden, as a Visiting Scholar, for two weeks.