3D-Subspace-Based Auto-Paired Azimuth Angle, Elevation Angle, and Range Estimation for 24G FMCW Radar with an L-Shaped Array

In this paper, a three-dimensional (3D)-subspace-based azimuth angle, elevation angle, and range estimation method with auto-pairing is proposed for frequency-modulated continuous waveform (FMCW) radar with an L-shaped array. The proposed method is designed to exploit the 3D shift-invariant structure of the stacked Hankel snapshot matrix for auto-paired azimuth angle, elevation angle, and range estimation. The effectiveness of the proposed method is verified through a variety of experiments conducted in a chamber. For the realization of the proposed method, K-band FMCW radar is implemented with an L-shaped antenna.


Introduction
Recently, many commercialized radar systems have become available for automotive, surveillance, anti-drone, medical, and personal security applications, as in [1,2]. As the need for inexpensive high-performance radar is increasing, frequency-modulated continuous waveform (FMCW) radar has become popular due to its inherent ability to utilize a large bandwidth with a low sampling rate [3][4][5]. Moreover, a homodyne receiver structure for the RF (radio frequency) frequency of FMCW radar can be easily implemented with a mixer. We call this transformation via a mixer de-chirping. A variety of estimation algorithms for FMCW radar have been developed, such as range estimation [6,7], range and Doppler estimation [8,9], and range and azimuth angle estimation [10,11]. Among the literature [6][7][8][9][10][11], there are no 3D-subspace-based algorithms for the joint estimation of range, azimuth angle, and elevation angle in FMCW radars. At the same time, there is no implementation for an FMCW radar system with an L-shaped receiving antenna array. In particular, concerning the joint estimation of azimuth and elevation angle, an L-shaped antenna structure has been proposed with two-dimensional (2D) estimation algorithms [12,13] for paired estimation of the two angles in recent studies. Although the effectiveness of the suggested algorithms [10][11][12][13] has been demonstrated through simulations, they have not been verified by the implemented system and corresponding experiments. However, the proposed 3D-subspace-based algorithm for joint estimation of range, elevation angle, and azimuth angle is verified through experiments using the implemented FMCW radar system with an L-shaped receiving array. A variety of experiments have been done to demonstrate the effectiveness of the proposed algorithm with L-shaped receiving antennas.
In this paper, a 3D-subspace-based joint azimuth, elevation, and range estimation algorithm is proposed with a 3D shift-invariant stacked Hankel matrix, which consists of one-dimensional (1D) Hankel matrices in a specific way to make use of the phase relationship between the receiving channels in horizontally placed antennas and vertically placed antennas, respectively. In addition to the 3D shift-invariant parameter estimation algorithm, a 24 GHz FMCW radar system has been implemented with transmitting lens antenna and receiving L-shaped antenna elements. The effectiveness of the proposed algorithm was verified through a variety of experiments with the implemented FMCW radar system.

System Model
The transmitted FMCW chirp signal can be modeled from [3] by: for 0 ≤ t < T sym elsewhere (1) where ω c denotes the carrier frequency, µ is the rate of change of the instantaneous frequency of a chirp signal, and T sym is the duration of chirp signals. Then, the bandwidth of the FMCW chirp signal is defined by B = µT sym /2π. Consider M far-field, non-coherent, narrowband sources impinging on the L-shaped array with 2K + 1 omnidirectional sensors as shown in Figure 1. The array on the xand z-axes consists of two uniform linear sub-arrays with element spacing d, each being composed of K + 1 elements. Let φ m , Θ m , and τ m denote the elevation angle, azimuth angle, and the time delay of the m-th target, respectively. Then, the received FMCW signals at each antenna element can be represented by where a m denote the complex amplitude for the m-th target, (p,q)∈{(K,0), ..., (1,0), (0,0), (0,1), ..., (0,K)}, M denotes the number of targets, and w p,q (t) is the additive white Gaussian noise (AWGN) signal. As in Equation (2), it is assumed that the received signals are perturbed by only an AWGN source. Thus, there is no consideration for the non-Gaussian noise sources as in [14,15]. However, in practice, non-Gaussian noises occur frequently in realistic outdoor environments. Since the proposed algorithm is developed from the assumption of AWGN noise, it is impossible for the proposed method to be applied with the non-Gaussian model directly. However, it is possible to make use of a pre-whitening technique as a preprocessing step for the proposed algorithm in the case of the non-Gaussian noise model as in [16].
The two electrical angles α m and β m in Equation (2) can be represented from [9] by where λ s denotes the wavelength of the carrier signal, d: element spacing. In FMCW radar, the received FMCW chirp signals can be easily transformed into a sinusoidal waveform via a mixer and a low-pass filter in the RF circuit [1]. We call these sinusoids beat signals. The beat signals, transformed from the received FMCW signals of (2), can be represented as in [1] by where w p,q (t) denotes the transformed AWGN signal, and z(t) denotes the local oscillator (LO) phase noise of the mixer, which is independent from receiving channel indexes (p,q). After analog-to-digital conversion, the discrete time model of (4) with sampling frequency f s = 1/T s (T s : the sampling duration) satisfying the Nyquist criterion can be derived by y p,q [n] = y p,q (nT s ) for n = 0, ..., N − 1 where N = T sym /T s . From [17], it was revealed that the LO phase noise is transferred to the output of the mixer. As depicted in Figure 2a, the LO signals of the implemented FMCW RF module are divided and transferred to the mixers of the five receiving channels. Thus, it can be said that the receiving five channels have the same phase noise z(t). Since the proposed method is designed to make use of phase difference between receiving channels for the estimation of elevation and azimuth angle, the LO phase noise cannot make any kind of disturbance on the estimation of azimuth and elevation angle due to the identical z(t) for all of the receiving channels. However, in respect of range estimation, z(t) can lead to some negative effects. In the implemented RF system, Temperature Compensate X'tal (crystal) Oscillator (TCXO), AST3TQ-50, which is used to generate the LO signal, shows the phase noise performance of −95 dBc/Hz at 10 Hz offset. Considering this phase noise characteristic, the perturbation of LO on range estimation is expected to be very small.

Proposed Algorithm
Since the proposed method has been developed for the joint estimation of elevation angle, azimuth angle, and range for FMCW radar with an L-shaped array, we propose a stacked Hankel matrix to exploit the 3D shift-invariant structure. Prior to explaining the 3D shift-invariant structure, the single shift-invariant structure in the temporal domain for range estimation is addressed with a mathematical factorization model. Then, the shift-invariant structure in the temporal domain is extended to the 3D shift-invariant structure for joint estimation of elevation, azimuth, and range in the spatial and temporal domains.

Shift-Invariant Structure for Range
For convenience in expressing the single shift-invariant structure for elevation, azimuth, and range, we assume that the received signals are not perturbed by AWGN. Noise perturbation is addressed later in conjunction with singular value decomposition (SVD).
Using the beat signals of the p-th and q-th antennas, y p,q [n] for n = 0, ..., N − 1, the Hankel snapshot matrix can be defined as where L r and L c = N − L r + 1 are the selection parameters, which satisfy the conditions L r ≥ M and L c ≥ M. Without considering the perturbation by AWGN, the Hankel snapshot matrix Y p,q can be factorized as in [18] by In (7), κ m is the delay-induced phase shift between the adjacent samples, such that In (6), the diagonal matrices H and R p,q are composed of the phase terms of the beat signals of (4). The Vandermonde structured matrices A and B are defined in terms of the phase shifts κ 0 , ..., κ M− 1 , which are not changed by array indices p and q.
From Y p,q of (5), the two sub-matrices Y p,q,0 and Y p,q,1 can be defined by Comparing Y p,q,0 with Y p,q,1 , it can be easily seen that Y p,q,1 is the shifted version of Y p,q,0 in the column direction. Based on the factorization model of (6), the sub-matrices Y p,q,0 and Y p,q,1 can also be factorized, such that In (11), it is noteworthy that factorizations for Y p,q,0 and Y p,q,1 are the same except for the diagonal matrix Σ of Y p,q,1 . We call this a shift-invariant structure. This single shift-invariant relationship between Y p,q,0 and Y p,q,1 in (11) is extended to a 3D shift-invariant structure for the joint estimation of range with elevation angle and azimuth angle.

Shift-Invariant Structure for Two Electrical Angles
The stacked Hankel snapshot matrix for the joint estimation of range and azimuth and elevation angles can be defined by In (13), the matrix Y 0,0 takes the role of reference for pairing afterwards. Prior to explaining the 3D shift-invariant structure in Y, the shift-invariant structure for two electronic angles α m and β m is addressed.
Concerning α m , the two stacked matrices, which are sub-matrices of Y X , can be defined as Then, based on the factorization of (6), Y X_0 and Y X_1 can be rewritten as Since the diagonal matrix R p,q is composed of the two electronics angles α m and β m , it is not convenient to carry out further factorization in (15). Thus, we define the two diagonal matrices Θ and Φ, whose elements are only composed of α m or β m , respectively, such that Then, we can derive the following relationship: R k,0 = R k−1,0 Θ. Using this, Y X_0 and Y X_1 of (15) can be reformulated as Applying the generalized eigenvalue decomposition (EVD) relationship for Y X_0 and Y X_1 such that (Y X_1 − λY X_0 )β = 0, where λ and β denote the eigenvalue and the corresponding eigenvector, respectively, this EVD can be redefined using (18) by In (19), it is straightforward that the rank of (Y X_1 − λY X_0 ) will be M unless λ = α m for m = 0, 1, ..., M − 1, since M sources are assumed in (2). However, when λ = α m for m = 0, ..., M − 1, the rank of (Y X_1 − λY X_0 ) will be M − 1 since one of the rows in [(Θ − λI M )] will be zero. This means that we can obtain α m as a result of the generalized EVD of Y X_1 and Y X_0 .
Similar to the derivation of the EVD relationship for α m in Equations (14)- (19), the EVD relationship for β m can be derived as follows. Two matrices, which are also sub-matrices of Y Y , can be defined as Concerning β m , it satisfies R 0,k = R 0,k−1 Φ. Then, Y Y_0 and Y Y_1 of (20) can be factorized and reformulated by Applying the generalized EVD relationship for Y Y_0 and Concerning the generalized EVD form about κ m for the given stacked matrix Y, two sub-matrices Y 0 and Y 1 can be defined to make use of the shift-invariant structure in (11) and (12) by where the selection matrices.
where I k denotes the identity matrix of k by k, and 0 m×n denotes the zero matrix of m by n. Then, the generalized EVD relationship for Y 0 and Y 1 such that (Y 1 − λY 0 )β = 0 can lead to λ = κ m for m = 0,1, ..., M − 1 based on (11).

Signal Subspace
In the preceding section, a noiseless environment in Y was assumed for the derivation of the EVD relationships in association with κ m , α m , and β m . Considering perturbation of AWGN on the received signals, first, singular value decomposition (SVD) should be applied with Y to separate the signal subspace from the noise subspace. In the presence of AWGN, the stacked matrix Y can be factorized by SVD into signal and noise subspaces as in [19], such that where U and V are the orthogonal matrix, Σ denotes the diagonal matrix with singular values in the diagonal elements, U s , Σ s , and V s are associated with the signal subspace, and U n , Σ n , and V n are associated with the noise subspace. Although U s and U n are modeled to be separated with each other after SVD in (24), the signal subspace and the noise subspace cannot be separated with each other through only SVD. Generally, all of the subspace-based estimation algorithms [12][13][14][15][16]18,19] essentially involve the subspace classification step, in which the signal subspace and the noise subspace are separated from each other. In this step, SVD or EVD is used for matrix decomposition and its singularor eigenvalues are used with the Akaike Information Criterion (AIC) or the Minimum Descriptive Length (MDL) [20] for subspace separation. At this time, a high signal-to-noise ratio (SNR) is required with the subspace-based algorithms, which use the signal subspace or the noise subspace for correct subspace classification as in [18]. Thus, in a low SNR situation, performance degradation can be incurred by incorrect subspace classification. Let us define a steering matrix Y S , which shows the shift-invariant structure in Y, such that Then, it satisfies [21] In (28), J X_0 and J X_1 are selection matrices, which can extract the sub-matrices from U s corresponding to Y X_0 and Y X_1 , respectively, such that In (29), J Y_0 and J Y_1 are selection matrices, which can extract the sub-matrices from U s corresponding to Y Y_0 and Y Y_1 , respectively, such that In the preceding section, noiseless samples are assumed and used to derive the EVD relationships, for κ m , α m , and β m , respectively. However, noiseless samples cannot be made in a real RF system, so the sub-matrices from the signal subspace can substitute for the noiseless sub-matrices in three kinds of EVD forms for κ m , α m , and β m , as follows:

Low-Complexity Pairing
The preceding section addressed how each of the three parameters α m , β m , and κ m can be estimated from U s in (32)-(34). However, how those estimates are paired to each other was not addressed, which must be handled for joint range, azimuth, and elevation estimation. In the literature [22][23][24][25][26][27], 2D parameter estimation algorithms have been proposed for angle and frequency in [22,23], range and Doppler in [24,25], and azimuth and elevation in [26,27]. All of the algorithms have their own mechanisms to pair the two kinds of estimated parameters for each target. In particular, to reduce the operation cost of paired estimation, computationally efficient algorithms have been proposed in [28,29]. However, all of the algorithms in [22][23][24][25][26][27][28][29] are limited to two-dimensional parameter estimation. In particular, when the number of parameters to be matched increases from two to three, the computational complexity greatly increases for 3D matching. In this paper, a low-complexity matching method using 1D EVD is proposed.
The generalized EVD forms in Equations (32)-(34) can be rewritten in the EVD form with a pseudo-inverse, such that G = U † s,0 U s,1 = BΣB −1 (35) where G, G X , and G Y are M-by-M full-rank matrices, the superscript † denotes the corresponding Hermitian conjugate of a matrix, and B, Using the obtainedB, we performed diagonalization on G and G Y : Through the derivation of (38)-(40), the values of estimatedκ m ,α m , andβ m for the m-th target can be found in the diagonal line of matricesΣ,Θ, andΦ in ascending order, respectively. That is, for each m, the estimated values {κ m ,α m ,β m } are matched well for the m-th target by low complexity. Based on (3) and (8), the estimated elevation angle, azimuth angle, and the time delay of the m-th target can be expressed asφ where asin ( angle(α m )λ s 2πd cos(φ m ) ) denotes the inverse sine function, acos ( angle(β m )λs 2πd ) denotes the inverse cosine function, and angle (β m ), angle (α m ), angle (κ m )denotes the angle in radian, respectively.

Complexity Analysis
Since the proposed method is designed to estimate 3D parameters, such as range, azimuth, and elevation, it is fair to compare the computational complexity of the proposed method with that of the conventional 3D parameter estimation methods of an FMCW radar system. However, it was difficult to find relevant research results about 3D parameter estimation algorithms for an FMCW radar system, while 2D parameter estimation algorithms have been suggested for FMCW radar as in [11,30]. Thus, in this article, the computational complexity of the proposed algorithm is analyzed and compared with that of 2D multiple signal classification (MUSIC) of an FMCW radar system in [11,30]. The costs for the required operations are summarized in Table 1.

Operation Description
Computational Complexity For the given data matrix Y of (13), the computational complexity costs for 2D MUSIC of [11,30] can be derived to be O( where b is the iteration number for two-dimensional searching. In general, the iteration number b is set to be much bigger than K, L c , and M for a high-resolution 2D pseudo-spectrum. Thus, using b >> K, L c , M, the derived complexity cost for 2D MUSIC can be simplified to O(b 2 K 2 L 2 c ). For the proposed method, the computational complexity can be derived to O(K 2 L c + KL 2 c + L c 3 + M 3 + 2M 2 (L c -1) + 2M 2 (K -1)), which can be approximated using of 2D MUSIC, we conclude that the cost of the proposed algorithm for 3D parameter estimation is much less than that of 2D MUSIC.

Implementation of 24-GHz L-Shaped Radar
To verify the effectiveness of the proposed method with L-shaped receiving antennas, a 24-GHz FMCW radar system with L-shaped receiving antennas was implemented. The implemented system consists of an antenna, an RF circuit, an IF circuit, a data-logging platform, and data-logging software. The data-logging platform takes the role of gathering received IF samples and passing them to the data-logging software in a personal computer (PC). The logged data is saved as files. Then, the proposed method, realized in MATLAB (The MathWorks, Natick, MA, USA), is applied with the saved data in the PC.

Transmitting and L-Shaped Receiving Antennas
The transmitting antenna is composed of two lens antennas, which have 14-dBi and 19-dBi antenna gain. The 19-dBi high-gain antenna is selected as the default, and it is possible to change to the low-gain antenna by using the data-logging platform. The L-shaped receiving antenna is composed of five micro-strip patch antennas, designed to have about 6-dBi antenna gain, fabricated on an Ro4003 substrate with 0.012 inch thickness. The five receiving antennas have feeding lines with the same electric length and a meander structure as shown in Figure 3.  Figure 3 shows the configuration of the transmitting and receiving antennas; the transmitting antenna is located on the upper side of system, and the receiving antenna is on the lower side. The RF front module is designed to be able to obtain the signals from the receiving antennas through a rectangular waveguide interface.
The receiving patch antennas are set in an L-shaped formation for simultaneous elevation and azimuth angle estimation. The two electrical angles α m and β m of Equation (3) can be estimated using receiving antennas 1 to 3 and receiving antennas 3 to 5, respectively, as shown in Figure 1. Each patch is located by 1 wave length distance vertically or horizontally, resulting in a field of view of ±30 • in azimuth and elevation.
The antenna operates over the frequency range of 24.025 to 24.225 GHz for the FMCW radar signal bandwidth as shown in Figure 4a. Figure 4b shows the simulated receiving antenna radiation pattern; in this case, the antenna is designed to face forward. The return loss and radiation pattern in Figure 4 are simulated together with the waveguide, transition, and meander feeding line. The antenna characteristic is optimized for radar performance with adjacent antennas and discontinuity of structure, although it is possible to improve the antenna characteristic. The transmitting antennas are designed by the lens to have a 12 • beam width for the high-gain antenna and a 36 • beam width for the low-gain antenna, and the side-lobe levels of both antennas are lower than 15 dB as shown in Figure 5. The pattern was measured from the effective isotropic radiated power and normalized by the calculated transmit power by the RF system output power. The measured equivalent isotropically radiated power (EIRP) values were 39.3 dBm for the high-gain antenna and 34.1 dBm for the low-gain antenna.

24-GHz Transceiver and IF
The 24-GHz FMCW RF system was implemented with a frequency synthesizer, a phase-locked loop (PLL) circuit, a voltage-controlled oscillator (VCO), an edge-coupled filter for band-pass filtering, a power amplifier (PA), and a low-noise amplifier (LNA). The overall structures of the RF and IF systems are described in Figure 2a.
Since we assume a saw-tooth wave as in (1), a PLL chip was set to generate a saw-tooth wave pattern with a 100-us period and a 200-MHz bandwidth over the range of 24.025 to 24.225 GHz with 10 dBm VCO output power. The reference clock of PLL is 20 MHz, which comes from TCXO. The VCO output is filtered with an edge-coupled filter to improve the side-lobe characteristic. As shown in Figure 2b, the edge-coupled filter is implemented with a center frequency of 24.125 GHz, a 10% bandwidth, and a Chebyshev-type band pass filter. The parameters of the implemented RF system are summarized in Table 2. The filtered VCO signal is divided for transmission (TX) and LO as shown in Figure 2b. To increase the power of the divided VCO signal, a PA with a 20-dBm output P1dB is used. Since a WR34 interface is used with transmitting and receiving antennas, a fin-line structure transition is implemented to transform between the micro-strip line and WR34. The front side of the fin-line transition is connected to WR34, i.e., the signals on the micro-strip line are rotated 90 • as shown in Figure 2b.

Data-Logging Platform
The radar-logging platform was implemented to transfer maximum 16 channel analog digital converter (ADC) input signals to the PC in real-time. This platform mainly consists of digital signal processor (DSP) and field programmable gate array (FPGA) chips, namely TI TMS320C6455 and Stratix EP4SE530. The ADC output signals of a maximum of 16 channels are first saved in the first in first out (FIFO) memory of the FPGA and then transferred to the double-data-rate two synchronous dynamic random access memory (DDR2 SDRAM) through a direct memory access mechanism. The saved data can be handled by the DSP chip or transferred to the PC through 1 G local area network (LAN) communication as shown in Figure 6. In this work, the proposed algorithm is implemented in MATLAB and applied with the data saved in the PC by the data-logging platform. The specifications of the developed data-logging platform are summarized in Table 3. Figure 6. Photograph of the developed data-logging platform. IO: input output, CPLD: complex programmable logic device, ADC: analog to digital converter, DSP: digital signal processor, DDR2: double data rate two, FPGA: field programmable gate array.

Experiments
This section describes several indoor experiments, which were conducted to detect one target, two targets, and four targets, respectively, in the microwave anechoic chamber of Daegu Gyeongbuk Institute of Science & Technology (Daegu, Korea). One outdoor experiment was also carried out to verify the performance of the proposed algorithm in a realistic scenario.
A radar module with an L-shaped antenna array was set up on one side of a pillar, and the location of the radar module could be adjusted by a motor at the bottom of the pillar as shown in the left part of Figure 7a. We configured the Cartesian coordinate system following the model in Figure 1 by arranging the pillar along the z-axis, the L-shaped array in the x-z plane, and the inflection point of the L-shaped antenna array as the origin. Four square iron blocks with a 10-cm side length were selected as targets, and every target was mounted on one railway. The scenarios of the experiments for detecting one target, two targets, and four targets are shown in Figure 7, and the settled positions of the targets for each experiment are shown as the form (x, y, z) in the following list: We conducted 100 trials for each experiment, and the estimated positions are illustrated in Figure 8a-e, respectively. Root mean square error (RMSE) was selected as the measurement of the experiment results, and we defined the RMSE for each target as 1 100 denotes the estimated position of the target (x, y, z) from the l-th trial. The measured RMSE values for all experiments are shown in Table 4.   In our experiments, it was firstly found that all of the five receiving channels have their own initial phases, which are different to each other due to implementation loss in the 24-GHz RF circuit. Therefore, without the estimation and compensation for the initial phases, the estimation results of the azimuth and elevation angles were biased a lot, resulting in a lot of RMSE. Thus, we tried to estimate the initial phases for each channel and use them for a calibration factor before azimuth and elevation estimation. After that, the estimation results shows less spread and less RMSE.
In the preceding section, the effect of low SNR is addressed in conjunction with subspace classification. Unlike the high SNR environment of the targets within 6 m in a chamber, the targets at a maximum of 30 m in outdoor environments are under a low SNR situation and can show performance degradation due to a low SNR and the effect of clutter.
As shown in Figure 9a, the radar module was set up on a holder, and the inflection point of the L-shaped antenna array was arranged as the origin at a height of 1 m from the ground. We then arranged three targets: (1) Iron target in (−0.8 m, 10 m, 0 m); (2) Human target in (1 m, 30 m, 0 m); (3) Human target in (1.5 m, 13 m, 0 m). We tried to arrange the center of the targets at around 1 m height from the ground, and hence we put '0 m' here for the three targets. It can be noted that the ground is besprinkled with little stones with various shapes, which contribute much interference for the estimations. We conducted 100 trials for the outdoor experiment, and the estimated positions are illustrated in Figure 9b. The corresponding estimation results for each target are encompassed by the corresponding ellipse. For the first and third targets, biases in estimation results occurred due to the effect of non-Gaussian noise in the real environment. For the estimation results of the second target at a range of 30 m, the miss detection rate really increased a lot due to a low SNR. Except for the estimation results from the three targets, many undesired results were obtained from the clutter, mainly from the stones on the ground in this experiment. In the realistic scenarios, the performance of the proposed algorithms degraded due to the non-Gaussian noises and the low SNR.

Conclusions
This paper has proposed an auto-matched range-angle-Doppler 3D estimation algorithm based on a 3D shift-invariant structure and presented an FMCW radar system with L-shaped receiving antennas. Our experimental results in a chamber demonstrated that the implemented L-shaped FMCW radar system with the proposed algorithm achieved high-quality joint estimation of range, azimuth angle, and elevation angle. However, the results of the outdoor experiments demonstrated that the performance of the proposed algorithm degraded under non-Gaussian noises and the low SNR situation. Moreover, since the proposed method is composed of a variety of matrix operations such as SVD and EVD, it was difficult to implement the algorithm on FPGA and DSP for real-time estimation. Instead, the proposed algorithm is realized in MATLAB and applied with the receiving data of the experiments.