A Non-Linear Correction Method for Terahertz LFMCW Radar

The nonlinearity in terahertz (THz) linear frequency modulated continuous wave (LFMCW) radar usually blurs the range profile and decreases the signal to noise ratio, hampering applications where high range-resolution is particularly emphasized. A software correction method, which comprises of transmitted nonlinearity estimation and nonlinear phase compensation for the beat signal, is proposed in this paper to drastically reduce the nonlinearity in THz LFMCW radar. Besides the commonly considered nonlinearity caused by voltage-controlled oscillator (VCO), the nonlinearity from other broadband hardware devices have also been included in our modified correction model, which gives the advantage of preciser compensation. Moreover, utilizing the phase gradient autofocus (PGA) method to estimate the transmitted nonlinear term and the residual video phase (RVP) removal method to remove the range dependency of the received nonlinearity, our method can uniformly compensate the nonlinearity in the whole range profile. In addition, no presupposed parametric model for the nonlinearity waveform is needed, which further strengthens the effectiveness of the proposed method in practical use. Both the simulated data and the real tested data, acquired by a 190 GHz radar with 60 GHz bandwidth, has been used to demonstrate the validity and the effectiveness of the method.

Apart from the hardware correction methods [12], [17], [22], [23], which solve the nonlinearity problem in the signalgenerating stage, many software solutions have also been proposed. They usually operate directly on the beat signal and offer the advantage of requiring no extra circuit, which helps to realize miniaturization, low power assumption, and low cost in industrial applications. There are two categories of software methods: direct beat signal correction methods [10], [12], [18] and methods based on the signal model [19]- [21], [24]. With the help of the reference target or some prior information, the direct beat signal correction methods [10], [12], [18] focus the echoes of the targets by uniformly eliminating a nonlinear phase term, regardless of the signal model, targets range, and the system's architecture. Therefore, the effectiveness of the methods is highly range dependent [25], [26]. On the other hand, with the consideration of the signal model, the nonlinearity in the beat signal is divided into two parts: the nonlinearity in the transmitted signal and in the received signal. Based on different kinds of simplifications, several methods [19]- [21], [24] have been proposed to solve the nonlinearity problem for the whole range profile. However, the effectiveness of this kind of methods is limited by their simplified signal models. In addition, these model-based methods only considered the dominant nonlinearity caused by voltage-controlled oscillator (VCO). The nonlinearity of other hardware devices in the system have not been included so far, which leads to range profile degradation especially for large bandwidth THz systems. Based on the signal model we proposed, these two kinds of methods will still be discussed in details in the section II. B.
In this paper, in order to solve the nonlinearity problem, we propose a novel and practical software correction method, which integrates the advantages of these two kinds of methods. Specifically, with the extra considerations of the nonlinearity caused by the THz broadband hardware devices, we firstly propose a modified signal model. Based on this model, we use the phase gradient autofocus (PGA) method to estimate the transmitted nonlinearity and then use the residual video phase removal (RVP) method to compensate the received nonlinearity in the whole range profile by making the nonlinearity effect range independent. Both the simulated data and the real tested data, which are acquired by a 60 GHz bandwidth THz radar, have been used to demonstrate the validity and the effectiveness of the method in qualitative and quantitative way.
Overall, the contributions of this paper are as follows: i) not only the dominant nonlinearity caused by VCO, but also the nonlinearity of the broadband hardware devices, such as the frequency multiplier, harmonic mixer, and antennas et al, have been included in the proposed signal model. This approach offers the advantage of precise compensation, especially for the wideband systems like THz LFMCW radars.
ii) for this modified signal model, we propose a new software correction method. Firstly, with the help of two reference targets, PGA method is used to estimate the trans- mitted nonlinearity. Then, by using the RVP method, the range dependent received nonlinearity in the signal can be uniformly eliminated. Therefore, the nonlinearity of undertest targets in the whole range profile can be corrected all at once.
iii) no pre-assumption of the specific parametric model for nonlinearity waveform is needed for our correction method, which further strengthens the effectiveness of the proposed method in practical use.
The rest of this paper is organized as follows. Section II firstly introduces the signal models. Then, based on this model to correct the nonlinearity, the PGA method is proposed to estimate the transmitted nonlinearity and the RVP method is proposed to correct the residual received nonlinearity. In section III, both the simulated data and the real test data are used to demonstrate the validity of the method. Finally, section IV draws the conclusion.

A. MODEFIED NONLINEAR SIGNAL MODEL
The diagram of our compact homodyne THz LFMCW sensor is shown in Figure 2. Driven by a linear ramp voltage, the VCO produces LFMCW signal in the microwave frequency band. The output signal of VCO is divided into two parts, one is fed into a frequency multiplier to reach the working THz frequency and the other is fed into a harmonic mixer for local oscillator (LO). The back-reflected signal from the objects return into the receiving antenna and is guided to the harmonic mixer. The received radio frequency (RF) signal is de-ramped by the LO, producing the homodyne frequency (zero intermediate frequency, Zero-IF) signal.
When the FMCW of RF and LO sweep linearly in time, the frequency of Zero-IF is proportional to the time delay (timeof-flight). Thus, the distances between the sensor and the targets are measured by the frequency components of Zero-IF signal.
An ideal LFMCW signal can be expressed as: where f 0 is the carrier frequency, t is the time varying within the pulse repetition interval T , and K is the frequency sweeping rate which equals to the ratio between the bandwidth B and the interval T . However, the linearity of the signal is often degraded by the non-linear characteristic of VCO. For example, the varactor-tuned method is often served as the voltage-variable oscillator to control the output frequency because of its fast tuning speed. But according to the varactor model, the relationship between the control voltage and its output frequency is not linear [27], and is also sensitive to load impedance variation over the bandwidth and even temperature [17]. Thus, the nonlinearity of VCO is degraded by the increase of sweeping bandwidth. Besides, the nonlinearity of THz broadband devices could not be neglected [25], [26], [28]- [30] because of their drastic group delay jitter [31]- [34], especially for electron devices like frequency multiplier and harmonic mixer in our system.
With the consideration of the nonlinearity in VCO, the harmonic signal of LO in the working THz band (the transmission path is illustrated with blue dashed line in Figure 2) can be expressed as the combination of ideal chirp signal with a nonlinear phase term ε(t): For the simplicity of the problem, the amplitude envelope of the signal is excluded in this paper. Then, after the propagation delay, the received RF signal (the transmission path is illustrated with red dashed line in Figure 2) can be expressed as: where τ = 2R/c is the time delay caused by the target and R is the distance between the target and the sensor. It worth mentioning that with the additional consideration of the influences caused by the cascaded broadband THz devices, like the frequency multiplier, antennas, and harmonic mixer in our system, we particularly include this extra nonlinearity in the nonlinear phase term ξ (t) in RF signal. Thus, the received nonlinear phase term ξ (t) (nonlinear phase term in RF) is different from the transmitted nonlinear phase term ε(t) (nonlinear phase term in harmonic signal of LO), which makes the problem more complicated. Then in the de-ramping procedure, the harmonic LO signal S LO−har (t) and the RF signal S RF (t) are mixed to get the Zero-IF signal: where the rectangular time window rect(t) guarantees the effectiveness of the beat signal and eliminates the errors caused by discontinuity of the modulation. According to (4), with the presence of the nonlinear phase term ε(t) − ξ (t − τ ), the beat signal is no longer a monochromatic signal. The nonlinearity spreads the target well-focused spectrum response to different frequencies, resulting in range resolution degradation, loss of SNR, side lobes raising, and main lobe dissymmetry, as illustrated in Figure 1. Thus, the nonlinearity compensation problem equals to the problem of eliminating the nonlinear phase term ε(t) − ξ (t − τ ) from the Zero-IF signal.

B. ALGORITHM DESCRIPTION
Many software methods have been proposed to eliminate the nonlinear phase term ε(t) − ξ (t − τ ) from the IF signal. One type of methods [10], [12] regard the nonlinear phase term as a fixed whole part ε(t) − ξ (t − τ ref ) (do not change with the time delay) and uniformly remove it from the beat signal. They are very effective in the limited range because there is no model assumption, but the correction effect decays sharply with the offset distance from the reference [25], [26]. The other type of methods only consider the dominant nonlinearity in VCO and simplify the problem with the hypothesis that the received RF signal is only a time delayed version of the harmonic signal of LO: (4)). Then, with the help of the RVP removal method [19], [20], the near range differential approximation ε(t) − ε(t − τ ) ≈ ε (t)τ [24], or with the further pre-assumed parametric model for term ε(t), like quartic polynomial in [21], the nonlinearity in the whole range profile can be compensated in different ways by making the nonlinearity effect range independent. Finally, the nonlinear phase term ε(t) − ε(t − τ ) from the IF signal can be eliminated.
With the extra considerations of the influence caused by the broadband THz devices, as illustrated in (4), we proposed a new correction method which integrate the advantages of these two kinds of methods. By the help of the reference targets, we firstly use the PGA method to estimate the transmitted nonlinear phase term ε(t). Then, the RVP method is used to compensate the nonlinearity in the whole range profile by making the received nonlinear phase term ξ (t − τ ) range independent.

1) TRANSMITTED NONLINEARITY ESTIMATION
In order to estimate the transmitted nonlinear phase term ε(t), we firstly use two strong reflected targets at different ranges as the references. To apply them in digital domain, the Zero-IF signals of these two references are discretized: where τ stdi is the time delay of each reference and i = 1, 2 represent the reference number. Then, the non-parametric PGA algorithm, one of the most typical azimuth phase compensation methods used in SAR imaging to remove motion-induced azimuthal blurs [35]- [37], is introduced here to estimate the nonlinear phase term ϕ stdi (n) = ε(n) − ξ (n − τ stdi f s ) for each reference, which gives the advantages of no parameters required (regardless of the nonlinearity waveform) and less sensitive to noise than the direct methods in [10], [15] (demonstrated in sec III. A). There are four main steps in the process of the PGA algorithm: circular shifting, windowing, phase estimation, and iterative phase correction [35]. The schematic of the PGA algorithm in frequency domain for an unfocused echo is illustrated in Figure 3 and the specific processing procedures for the reference are provided in details: STEP 1, (Circular Shifting): The first step in PGA is to shift the isolated scatter echo from its original spot in the frequency domain to zero frequency, as illustrated in Figure 3 by frequency shifting the echo (blue line) to the zero frequency (black line). Then the frequency offset due to the target time delay is roughly eliminated from S q ZF−stdi (n) (trying to eliminate the term 2πKnτ stdi /f s in (6)), with the purpose of leaving only the nonlinear phase information and the noise. STEP 2, (Windowing): The next step is windowing the signal in frequency domain. The rectangular window filter is used in this paper to preserve the dominant unfocused echo while discarding frequency components that cannot contribute to the phase estimation. Thus, noise can be dramatically eliminated while preserving most of the scatter's information.
where s q stdi (n) is the signal in time domain after circular shifting and windowing.
where m represents the cumulative observations for the reference in order to utilize the repetitive nature of the nonlinear phase term and to reduce the impact of noise [24].
Then the nonlinear phase term in s q stdi (n) is coarsely estimated by summing over φ q stdi (n): It is noted that the bias and the linear trend should be firstly removed from theφ q stdi (n) prior to correction. STEP 4, (Iterative Phase Correction): firstly, phase correction is implemented by: is the measured Zero-IF signal S ZF−stdi (n). Then, the estimation (the above three steps) and the phase correction process should be repeated iteratively, where q represents the iteration number. As the reference's echo in the range profile becomes more focused and the side lobes more reduced after phase correction, the circular shifting more precisely removes the frequency offset and the diminishing window more effectively eliminates the noise while preserving the scatter's information. Eventually, the algorithm is driven toward convergence. Therefore, the nonlinear phase term ϕ stdi (n) = ε(n) − ξ (n − τ stdi f s ) for each reference is eliminated and thus be estimated by summing all theφ q stdi (n) in each iteration. Our experience has shown that when the standard deviation of the estimated phase term ϕ q stdi (n) drops to a few hundredth of a radian, the results will improve little with additional iterations.
Finally, with the estimation of the two references' nonlinear phase termsφ stdi (n) and the derivative approximation, the transmitted nonlinear phase term ε(n) is estimated:

2) NONLINEAR PHASE COMPENSATION FOR THE BEAT SIGNAL
With the estimation of the transmitted nonlinear phase term ε(n), the nonlinear phase term ϕ(n) = ε(n) − ξ (n − τ f s ) for the beat signal can be separated into two parts: ε(n) and ξ (n − τ f s ). Then, with the help of the RVP removal [19], [20], the dependence of the received nonlinearity term ξ (n − τ f s ) on the target's time delay τ can be removed. Finally, the nonlinearity in the whole range profile can be corrected with a unique function. The specific processing procedures for the beat signal's nonlinear phase compensation are provided in the following 3 steps, and the flowchart of the whole proposed algorithm is shown in Figure 4.

STEP 1, (Transmitted Nonlinear Phase Term Removal):
with the estimated nonlinear phase termε(n), its contribution to the beat signal nonlinearity can be eliminated by the following multiplication: (12)

STEP 2, (Range Dependent Time Shift)
: the same method used for Residual Video Phase removal [39] can be used here to introduce a range dependent time shift. The RVP filter (exp jπk 2 /K ) induces the frequency (range/ time delay) dependent phase shift, which compensates for the time delay in the residual nonlinearity (term ξ (n−τ f s ) in (12)) and makes it range independent [19,20]: where k is the discrete frequency and the range independent term ζ RVP (n) represents the response after the received nonlinearity passes through the RVP filter: STEP 3, (Nonlinearity Correction): as the range dependency of the nonlinearity has been removed, a simple multiplication will eliminate the nonlinearity in the whole range profile: After the correction, the ideal response of the target's beat signal has been recovered in S ZF4 (n). Noted that for the simplicity of the demonstration, the signal model and the correction procedure are only illustrated with a single undertest target. The method can also be used when there are multiple undertest targets in the range profile.

III. RESULTS
In this section, not only the simulated data, but also the real measured data are used to demonstrate the superiority of our method both quantitatively and qualitatively. In addition, both the direct beat signal correction method [10] and the correction method based on the signal model [19] are also shown for comparison.
The received nonlinear phase term ξ (n) was introduced as: We manipulated the difference between the transmitted and the received nonlinear phase term by adjusting the coefficient d. The reference target1 was located at 0.4 m and the reference target2 was located at 0.6 m ahead of the sensor. It worth mentioning that the direct method [10] (labeled as JET in the following of this paper) and the correction method based on the signal model [19] (labeled as RVP) also used target2 as their references.
Firstly, the nonlinear phase estimation for reference tar-get2 (ε(n) − ξ (n − τ std2 f s )) by the PGA method is demonstrated, as shown in Figure 5. The SNR (Gaussian white noise) was set to 26 dB for the beat signal in time domain and d was set to 2.5 in this case. In addition, the direct nonlinear phase measurement method used in [10], [15], where the estimated nonlinear phase equals to the measurement phase of the beat signal minus the ideal phase of the target2 (φ std2 (n) = phase [S ZF−std2 (n)] − 2π f 0 τ std2 + Knτ std2 /f s − 1/2K τ 2 std2 ), is also shown for comparison. It is noted that no repetitive observation (m = 1 in (8)) was used by the two methods for a fair comparison. Both methods had effectively estimated the nonlinear phase of target2, but due to the circular shifting and windowing procedure, the PGA method was more insensitive to the noise, as can be seen from the enlarged part in Figure 5.
As for the correction effect, we tested 3 conditions when d = 0, 2.5, and 10, as shown in Figure 6. The first column of Figure 6 show the correlation coefficient between the compensated echo and the ideal echo when the range of the single undertest target changed from 0.01 m to 1 m (0.01 m interval). In other words, the more similar with the ideal focused echo (the better the correction effect), the closer the value is to 1. The second and third column show the corrected echoes when the target located at 0.6 m and 0.9 m, respectively. It is noted that in order to evaluate the correction effect for nonlinearity, the phase noise and the thermal noise, which are the other two main factors that degrade the range profile [30], were not included in this simulation. And that, in order to reduce the side lobes, Hanning window was added before the FFT in this paper.
As we discussed in sec II. B, the direct method [10] (JET) treats the nonlinear phase term uniformly as a fixed whole part, so it is relatively insensible to the difference between ε(n) and ξ (n). However, the effectiveness of the method decays sharply with the offset distance from the reference (0.6 m for this case), as can be seen from Figure 6. The correction method based on the signal model [19] (RVP) solves the problem under the pre-assumption that ε(n) = ξ (n). Thus, the effectiveness of the method dramatically declines with the difference between ε(n) and ξ (n).
Based on the precise signal model, our method reached the outstanding results regardless of the range of the target and the nonlinearity differences. However, errors are still introduced by the inherently biased derivative approximation [20] in (11), which deserves further study.

B. REAL THz DATA VALIDATION
In this section, real data of the THz radar was used to validate the proposed correction method. The THz system was established based on the microwave up-conversion technology  with the transmitted bandwidth of 60 GHz, a center carrier frequency of 190 GHz, the PRI of 7 kHz, and the sampling rate of 250 MHz.
We firstly tested a metal plate at different ranges. The 0.4 m data (metal plate located 0.4 m ahead of the sensor) was used as the reference1 and 0.6 m data was used as the reference2 for our method. The reference2 was also used as the reference for the JET method [10] and the RVP method [19]. Figure 7 shows the original range profiles and the corresponding corrected range profiles by different methods when the range of the undertest metal plate changed from 0.38 m to 1.08 m with 0.02 m interval (36 test points). Figure 8 shows the comparison of the main echoes at several typical ranges (0.38 m, 0.6 m, and 1.08 m, respectively). In order to further quantitively evaluate the correction effect, the correlation coefficient between the calibrated echoes and the ideal echoes (simulated) is shown in Figure 9 (a), and the peak values of the echoes is shown in Figure 9 (b). Furthermore, for the ultra-wideband system, the range resolu-tion is specially emphasized. Thus, the width of the echo is evaluated by the parameter c in the fitted Gaussian function a exp − ((z − b) /c) 2 , as shown in Figure 9 (c). It is noted that some of the points (some JET and RVP points and all the ORI points) are missing in Figure 9 (c) because of the poor fitting results (R-square < 0.85: the echo is so unfocused that it is inappropriate to be evaluated by a single Gaussian function).
From these results, it can be seen that the effectiveness of all these correction methods decays with the offset distance from the reference2. However, due to the modified signal model and the correction for the range dependent nonlinearity, the decay rate of our method is much smaller than the other two methods. Moreover, the sensitivity, range resolution, and response shape of the echoes have all been improved for our method. Though the JET method reaches the best effect around the reference location because there is less approximation, the superiority is negligible, as can be shown in Figure 8 (b) and Figure 9.   In order to test the correction effect when there are multiple targets, a 4 mm thick ABS plate was inserted in front of the metal plate. The front surface of the ABS plate was located at 0.445 m and the metal plate moved from 0.48 m to 1.08 m with the interval of 0.02 m (31 test points), as shown in Figure 10. The front and back surface of the ABS plate together with the metal plate can be clearly seen in the corrected range profiles. Figure 11 shows the range profiles when the metal plate locates at 0.5 m and 1.08 m, respectively. The superiority of our correction method when the offset distance increases has also been demonstrated. It should be mentioned that due to the refractive index of the ABS sample (n≈1.56 at 200 GHz [41]) and the unperfect orthogonality between the radar and the ABS plate, measured locations of the ABS back surface and the metal plate were 5.4 mm larger than the real distance.

IV. CONCLUSION
In this paper, an effective software correction method, which comprises of transmitted nonlinearity estimation and nonlinear phase compensation, has been proposed to eliminate the nonlinearity in the beat signal of THz LFMCW radar. Compared with the two categories of software correction methods in literature [10], [18]- [21], the proposed method compensates the nonlinearity with three advantages: (1) with the extra considerations of the nonlinearity caused by THz broadband hardware devices, a modified signal model has been proposed. (2) based on this model, we proposed a practical software correction method, which integrates the advantages of these two kinds of methods, and more precisely correct the nonlinearity in the whole range profile. (3) no pre-assumption of the parametric model for nonlinearity waveform is needed in our method, which further strengthens the effectiveness of the proposed method in practical use. Finally, both the simulated data and the real THz data have been used to demonstrate the superiority of the method in qualitative and quantitative way.