Analyze the FMCW Waveform Skin Return of Moving Objects in the Presence of Stationary Hidden Objects Using Numerical Models

: In this paper, a high-performance antenna array system model is presented to analyze moving-object-skin-returns and track them in the presence of stationary objects using frequency modulated continuous wave (FMCW). The main features of the paper are bonding the aspects of antenna array and electromagnetic (EM) wave multi-skin-return modeling and simulation (M&S) with the aspects of algorithm and measurement/tracking system architecture. The M&S aspect models both phase and amplitude of the signal waveform from a transmitter to the signal processing in a receiver. In the algorithm aspect, a novel scheme for FMCW signal processing is introduced by combining time-and frequency-domain methods, including a vector moving target indication ﬁlter and a vector direct current canceller in time-domain, and a constant false alarm rate detector and a mono-pulse digital beamforming angle tracker in frequency-domain. In addition, unlike previous designs of using M × N fast Fourier transform (FFT) for an M × N array, only four FFTs are used, which tremendously save time and space in hardware. With the presented model, the detection of the moving-target-skin-return in stationary objects under a noisy environment is feasible. Therefore, to track long range and high-speed objects, the proposed technique is promising. Using a scenario having (1) a target with 17 dBm 2 radar cross section (RCS) at about 40 km range with 5.936 Mach speed and 11.6 dB post processing signal to noise ratio, and (2) a strong stationary clutter with 37 dBm 2 RCS located at the proximity of the target, it demonstrates that the root-mean-square errors of range, angle, and Doppler measurements are about 26 m, 0.68 degree, and 1100 Hz, respectively.


Introduction
Detecting and tracking moving targets' locations, Doppler frequencies, and velocities based on the target-skin-return is the main function of modern electromagnetic (EM) based tracking and measurement systems, e.g., radars [1]. Successful EM propagation and antenna array modeling for this detection and tracking process in software is important for the designing of such systems for radar electronic warfare studies. In this paper, we present the use of numerical models for analyzing the EM skin return of moving targets in the presence of stationary clutters and receiver noise with Simulink software. The main models include the "point-reflector" target simulation model, the simplified radio frequency (RF) frontend, and the digital beamforming (DBF) receiver. The "point-reflector" target simulation model simulates a target at certain three-dimensional (3D) coordinates using a radar cross section (RCS) number, the 3D coordinates, arbitrary 3D velocities, and arbitrary 3D accelerations as inputs to this model. The simplified RF frontend models transmit-to-receive (TX-to-RX) coupling (or leakage) and the additive white Gaussian noise (AWGN). The features of the paper are to tie the aspects of antenna array and EM wave multi-skin-returns (target and clutter) modeling and simulation (M&S) with the The radar equation used to model the power-level propagation effect is the bi-static propagation equation. This was used to model the radar received power, as a function of the distance between the radar and its target. The equation is where is the received signal power in dBm, is the transmitted signal power in dBm, is the transmit antenna gain of the radar in dBi, is the receive antenna gain of the radar in dBi, is in dBm , is the wavelength of the carrier frequency, is the distance between the radar transmitter and the target in meters, is the distance between the radar receiver and the target in meters, and is the loss in dB. Note that, since targets move, all the parameters in Equation (1) can change with time.
It is also assumed that the RF propagation channel between the TX antenna and each receiving antenna element in the phased array has all-pass characteristics with additive white Gaussian noise (AWGN). This assumption implies that we can simply model the RF path with unit gain, zero delay from the phase reference point of TX antenna to the phase reference center of the RX phased array, and also, the receiving phased array is calibrated to remove all the differences of extra delay, gain, and phase of the RF paths among × receiving channels.
To design and verify the system in computer simulations, we need to model and simulate continuous target movement effects. This is done by using variable delay filters operated on the transmitted waveform, based on the observations from Equation (3). Because of the band-limited property of the transmitted waveform, we only need to delay the in-band signal by a variable delay corresponding to each sample based on the 3D relationship between TX/RX antenna phase centers and the target location, which changes with time.

RX Coordinate System
The coordinate system in the M&S of the FMCW tracking and measurement system is presented in Figure 2. The panel of the RX phased array is in the XZ-plane, and its phase reference center is at the origin of the Cartesian coordinate. A target is located at ( , , ), which is also known as azimuth-angle (Az) (in degrees), elevation-angle (El) (in degrees), and range ( ) (in meters). Since the target moves, the ( , , ) values change continuously with time. The radar equation used to model the power-level propagation effect is the bi-static propagation equation. This was used to model the radar received power, as a function of the distance between the radar and its target. The equation is P r = P t + G tR + G rR + σ + 10·log 10 λ 2 (4π) 3 where P r is the received signal power in dBm, P t is the transmitted signal power in dBm, G tR is the transmit antenna gain of the radar in dBi, G rR is the receive antenna gain of the radar in dBi, σ is in dBm 2 , λ is the wavelength of the carrier frequency, R tx is the distance between the radar transmitter and the target in meters, R rx is the distance between the radar receiver and the target in meters, and L is the loss in dB. Note that, since targets move, all the parameters in Equation (1) can change with time.
It is also assumed that the RF propagation channel between the TX antenna and each receiving antenna element in the phased array has all-pass characteristics with additive white Gaussian noise (AWGN). This assumption implies that we can simply model the RF path with unit gain, zero delay from the phase reference point of TX antenna to the phase reference center of the RX phased array, and also, the receiving phased array is calibrated to remove all the differences of extra delay, gain, and phase of the RF paths among M × N receiving channels.
To design and verify the system in computer simulations, we need to model and simulate continuous target movement effects. This is done by using variable delay filters operated on the transmitted waveform, based on the observations from Equation (3). Because of the band-limited property of the transmitted waveform, we only need to delay the in-band signal by a variable delay corresponding to each sample based on the 3D relationship between TX/RX antenna phase centers and the target location, which changes with time.

RX Coordinate System
The coordinate system in the M&S of the FMCW tracking and measurement system is presented in Figure 2. The panel of the RX phased array is in the XZ-plane, and its phase reference center is at the origin of the Cartesian coordinate. A target is located at (φ, θ, R rx ), which is also known as azimuth-angle (Az) (in degrees), elevation-angle (El) (in degrees), and range (R rx ) (in meters). Since the target moves, the (φ, θ, R rx ) values change continuously with time.

Received Signal Model
The observable range for radar measurement at time is defined as the average of ' and ' , i.e., ( ) = and = + , where ' is the instance of EM wave that hits the target after the wave leaves the TX phase-center, and is the time interval that is needed when the target echo signal travels backs to the RX phasecenter. Where, is the free space speed of light. Note that, here we assume that the phasecenters of the TX and RX antennas are not collocated. After the two-way propagation, the target echo signal ( ) is intercepted by the elements in the RX phased array. At the ( , ) th element, the received analog signal in the complex format can be expressed as where is the complex gain coefficient of the ( , ) th element, and ( ) is the transmitted waveform. The received signal is then directly down-converted by the RF front-end, and digitized by ADC at the sampling rate of . The sampled signal is where is the time sample step, = 1/ . The complex gain coefficient can be detailed as where is the target RCS, is the propagation loss defined by the radar equation in Equation (1), ℎ is the RX RF channel gain, and the beamforming coefficient is in Equation (13), which defines the required phase from the ( , ) th element in the array in order to form the antenna beam in ( , ) or [ (0), (1)] = [cos( ) sin( ) , cos( )] direction, as described in Equation (12). We have assumed that oversampling is done at the ADC. At the receiver, the ADC I/Q samples are down-sampled by a factor of to the baseband sampling frequency , where = .

Waveform Interpolation
As shown Equation (3), the round-trip time delay

Received Signal Model
The observable range for radar measurement at time t is defined as the average of c , where t is the instance of EM wave that hits the target after the wave leaves the TX phase-center, and is the time interval that is needed when the target echo signal travels backs to the RX phase-center. Where, c is the free space speed of light. Note that, here we assume that the phase-centers of the TX and RX antennas are not collocated. After the two-way propagation, the target echo signal r(t) is intercepted by the elements in the RX phased array. At the (m, i) th element, the received analog signal in the complex format can be expressed as where g c mi is the complex gain coefficient of the (m, i) th element, and x(t) is the transmitted waveform. The received signal is then directly down-converted by the RF front-end, and digitized by ADC at the sampling rate of f S . The sampled signal is where k is the k th time sample step, T S = 1/ f S . The complex gain coefficient can be detailed as where σ is the target RCS, g air is the propagation loss defined by the radar equation in Equation (1), h c r f is the RX RF channel gain, and the beamforming coefficient b c mi is in Equation (13), which defines the required phase from the (m, i) th element in the array in order to form the antenna beam in (φ, θ) or [a(0), a(1)] = [cos(φ) sin(θ), cos(θ)] direction, as described in Equation (12). We have assumed that oversampling is done at the ADC. At the receiver, the ADC I/Q samples are down-sampled by a factor of K to the baseband sampling frequency f B , where A filter with the approximated all pass characteristics in the passband of − 1 2K , 1 2K means that the filter has unit gain and linear phase slope within [−π, π] in its passband. Using the filter in the oversampling domain, we can do interpolation to obtain the required signal phase that relates to the time delay of τ f . This technique is called Lagrange interpolation filter, which can be expressed as where D = τ f + floor P 2 , k is the k th tap of the Lagrange filter, and P is the order of the Lagrange filter. For example, when K = 2, the interpolation filter for a fractional delay of 0.25 has a gain and phase response as shown in Figure 3.
A filter with the approximated all pass characteristics in the passband of − , means that the filter has unit gain and linear phase slope within [− , ] in its passband. Using the filter in the oversampling domain, we can do interpolation to obtain the required signal phase that relates to the time delay of . This technique is called Lagrange interpolation filter, which can be expressed as where = + floor , ' is the ' th tap of the Lagrange filter, and P is the order of the Lagrange filter. For example, when = 2, the interpolation filter for a fractional delay of 0.25 has a gain and phase response as shown in Figure 3. Before applying the reflection/scattering gain to the skin return, let us illustrate the post Lagrange interpolation waveform amplitude and phase vs. the fractional delay numbers due to the target movement in Figures 4 and 5. We compare these plots with the post nearest sample interpolation waveform amplitude and phase vs. the fractional delay numbers in Figures 6 and 7. The amplitude and phase in Figures 4 and 5 are computed from the complex differential waveform ( ) based on the following equation: where (•) * stands for conjugate operation. The complex differential waveform ( ) represents the phase and amplitude changes of − . Clearly, Before applying the reflection/scattering gain to the skin return, let us illustrate the post Lagrange interpolation waveform amplitude and phase vs. the fractional delay numbers due to the target movement in Figures 4 and 5. We compare these plots with the post nearest sample interpolation waveform amplitude and phase vs. the fractional delay numbers in Figures 6 and 7. The amplitude and phase in Figures 4 and 5 are computed from the complex differential waveform r di f f (kT S ) based on the following equation: where (·) * stands for conjugate operation. The complex differential waveform r di f f (kT S ) represents the phase and amplitude changes of x kT S − 2R(kT s ) c over one T S duration.
Using r di f f (kT S ), we can more clearly observe the continuity of the signal x kT S − 2R(kT s ) c Electronics 2021, 10, 28 7 of 30 in amplitude and phase. If Lagrange interpolation is used, we reconstruct x kT S − 2R(kT S ) c using the Lagrange filtered x(kT S ). If the nearest sample interpolation is used, we approxi- T S c T S . Clearly, the Lagrange interpolation functions well, with no discontinuity in phase and amplitude after interpolation. The nearest sample interpolation shows an occasional phase discontinuity at 1332.5 T S delay time. Note that the phase discontinuity implies that the simulation with the nearest sample interpolation is inaccurate at the phase jump instance.
Electronics 2021, 10, x FOR PEER REVIEW 7 of 29 after interpolation. The nearest sample interpolation shows an occasional phase discontinuity at 1332.5 delay time. Note that the phase discontinuity implies that the simulation with the nearest sample interpolation is inaccurate at the phase jump instance.   Electronics 2021, 10, x FOR PEER REVIEW 7 of 29 after interpolation. The nearest sample interpolation shows an occasional phase discontinuity at 1332.5 delay time. Note that the phase discontinuity implies that the simulation with the nearest sample interpolation is inaccurate at the phase jump instance.

Digital Phased Array
It is assumed that the phased array has and elements in X-and Z-direction, respectively, as shown in Figure 8. In our design, the array has uniform element spacing ( ) in both X-and Z-directions. Without loss of generality, it is also assumed and are even numbers. Based on the microwave antenna array theory, the array factor in the direction of ( , ) can be written as the following, if all the elements are uniformly excited with unit amplitude: where ( , ) is the wave-vector in the direction of ( , ), and ( , ) is the displacement from the ( , ) element to the origin of the coordinate system, they can be expressed as

Digital Phased Array
It is assumed that the phased array has and elements in X-and Z-direction, respectively, as shown in Figure 8. In our design, the array has uniform element spacing ( ) in both X-and Z-directions. Without loss of generality, it is also assumed and are even numbers. Based on the microwave antenna array theory, the array factor in the direction of ( , ) can be written as the following, if all the elements are uniformly excited with unit amplitude: where ( , ) is the wave-vector in the direction of ( , ), and ( , ) is the displacement from the ( , ) element to the origin of the coordinate system, they can be expressed as

Digital Phased Array
It is assumed that the phased array has M and N elements in X-and Z-direction, respectively, as shown in Figure 8. In our design, the array has uniform element spacing (d c ) in both X-and Z-directions. Without loss of generality, it is also assumed M and N are even numbers. Based on the microwave antenna array theory, the array factor in the direction of (φ, θ) can be written as the following, if all the elements are uniformly excited with unit amplitude: where K(φ, θ) is the wave-vector in the direction of (φ, θ), and r (m, i) is the displacement from the (m, i) th element to the origin of the coordinate system, they can be expressed as Electronics 2021, 10, 28 where λ is the free space wavelength, and since the array is in the XZ-plane, where x(m) = d c m − M 2 + 0.5 and z(i) = d c i − N 2 + 0.5 . So, Equation (8) can be rewritten as where [a(0), a(1)] = [cos(φ)sin(θ), cos(θ)].
Electronics 2021, 10, x FOR PEER REVIEW 9 of 29 where λ is the free space wavelength, and where The term of 2 (0) • − + 0.5 + (1) • − + 0.5 in Equation (11) can also be considered as the phase needed on the ( , ) th element in order to let the array form a beam in the ( , ) direction in digital-domain.

FMCW Triangular Waveform
The chirp signal uses the triangular waveform. The triangular shape depicts the change of the chirp signal's frequency component over time. The frequency of transmitted waveform versus time is: where is the carrier frequency of the chirp signal, and stands for the slope of the chirp signal, while stands for the upbeat or downbeat time duration as illustrated in Figure 9. The transmitted waveform in a closed form is The term of 2π d c λ a(0)· m − M 2 + 0.5 + a(1)· i − N 2 + 0.5 in Equation (11) can also be considered as the phase needed on the (m, i) th element in order to let the array form a beam in the (φ, θ) direction in digital-domain. We define the beamforming coefficient of the (m, i) th element in the array as

FMCW Triangular Waveform
The chirp signal uses the triangular waveform. The triangular shape depicts the change of the chirp signal's frequency component over time. The frequency of transmitted waveform versus time is: where f c is the carrier frequency of the chirp signal, and b stands for the slope of the chirp signal, while T stands for the upbeat or downbeat time duration as illustrated in Figure 9. Electronics 2021, 10, x FOR PEER REVIEW 10 of 29 Figure 9. The triangular FMCW waveform.
When this RF signal hits a moving target with the radial speed ( ) (m/s) at a distance ( ) relative to the transmitter, it bounces back to the receiver with the radial speed ( ) (m/s) at a distance ( ) relative to the receiver. The received signal can be expressed in a classical physics sense as: where and are the gain and phase introduced by the EM wave reflection/scattering and other losses in the two-way path of TX (antenna)-target-RX (antenna), respectively.
is the time of receiving the signal at the radar receiver, = + and ' is the time instance that reflection/scattering occurs, and where , is the initial distance between the transmitter and the target, and , is the initial distance between the receiver and the target. Note that the observable range for radar measurement is defined as ( ) = , and we assume ( ) < , so that the range ambiguity is not considered. We choose the FFT windows centered at and as shown in Figure 9 to estimate the Doppler values and ranges with both the upbeat and the downbeat waveform. We assume that the transmitter and the receiver are powered up at the same time, i.e., time 0, thus and are with respect to this power-up time. In this paper, we choose = 2 + and = 2 + as the centers of the FFT windows, and assume that the FFT window has a width of . Here, stands for the 2 period. With this parameter selection, the guard time is and can guard max one-way propagation time of .
By taking derivatives to the phases at = 2 + and = 2 + in Equation (16) using the chain rule of derivatives, where stands for the 2 period, we can estimate the frequency of the received waveform in the selected FFT windows as: The transmitted waveform in a closed form is When this RF signal hits a moving target with the radial speed v tx (t ) (m/s) at a distance R tx (t ) relative to the transmitter, it bounces back to the receiver with the radial speed v rx (t ) (m/s) at a distance R rx (t ) relative to the receiver. The received signal can be expressed in a classical physics sense as: where g and ϕ are the gain and phase introduced by the EM wave reflection/scattering and other losses in the two-way path of TX (antenna)-target-RX (antenna), respectively.
t is the time of receiving the signal at the radar receiver, t = t + R rx (t ) c and t is the time instance that reflection/scattering occurs, and (17) where R tx,0 is the initial distance between the transmitter and the target, and R rx,0 is the initial distance between the receiver and the target. Note that the observable range for radar measurement is defined as , and we assume R(t) < cT 4 , so that the range ambiguity is not considered.
We choose the FFT windows centered at t 1 and t 2 as shown in Figure 9 to estimate the Doppler values and ranges with both the upbeat and the downbeat waveform. We assume that the transmitter and the receiver are powered up at the same time, i.e., time 0, thus t 1 and t 2 are with respect to this power-up time. In this paper, we choose t 1 = 2nT + 2 3 T and t 2 = 2nT + 5 3 T as the centers of the FFT windows, and assume that the FFT window has a width of T 3 . Here, n stands for the n th 2T period. With this parameter selection, the guard time is T 2 and can guard max one-way propagation time of T 4 .
By taking derivatives to the phases at t 1 = 2nT + 2 3 T and t 2 = 2nT + 5 3 T in Equation (16) using the chain rule of derivatives, where n stands for the n th 2T period, we can estimate the frequency of the received waveform in the selected FFT windows as: The FFT windows centered at t 1 and t 2 are properly guarded before and after the windows in time so that they only contain the skin return signals in upbeat or downbeat, respectively. This avoids the difficulty in spectrum analysis due to having received signals in both upbeat and downbeat.
We also define the observable radial speed for radar measurement at time t as The quantity R(t) and v(t) will be estimated in the simulation. Since the transmitter and the receiver are very close to each other and the speed of the target is much less than the speed of light, we have Subtracting the frequencies of the transmitted waveforms from f u (t) and f d (t), we obtain: After plugging in t 1 = 2nT + 2 3 T and t 2 = 2nT + 5 3 T, we obtain: The Doppler frequency is thus estimated by By properly choosing the parameters b and T, and by ignoring the terms of 6c since they are very small, as v(t 1 ) and v(t 2 ) are much smaller than the speed of light c, we obtain:

The Receiver Diagram
The digital portion of the receiver diagram is shown in Figure 10. It takes complex I/Q samples after ADCs from a 4 × 4 array and filters them with a root-raised-cosine (RRC) filter for each of the ADC input stream. The RX RRC filter is the same as the TX RRC filter. Both RRC filters are used to filter out the out of band spurious signals and to improve the flatness of the in-band frequency domain response. Then de-chirp operations are applied using the baseband version of the TX chirp waveform. After de-chirping, the DC component of each input stream, which retains the TX-to-RX leakage signal, is removed. Then DBF is applied based on angle quantities a(0) and a(1). After DBF, 16 signal-streams from vector MTI filter form four partial sums in four quadrants, labeled as A, B, C, and D in Figure 17, for one target. These four signals are sent to four Chebyshev window function units to reduce the sidelobes of the main spurs.
FFT outputs are added together to form the sum signal after DBF. Note that the four FFT outputs are serialized from FFT index 0 to FFT index − 1 and form four time series at the baseband sampling rate. This is because -point time-domain I/Q signals fed in the FFT and bins are produced by the FFT engine. Since these FFT indices are proportional to the two-way propagation delays, these time series are called "post-FFT time series". The pre-FFT time series at the baseband sample rate is indexed by and the post-FFT time series at the baseband sample rate is indexed by . The key notations are summarized in Table 1.   Tone detection is performed within a specified range gate. The definition of the range gate will be discussed in a later section. The tone location in the post-FFT time series and the range value corresponding to the desired tone's FFT index are outputted. The error signals and for the measured target are formed based on the tone location and four FFT outputs of the four quadrants A, B, C, and D. Then the PI controller uses these error signals and the initial angle information to form the output angles quantities (0) and (1). Since the target moves continuously, the range gate needs to be adjusted continuously. This is done by adjusting the delay value of the range gate according to the estimated propagating time based on the range estimates. After the window function, the signals are fed into four FFT engines. Note that the time window of the FFT centers at t 1 and t 2 specified are shown in Figure 9. These four FFT outputs are added together to form the sum signal after DBF. Note that the four FFT outputs are serialized from FFT index 0 to FFT index N FFT − 1 and form four time series at the baseband sampling rate. This is because N FFT -point time-domain I/Q signals fed in the FFT and N FFT bins are produced by the FFT engine. Since these FFT indices are proportional to the two-way propagation delays, these time series are called "post-FFT time series". The pre-FFT time series at the baseband sample rate f B is indexed by k and the post-FFT time series at the baseband sample rate f B is indexed by p. The key notations are summarized in Table 1. Table 1. Notations of different indices used in the paper.

Index Name Definition m
The index of the receive antenna array along the x-axis. i The index of the receive antenna array along the z-axis. k The pre-Fast Fourier Transform (FFT) time index at the sample rate of f B . p The post-FFT time index at the sample rate of f B . n The time index for the 2T time duration intervals.

I g
The post-FFT time index of the range gate center.

I tone
The post-FFT time index of the desired tone location of the target.
Tone detection is performed within a specified range gate. The definition of the range gate will be discussed in a later section. The tone location in the post-FFT time series and the range value corresponding to the desired tone's FFT index are outputted. The error signals e 0 and e 1 for the measured target are formed based on the tone location and four FFT outputs of the four quadrants A, B, C, and D. Then the PI controller uses these error signals and the initial angle information to form the output angles quantities a(0) and a(1). Since the target moves continuously, the range gate needs to be adjusted continuously. This is done by adjusting the delay value of the range gate according to the estimated propagating time based on the range estimates.

De-Chirp of the Received Signal
The digital version of the transmitted chirp signal from the transmitter is written as x[k]. The output signal from de-chirp signal processing can be expressed as follows: where () * is the complex conjugate operation and r mi [k] is equivalently to r mi (kT B ), which is down-sampled from r mi (kT S ) by a down-sampling factor K = f S f B . As explained in Equation (16), a propagation delay in the round-trip path of the waveform is translated to a frequency shift from the center frequency. By analyzing the de-chirped signal r d,mk [k] in the baseband, we can see a distinct frequency component, away from the origin of the frequency axis. The amount of frequency shift is in proportion to the signal delay.

The Vector DC Canceller
It is well known that the useful signal of the target skin return can be buried by the TX-to-RX leakage signal at the receiver with the continuous wave. In some prior arts, the TX-to-RX leakage is attenuated by a DC canceller after the beamforming operation. This imposes some restrictions to the beamforming operation. For example, it requires a stepping response of the angle update for the beamforming operation and requires tight synchronization between the angle update and the DC cancellation. In our design, we use a per antenna DC canceller before the beamforming operation as illustrated in Figures 10 and 11. This is a first order IIR filter. In this manner, the strong DC of the TX-to-RX leakage is removed completely. We illustrate the filter response of the digital DC filter in Figure 12. is removed completely. We illustrate the filter response of the digital DC filter in Figure  12.
Although we can remove the DC component using the DC canceller, relatively weak near DC spurs are in the presence as illustrated later in Figure 18. They are caused due to the non-ideality of the RRC filter. The amplitude response of the RRC filter is not flat. Thus, the FMCW waveform amplitude is modulated by this amplitude response after passing the ideal FMCW waveform through the RRC filter.

The Vector MTI Filter
In order to remove the effects of the strong echoes of the clutters, an MTI filter is introduced. Since the antenna panel does not move with respect to the ground, the MTI filter can be applied in each receiving channel after each RX antenna before DBF. Normally, the MTI filter can be done in the post-FFT domain. However, it requires per-antenna FFT operations. As illustrated in Figure 10, the FFT operation is after the MTI filters, hence the time domain vector MTI filter is introduced in our design and illustrated in Figure 13. This is because the FMCW waveform is periodic with a period of 2 baseband samples. The MTI filter is FIR and has an order of 4. An IIR version can be introduced to further improve the filter response as a design choice. The MTI filter amplitude response versus the target radial speed is shown in Figure  14 as an example. The nulls of the filter correspond to the blind speeds of the target. The first null corresponding to the zero speed can remove the stationary or very low speed clutter echoes from the received signal. The second null corresponding to about 275.5 Although we can remove the DC component using the DC canceller, relatively weak near DC spurs are in the presence as illustrated later in Figure 18. They are caused due to the non-ideality of the RRC filter. The amplitude response of the RRC filter is not flat. Thus, the FMCW waveform amplitude is modulated by this amplitude response after passing the ideal FMCW waveform through the RRC filter.

The Vector MTI Filter
In order to remove the effects of the strong echoes of the clutters, an MTI filter is introduced. Since the antenna panel does not move with respect to the ground, the MTI filter can be applied in each receiving channel after each RX antenna before DBF. Normally, the MTI filter can be done in the post-FFT domain. However, it requires per-antenna FFT operations. As illustrated in Figure 10, the FFT operation is after the MTI filters, hence the time domain vector MTI filter is introduced in our design and illustrated in Figure 13. This is because the FMCW waveform is periodic with a period of 2T T B baseband samples. The MTI filter is FIR and has an order of 4. An IIR version can be introduced to further improve the filter response as a design choice.
filter can be applied in each receiving channel after each RX antenna before DBF. Normally, the MTI filter can be done in the post-FFT domain. However, it requires per-antenna FFT operations. As illustrated in Figure 10, the FFT operation is after the MTI filters, hence the time domain vector MTI filter is introduced in our design and illustrated in Figure 13. This is because the FMCW waveform is periodic with a period of baseband samples. The MTI filter is FIR and has an order of 4. An IIR version can be introduced to further improve the filter response as a design choice. The MTI filter amplitude response versus the target radial speed is shown in Figure  14 as an example. The nulls of the filter correspond to the blind speeds of the target. The first null corresponding to the zero speed can remove the stationary or very low speed clutter echoes from the received signal. The second null corresponding to about 275.5 km/h radial speed is called the blind speed. When the target speed is about an integer multiple of the blind speed, the MTI filter removes the target skin return as well. In the later part of this paper, we introduce the CFAR detector and range extrapolation as a work-around to this issue. The MTI filter amplitude response versus the target radial speed is shown in Figure 14 as an example. The nulls of the filter correspond to the blind speeds of the target. The first null corresponding to the zero speed can remove the stationary or very low speed clutter echoes from the received signal. The second null corresponding to about 275.5 km/h radial speed is called the blind speed. When the target speed is about an integer multiple of the blind speed, the MTI filter removes the target skin return as well. In the later part of this paper, we introduce the CFAR detector and range extrapolation as a work-around to this issue.

Range and Doppler Estimation
Let [ ] be the picked tone location from the sum signal in the FFT domain of [0, − 1] for the upbeat, we have: Similarly, for the downbeat, we can detect a tone at FFT index [ ]. We have: Based on Equation (26), we can estimate the range as: Given the baseband bandwidth , the normalized Doppler is defined as:

Range and Doppler Estimation
Let S 1 [n] be the picked tone location from the sum signal in the FFT domain of [0, N FFT − 1] for the n th upbeat, we have: Similarly, for the n th downbeat, we can detect a tone at FFT index S 2 [n]. We have: Based on Equation (26), we can estimate the range as: Given the baseband bandwidth f B , the normalized Doppler is defined as: Based on Equation (25), we can estimate the normalized Doppler value at the n th FFT window pair (the FFT window for the upbeat and the FFT window for the downbeat) as, Note that the estimated range and Doppler are updated at the 1/T sampling rate.

Range Gate
The range gate is denoted by I g − G, I g + G , where I g − G and I g + G lie within the FFT index range of [0, N FFT − 1], and I g is the center of the range cell under test in either upbeat or downbeat of the echo signal. The tone location I tone of the intended target needs to lie within the range gate. Note that the desired tone index I tone and range gate center index I g are not always the same. The quantity G allows an uncertainty from the range estimation process and accounts for the sidelobe region in the spectrum due to the window function.
Since the FMCW waveform is modulated, the post-FFT tone is strong in the spectrum main-lobe and it is very weak in the spectrum sidelobe (due to the window function), the conventional range gate control mechanism for the un-modulated signal [4] does not apply. In our design, the gate center control simply uses a moving average filter. The feedback part of the range gate tracking loop consists of "form gate", "smoothing", and "variable delay", as shown in Figure 10. The range gate is first formed using a square waveform generation circuitry with the properly generated gate width = 2G. Then the range gate tracking loop utilizes a variable delay function, which delays the range gate by the detected FFT tone location index in the post-FFT time series. To improve the tracking accuracy, the FFT tone location indices S 1 [n] and S 2 [n] are smoothed by smoothing R[n] and D dop [n]. This is the smoothing block illustrated in Figure 10. This moving average tone location is translated into a delay time and is used to delay the range gate accordingly.

The FMCW CFAR Detector
The standard cell-averaging CFAR detector [26] is applied to the post-FFT domain of the sum signal of DBF as illustrated in Figure 15. From [26], we adopt the CFAR detector figure here for illustration purposes. The cell under test is governed by the range gate of 2G + 1 post-FFT frequency indices centered at I g . In the range cell under test, we pick the index that corresponds to the maximum power of the 2G + 1 frequency indices, i.e., P tc max . Ideally, if the range gate is tracked properly, the center of the test cell corresponds to the maximum. Outside the range cell under test, we collect the average power of the 2G a post-FFT indices to determine the noise and residual clutter power floor, namely, P f loor . We define γ as the power ratio between P tc max and P f loor . When γ > γ th , the CFAR detector threshold, the index location with the maximum power will be recorded. Otherwise, the location of the maximum will not be discarded.
index that corresponds to the maximum power of the 2 + 1 frequency indices, i.e., . Ideally, if the range gate is tracked properly, the center of the test cell corresponds to the maximum. Outside the range cell under test, we collect the average power of the 2 post-FFT indices to determine the noise and residual clutter power floor, namely, . We define as the power ratio between and . When , the CFAR detector threshold, the index location with the maximum power will be recorded. Otherwise, the location of the maximum will not be discarded. If this tone location is recorded, the mono-pulse "Form error signals" in Figure  10 is activated at this tone location of the post-FFT time series. The mono-pulse error signals are then enabled and are fed into the PI controller. More on mono-pulse error signals will be given later. Otherwise, the mono-pulse angle tracking follows the angle memory stored in the PI controller in Figure 16 in the current time duration interval of .

Range Extrapolation Using Doppler Estimate for Low Signal-to-Noise Ratio (SNR)
When the strength of the signal tone is below the CFAR detection threshold, we use a scheme to extrapolate the range estimation based on the last known range estimate [ ] and last known Doppler estimate [ ]. Since the signal tone can be significantly attenuated when the radial velocity is close to the integer multiple of the MTI blind Figure 15. The example of cell-averaging constant false alarm rate (CFAR) detector (credited to [26]), where G = 4, G a = 12, and threshold = γ th 2G a . The tone location masks will be further classified as upbeat tone location masks and downbeat tone location masks as shown in Figure 20.
If this tone location I tone is recorded, the mono-pulse "Form error signals" in Figure 10 is activated at this tone location of the post-FFT time series. The mono-pulse error signals are then enabled and are fed into the PI controller. More on mono-pulse error signals will be given later. Otherwise, the mono-pulse angle tracking follows the angle memory stored in the PI controller in Figure 16 in the current time duration interval of T. speed, under this situation, the range extrapolation scheme helps to follow the actual range. The extrapolation uses the following simple formula: The range extrapolation can be effective only when the last known range and Doppler estimation are relatively accurate.

Four-Quadrant Digital Mono-Pulse Tracking
The four-quadrant mono-pulse tracking algorithm for pulse Doppler radar with mechanical steering can be found in [28]. Here, we extend this algorithm to the FMCW phased array angle tracking scenario. The antenna panel is equally divided into four quadrants, labeled as A, B, C, and D, as illustrated in Figure 17. The partial sum of received signals after digital beam forming in the quadrant = , , has the following form:

Angle Tracking Error Signals
The tone location in the post-FFT time series will be used to read the FFT outputs from the A, B, C, and D quadrants, denoted Then we can form the angle tracking error signals and as follows, Figure 16. PI controller. Input is the error signal and output is the estimate of angle tracking quantities.

Range Extrapolation Using Doppler Estimate for Low Signal-to-Noise Ratio (SNR)
When the strength of the signal tone is below the CFAR detection threshold, we use a scheme to extrapolate the range estimation based on the last known range estimate R[n last ] and last known Doppler estimate d dop [n last ]. Since the signal tone can be significantly attenuated when the radial velocity is close to the integer multiple of the MTI blind speed, under this situation, the range extrapolation scheme helps to follow the actual range. The extrapolation uses the following simple formula: The range extrapolation can be effective only when the last known range and Doppler estimation are relatively accurate.

Four-Quadrant Digital Mono-Pulse Tracking
The four-quadrant mono-pulse tracking algorithm for pulse Doppler radar with mechanical steering can be found in [28]. Here, we extend this algorithm to the FMCW phased array angle tracking scenario. The antenna panel is equally divided into four quadrants, labeled as A, B, C, and D, as illustrated in Figure 17. The partial sum of received signals after digital beam forming in the quadrant q = A, B, C or D has the following form:

Angle Tracking Error Signals
The tone location in the post-FFT time series will be used to read the FFT outputs from the A, B, C, and D quadrants, denoted Then we can form the angle tracking error signals and as follows, Figure 17. Four quadrants on the antenna plane, note that the view direction is from the receiver to the target, i.e., Y-direction points into the paper.
After the N FFT point FFT, the output from each quadrant is Y q [p], where p is the post-FFT time-series index. If target tone locations from CFAR detection are recorded in the upbeat and downbeat of the n th FFT window pair during the n th 2T time-period, the signal of these tones can be represented as Y u q [n] and Y d q [n], respectively. The sum signal is where ∠ stands for evaluating the angle of a complex number.

Angle Tracking Loop
The angle tracking loop relies on a closed loop control feedback mechanism. The classical PI controller is used to drive the digital beamformer. Starting from the initial a(0) and a(1) values in the proximity of true a(0) and a(1), the controller uses e 0 and e 1 to update the new estimated a(0) and a(1) values for the digital beamformer. The new estimated a(0) and a(1) will affect the estimation of the new e 0 and e 1 values.
The PI controller implementation in the baseband sampling domain has the structure shown in Figure 16.
The constant value c 0 = 0.0188. The constant c 1 = 1 256 . c 0 and c 1 are fine-tuned for a medium range scenario. They can be scaled by the inverses of the estimated range values.

Velocity Estimation
Based on two range and angle estimations separated by dd· 2T time duration, we can estimate the velocity on the x, y, z axes (38)

Results
In this section, we present the Simulink simulation results. The simulation setup is first presented. Then, in this controlled environment, we compare the angle, range, Doppler, and velocity estimation results of three configured scenarios with the true angle, range, Doppler, and velocity values.

Simulation Setup
We first simulate the effect of line-of-sight target skin return to lower the computational complexity in order to demonstrate a fully functional FMCW-based moving target tracking and measurement system. Then, we focus on the results in the scenario with a target, where it is fast moving and a strong clutter with zero velocity. Additionally, for simplicity, the simulation considers none of the effects of target glint, angle scintillation, atmospheric propagation, and multi-path propagation. The Simulink parameters are shown in Table 2.
The initial azimuth and elevation angle for the RX beam.

aa0 variable
The initial azimuth angle of the moving target.

taa0 variable
The initial value of the azimuth tracking angle.

ae0 variable
The initial elevation angle of the moving target.

tae0 variable
The initial value of the elevation tracking angle.

R rx,0 variable
The initial distance between the receiver and the target. Additionally, it is the range corresponding to the center of the initial range bin.

d TR 3 m
The distance between the transmitter and the receiver phase reference centers.

aa1 variable
The azimuth angle of the clutter.

ae1 variable
The elevation angle of the clutter.
The distance between the tracking and measurement system and the clutter.   rrcord 32 Root raised cosine pulse shaping filter order in the baseband sampling rate. This filter is applied to the transmitter. rrcord1 32 Root raised cosine pulse shaping filter order in the baseband sampling rate. This filter is applied to the receiver.  Note that the quantities like initial angles and initial ranges of the moving target and the clutter are variables changing based on the scenarios. This table shows the common parameters used for the simulation scenarios. We will define the three scenarios later. In each scenario, we apply the standard Chebyshev window with 80 dB sidelobe suppression.
With the 1024 point FFT, the resolution range of the radar in the frequency domain = 0.25·c·sweep time for one FFT index = 0.25·c·( f B /n f f t/2b) = 59.8 m. The range gate spans 9 frequency indices. The guard range of the range gate is translated from the guard frequency indices from the left or the right of the center frequency index in the range gate, and it is equal to 538 m. Due to the blind speed effect, the target tone is fully buried under the noise floor for a certain period of time, and even when the target moves out of the blind speed, the target tone is no longer at the center of the range gate. There exists a range offset between the times before and after the blind speed zone when the target is detected. Therefore, we need a wider range gate to accommodate this range offset so that the target tone can be re-detected when the target moves out of the blind speed.
The radar operates at 1.5 GHz carrier frequency. The triangular chirp signal has a period of 6528 T B samples. The blind speed is 275.5 km/h.
There are three simulation scenarios. Scenario I and II are used to demonstrate both the CFAR detector and the range extrapolation functions. Scenario III is used to demonstrate how the long range and high-speed tracking work.
The initial range and angle and the TX antenna gain information are common for the scenarios I and II. They are shown in Table 3. The initial range of the moving target is 20 km. The mean square of the velocity vector of the moving target is 3402 km/h, equivalently 2.74 Mach. The initial radial speed is about 1966 km/h, i.e., 7.14· blind speed. The target moves into the blind speed Doppler zone, then moves out of it in these scenarios. In the scenario I, we set the CFAR detector threshold = 10 dB and turn off the range extrapolation. In the scenario II, we set the CFAR detector threshold = 14.77 dB and turn off the range extrapolation. We also illustrate the matching of the SNR measurement and the SNR calculation in the simulations. We also illustrate the CFAR detector works with its detection diagram over time. Table 3. Simulation parameters for Scenario I and II.

Variable Name
Value Definition aa0 1.2490 radian The initial azimuth angle of the moving target.   R rx,0 20 km The initial distance between the receiver and the target. Additionally, it is the range corresponding to the center of the initial range bin.   Scenario III is a long-range and high-speed target scenario. The parameters like initial range and angle information are shown in Table 4. The initial range of the moving target is 40 km. The mean square of the velocity vector of the moving target is 7360 km/h, equivalently 5.93 Mach. The initial radial speed is about 3657 km/h, i.e., 13.28·blind speed. The target moves into the blind speed Doppler zone, then moves out of it in this scenario. Due to relatively low SNR in Scenario III, we just set the CFAR detector threshold = 10 dB and turn on the range extrapolation. ae0 π/2 radian The initial elevation angle of the moving target.
tae0 π/2 radian The initial value of the elevation tracking angle.
R rx,0 40 km The initial distance between the receiver and the target. Additionally, it is the range corresponding to the center of the initial range bin.   The TX antenna gain for the clutter assuming a uniform beam pattern on the azimuth and a very narrow beam pattern on the elevation. The gain on the azimuth plane is 21.46 dBi uniformly.

Unit Test for the Vector MTI Filter
The spectrum plot at 0.009 s of the simulation with the vector MTI filter turned off is shown in Figure 18. The simulation setup is using scenario II. The post beamforming noise floor in the simulation is normalized to 0 dBm (equivalent to 1 mW) with a gain normalization. This gain is introduced at the receiver and the signal power is scaled accordingly. Based on Table 2  The difference is small. The actual SNR of the target is 20.81 dB. Similarly, the signal power of the clutter after the FFT and the DBF is 38.95 dBm based on the calculation. The measured signal power of the clutter spur is 38.23 dBm. The difference is caused by the signal leakage that happens when the exact spur location in the frequency domain does not align with the FFT index grid.
The spectrum plot at the end of the 0.914 s simulation with the vector MTI filter turned on is shown in Figure 19. In this case, the noise floor is 0.97 dBm because the MTI filter enhances the noise power by 1.

Unit Test for the CFAR Detector
For scenario II, we plot the unit test results for the CFAR detector. In Figure 20, we plot the mask values of the tone locations, the mask values of the range gates, and the SNRs in dB scale. The mask takes a Boolean value of 1 or 0, representing whether or not the tone is detected in the second sub-figure or whether or not the range gate is selected in the first sub-figure. We observe that the CFAR detector cannot detect any tone location when the SNR is lower than 11 dB. Note that in this scenario, the CFAR detector threshold is 14.77 dB. This means that the cut-off of the successful detection does not happen at 14.77 dB, but lower than 14.77 dB. The main reason is that in CFAR detector, the cell averaging noise floor estimator uses only 24 cells, hence the noise floor estimation can be underestimated. In Figure 21, we zoom in the plot of the CFAR detection time diagram. The tone locations in upbeat and downbeat are at the centers of the range gates, respectively.

Unit Test for the CFAR Detector
For scenario II, we plot the unit test results for the CFAR detector. In Figure 20, we plot the mask values of the tone locations, the mask values of the range gates, and the SNRs in dB scale. The mask takes a Boolean value of 1 or 0, representing whether or not the tone is detected in the second sub-figure or whether or not the range gate is selected in the first sub-figure. We observe that the CFAR detector cannot detect any tone location when the SNR is lower than 11 dB. Note that in this scenario, the CFAR detector threshold is 14.77 dB. This means that the cut-off of the successful detection does not happen at 14.77 dB, but lower than 14.77 dB. The main reason is that in CFAR detector, the cell averaging noise floor estimator uses only 24 cells, hence the noise floor estimation can be underestimated. In Figure 21

Angle Estimation
In Figure 22, the estimated angle values in Scenario I and Scenario II are compared against the true angle values in terms of (0) and (1). Due to the blind speed effect, in Scenario I, there is about 0.4 s of time during which the intended tone has a low SNR, which is below or close to the CFAR detector threshold of 10 dB. The low CFAR detector threshold value causes a larger tracking error during this period of time than the result of Scenario II. In Scenario II, there is also about 0.4 s of time during which the intended tone has a low SNR below the CFAR detector threshold of 14.77 dB. In this case, the angle tracking is purely driven by the memory of the PI controller. When the radial speed is not close to an integer multiple of the blind speed, the tracking error is smaller than 0.01 radian (0.57 degree).
In Figure 23, the estimated angle values in Scenario III are compared against the true (0) and (1). During the 0.4 s of time when the SNR is below the CFAR detector threshold, the angle tracking is driven by the memory of the PI controller and occasional false alarms. When the target is out of the blind speed Doppler zone, the tracking error is smaller than 0.025 radian (1.43 degree).

Angle Estimation
In Figure 22, the estimated angle values in Scenario I and Scenario II are compared against the true angle values in terms of (0) and (1). Due to the blind speed effect, in Scenario I, there is about 0.4 s of time during which the intended tone has a low SNR, which is below or close to the CFAR detector threshold of 10 dB. The low CFAR detector threshold value causes a larger tracking error during this period of time than the result of Scenario II. In Scenario II, there is also about 0.4 s of time during which the intended tone has a low SNR below the CFAR detector threshold of 14.77 dB. In this case, the angle tracking is purely driven by the memory of the PI controller. When the radial speed is not close to an integer multiple of the blind speed, the tracking error is smaller than 0.01 radian (0.57 degree).
In Figure 23, the estimated angle values in Scenario III are compared against the true (0) and (1). During the 0.4 s of time when the SNR is below the CFAR detector threshold, the angle tracking is driven by the memory of the PI controller and occasional false alarms. When the target is out of the blind speed Doppler zone, the tracking error is smaller than 0.025 radian (1.43 degree).

Angle Estimation
In Figure 22, the estimated angle values in Scenario I and Scenario II are compared against the true angle values in terms of a(0) and a(1). Due to the blind speed effect, in Scenario I, there is about 0.4 s of time during which the intended tone has a low SNR, which is below or close to the CFAR detector threshold of 10 dB. The low CFAR detector threshold value causes a larger tracking error during this period of time than the result of Scenario II. In Scenario II, there is also about 0.4 s of time during which the intended tone has a low SNR below the CFAR detector threshold of 14.77 dB. In this case, the angle tracking is purely driven by the memory of the PI controller. When the radial speed is not close to an integer multiple of the blind speed, the tracking error is smaller than 0.01 radian (0.57 degree). Electronics 2021, 10, x FOR PEER REVIEW 24 of 29

Range Estimation
We plot the range estimate versus the true distance between the receiver (the origin of the coordinate system) and the target in Figures 24 and 25. In Figure 24, due to the higher CFAR detector threshold and the range extrapolation scheme in Scenario II, the range tracking error in Scenario II is much smaller than that of Scenario I during the "blind speed" period of about 0.4 s. During rest of the time, the tracking errors of the two scenarios are similar to each other. In Figure 25, the maximum range estimation error of Scenario III is lower than 200 m during the blind speed Doppler zone of 0.4 s time. In Figure 23, the estimated angle values in Scenario III are compared against the true a(0) and a(1). During the 0.4 s of time when the SNR is below the CFAR detector threshold, the angle tracking is driven by the memory of the PI controller and occasional false alarms. When the target is out of the blind speed Doppler zone, the tracking error is smaller than 0.025 radian (1.43 degree).

Range Estimation
We plot the range estimate versus the true distance between the receiver (the origin of the coordinate system) and the target in Figures 24 and 25. In Figure 24, due to the higher CFAR detector threshold and the range extrapolation scheme in Scenario II, the range tracking error in Scenario II is much smaller than that of Scenario I during the "blind speed" period of about 0.4 s. During rest of the time, the tracking errors of the two scenarios are similar to each other. In Figure 25, the maximum range estimation error of Scenario III is lower than 200 m during the blind speed Doppler zone of 0.4 s time.

Range Estimation
We plot the range estimate versus the true distance between the receiver (the origin of the coordinate system) and the target R rx in Figures 24 and 25. In Figure 24, due to the higher CFAR detector threshold and the range extrapolation scheme in Scenario II, the range tracking error in Scenario II is much smaller than that of Scenario I during the "blind speed" period of about 0.4 s. During rest of the time, the tracking errors of the two scenarios are similar to each other. In Figure 25

Doppler Estimation
The Doppler tracking results are presented in Figures 26 and 27. In Figure 26, during the "blind speed" period, the estimated Doppler value fluctuates more often in Scenario I than the estimated value in Scenario II due to the lower CFAR detector threshold value in Scenario I. The normalized Doppler tracking error is smaller than 0.002. The Doppler tracking error is thus smaller than 10000 Hz. In Figure 27, the Doppler tracking error of Scenario III is also lower than 10000 Hz.

Doppler Estimation
The Doppler tracking results are presented in Figures 26 and 27. In Figure 26, during the "blind speed" period, the estimated Doppler value fluctuates more often in Scenario I than the estimated value in Scenario II due to the lower CFAR detector threshold value in Scenario I. The normalized Doppler tracking error is smaller than 0.002. The Doppler tracking error is thus smaller than 10000 Hz. In Figure 27, the Doppler tracking error of Scenario III is also lower than 10000 Hz.

Doppler Estimation
The Doppler tracking results are presented in Figures 26 and 27. In Figure 26, during the "blind speed" period, the estimated Doppler value fluctuates more often in Scenario I than the estimated value in Scenario II due to the lower CFAR detector threshold value in Scenario I. The normalized Doppler tracking error is smaller than 0.002. The Doppler tracking error is thus smaller than 10, 000 Hz. In Figure 27, the Doppler tracking error of Scenario III is also lower than 10, 000 Hz. Electronics 2021, 10, x FOR PEER REVIEW 26 of 29

Velocity Estimation
The velocity tracking results are presented in Figures 28 and 29. The estimated velocities the in , , and axes during the "blind speed" period are not in agreement with the true velocities for all three scenarios, mainly due to the angle tracking error in the blind speed Doppler zone.

Velocity Estimation
The velocity tracking results are presented in Figures 28 and 29. The estimated velocities the in , , and axes during the "blind speed" period are not in agreement with the true velocities for all three scenarios, mainly due to the angle tracking error in the blind speed Doppler zone.

Velocity Estimation
The velocity tracking results are presented in Figures 28 and 29. The estimated velocities the in x, y, and z axes during the "blind speed" period are not in agreement with the true velocities for all three scenarios, mainly due to the angle tracking error in the blind speed Doppler zone. Electronics 2021, 10, x FOR PEER REVIEW 27 of 29

The RMS Errors When Normal Tracking Takes Place
From the angle, range, Doppler, and velocity plots, both scenario I and II have the common normal tracking period beyond 0.5 s of simulation time. Normal tracking simply means that the SNR values in each scenario after the MTI filter are above each of the CFAR detector thresholds of each scenario. Since all random generators for noise use the same random seeds, the simulation results are roughly the same after 0.5 s of simulation time. The calculated RMS errors are shown in Table 5.

The RMS Errors When Normal Tracking Takes Place
From the angle, range, Doppler, and velocity plots, both scenario I and II have the common normal tracking period beyond 0.5 s of simulation time. Normal tracking simply means that the SNR values in each scenario after the MTI filter are above each of the CFAR detector thresholds of each scenario. Since all random generators for noise use the same random seeds, the simulation results are roughly the same after 0.5 s of simulation time. The calculated RMS errors are shown in Table 5.

The RMS Errors When Normal Tracking Takes Place
From the angle, range, Doppler, and velocity plots, both scenario I and II have the common normal tracking period beyond 0.5 s of simulation time. Normal tracking simply means that the SNR values in each scenario after the MTI filter are above each of the CFAR detector thresholds of each scenario. Since all random generators for noise use the same random seeds, the simulation results are roughly the same after 0.5 s of simulation time. The calculated RMS errors are shown in Table 5. The total recorded simulation time is roughly 5000 s for 0.914 s of real time simulation. To speed up the simulation and the processing, the proposed hybrid time and frequency domain signal processing with target moving scenario will be planted on a commercial FPGA board, so that much faster M&S can be achieved.

Discussion
The simulation results presented in Section 5 agree with the theoretical predictions given by different CFAR detector threshold values. On the other hand, the hybrid timedomain and frequency-domain processing design has lower complexity than the per antenna FFT design. The number of FFT engines used is 4 instead of 16 for the case of 16 antenna. In future work, we would like to build a moving target tracking and measurement test bed with the presented receiver processing unit to collect some experimental data. This testbed will be set up when the FPGA implementation, the RF frontend, and the antenna system are ready.

Conclusions
The system model of a digital beamforming phased array system for fast moving object skin-return signal processing is presented in this paper with numerical modeling in Simulink. Using FMCW modulated EM wave, the proposed array system interrogates a moving target in a clutter environment and measures object in range, angle, and Doppler using joint time and frequency domain signal processing. Comparing with existing FMCW based tracking systems, the new system has much lower complexity, especially for a DBF array with large number of elements. Even at 40 km range and 5.93 Mach radialspeed, simulation results show that the system root-mean-square errors of range, angle, and Doppler tracking are 26 m, 0.68 degree, and 1100 Hz, respectively. To speed up the simulation and the design process, the proposed hybrid time and frequency domain signal processing with target moving scenario will be planted on a commercial FPGA board. This way, much faster M&S can be achieved. In addition, the actual baseband processing chain on FPGA would allow real measurements in field experiments.