Theoretical background

A common approach for isolating reflections leading to target detections are constant false alarm rate algorithms. They perform well against background noise and interference [Reference Gamba13]. However, the approach of this paper is the recreation of the sensor receiver values including effects which may challenge the target detection of a radar. In order to form a digital twin of the radar data, the information is filtered using CLEAN. The goal is an approximation of the sensor receiver values with a limited set of point targets. CLEAN assumes the data can be described by a superposition of point sources which are affected by a point spread function (PSF). By decomposing the original data into a set of point sources and removing the PSF, a list of reflections is created which can be used as target list. Figure 2 describes the algorithm. The terminology assumes a “Dirty Map” or “Dirty Image” as initial dataset to be filtered. The Dirty Map is scanned for its maximum located in bin I. The peak multiplied by the PSF is considered the “Dirty Beam.” The Dirty Beam is then subtracted and the residual becomes the “Residual Map” to be scanned for the next iteration. The algorithm continues until a defined threshold is reached.

Figure 2.

CLEAN algorithm flow chart.

Figure 2 Long description

Begins with ′Start′ leading to ′Initiate Dirty Map′. The next step is ′Identify peak in bin I′, followed by a decision point ′Stop criteria reached?′. If ′No′, it proceeds to ′Form Dirty Beam around Peak′, then ′Subtract Dirty Beam from Residual Map′ and ′Update Residual Map, Update Clean Map(I)′. This loops back to ′Identify peak in bin I′. If ′Yes′, the process ends with ′End′.

Each RD map $R$ can be sampled by a sum of weighted and shifted PSFs $A\cdot PSF(x-I)$ according to:

(1)\begin{equation} R = \sum_{i=1}^k\kappa \cdot A_i \cdot PSF(x-I_i) \end{equation}

$\kappa$ is the loop gain and can be chosen to (0,1]. With a loop gain $ \lt 1$, a reflection can be sampled by multiple smaller ones resulting in a better recreation of the signal. Especially sampling reflections from extended objects could be improved by this. As the number of available point reflections scales with hardware in an RTS, it is beneficial to synthesize the environment with the lowest amount of point targets. Therefore, a loop gain of $\kappa = 1$ is chosen. $A$ scales the PSF with the amplitude of the measured peak.

Field data collection

Radar field data are recorded using an automotive radar evaluation kit TI AWR1642 with DCA1000 evaluation board mounted at the front of a car. While driving, analog-to-digital converter (ADC) values of the intermediate frequency signal are recorded. The collected data consist of radar cubes with 272 ADC samples, 255 chirps, and 4 receive channels. Every 500 ms, a new radar cube is fetched as radar snapshot. A total of 600 snapshots are recorded, resulting in a scenario duration of 5 minutes. For further analysis, a second dataset of the same environment is recorded with 5999 snapshots in the same time using the radar evaluation kit TI AWR1843. The datasets will be referred to as dataset 1 (AWR1642) and 2 (AWR1843). Parameters regarding the sensor configurations are summarized in Table 1.

Table 1.

Radar parameters of dataset 1 and 2

Table 1 Long description

Radar parameters are listed for two devices, with columns labeled Value, AWR1642 and Value, AWR1843. Start Frequency is 77 GHz for AWR1642 and 76 GHz for AWR1843. Range Resolution is 0.1526 m for AWR1642 and 0.2891 m for AWR1843. Doppler Resolution is 0.0919 m per second for AWR1642 and 0.0885 m per second for AWR1843. Radar Cube is 272 by 255 by 4 for AWR1642 and 108 by 128 by 4 for AWR1843. Stored Sample Width is 16 bits for both, and Sampling Rate is 5.535 MSPS for AWR1642 and 12.5 MSPS for AWR1843. RX Channels are 4 for both, and Field of View is plus or minus 60 degrees for both.

As driving environment, a parking garage was chosen. The scenario is common in the automotive domain and features many large close-by scatterers like cars and metallic bars. Figure 3 shows a picture of the environment. The receiver noise floor is estimated by a background measurement against sky to avoid the measurement of any reflection. The mean and standard deviation $s$ of the background data are calculated. Based on this, 3 $s$ environment around the mean noise level is excluded from the scenario data to ensure that noise will not contribute to the measurements. Furthermore, the field of view of the sensor is limited to $\pm$60 $^{\circ}$. Data with wider estimated AoA are excluded.

Figure 3.

Parking garage as measurement environment.

ADC samples processing

The radar cube of each snapshot can be processed by applying an FFT. Its input data are multiplied by a Blackman window, zero-padding of factor 4 was applied. Reflections can be identified in the RD map. For AoA estimation, the Bartlett beamformer is used [Reference Gamba13]. As output, receiver amplitudes normalized to the maximum and given in dBFS (decibels relative to maximum level or full scale) are available along range, Doppler, and AoA. Receiver values of the scenario are further regarded for the dataset 2 from TI AWR1843. Histograms of all values above the noise floor are summarized in Fig. 4 showing characteristic details about the environment and the applied processing. The vehicle was manually driven with a targeted speed around 10 km/h agreeing with 2.78 m/s. Due to manual driving and accelerating, the targeted speed can be slightly exceeded. The range from 6 to 20 m was further analyzed. An RD heat map shows maxima occurring in the scenario. It shows most reflections at the speed of the ego-vehicle and features angle-dependent Doppler migration effects from the detection of static objects and curve driving. The histogram of recorded velocities shows its peak at the measured negative ego-vehicle speed of $-$3.281 m/s and rolls off toward 0 m/s. The detected AoA histogram shows a maximum at 0 $^{\circ}$, which can be explained by the architecture of a garage forming a tunnel and reflecting energy back preferably at center axis. Obstacles along the driving way are present at all angles but appear dominantly as parking cars to the left and right. In addition, the sensor antenna pattern decreases toward larger angles resulting in fewer received values above the noise floor.

Figure 4.

Histograms of recorded scenario.

Figure 4 Long description

The A showing a heat map with the horizontal axis labeled, Velocity, bracket m slash s bracket, ranging from minus 4 to 4. The vertical axis is labeled, Range, bracket m bracket, ranging from 6 to 18. A vertical scale bar is labeled, Amplitude, bracket dBFS bracket, with tick labels minus 55, minus 65 and minus 75. The B showing a histogram with the horizontal axis labeled, Angle, bracket degree bracket, ranging from minus 60 to 60. The vertical axis is labeled, Counts times 10 superscript 6, ranging from 0 to 1. Bars form a distribution with a tall narrow peak near 0 and broader groups of bars around negative 30 to negative 10 and around 20 to 40. The C showing a histogram with the horizontal axis labeled, Range, bracket m bracket, ranging from 6 to 18. The vertical axis is labeled, Counts times 10 superscript 6, ranging from 0 to 1. Bars are highest near 8 to 10 and then decrease toward 18. The D showing a histogram with the horizontal axis labeled, Velocity, bracket m slash s bracket, ranging from minus 4 to 2. The vertical axis is labeled, Counts times 10 superscript 6, ranging from 0 to 4. Bars peak between about minus 3 and minus 2, then decrease toward 2.

For simplification, the 2D RD map is compressed to a 1D range plot with the maximum receiver amplitudes across the Doppler dimension. Doppler shift is introduced by the only moving object of the scenario, which is the ego-vehicle itself. Furthermore, due to the compression and low antenna gain in the elevation plane, an AoA in the azimuth plane is assumed and each range bin in the 1D plot is associated with a single AoA. These simplified processing steps are expected to affect a target list and clutter statistics as, e.g., weak scatterers will be masked by stronger ones in a 1D plot. For converting the received values into RCS, a Keysight E8719A RTS is used. The received amplitude from the sensor is corrected against a reference target from the RTS. Any type of RTS or reference target with a fixed known RCS could be used. The gain pattern of the radar antenna is known; therefore, each AoA can be mapped to an antenna gain. The setup is shown in Fig. 5. By applying the radar equation [Reference Gamba13], an RCS and the according received power for each range bin can be calculated.

Figure 5.

RCS correction factor measurement using an RTS.

Radar target emulation

An RTS modifies the radar signal emulating a fixed amount of point reflections with different delays, Doppler shifts, and amplitudes. This results in reported point targets with range, relative velocity and RCS [14]. Some architectures also allow the emulation of AoA [Reference Dallmann, Mende and Wald15, 16]. Advanced RTS solutions feature an arrangement of antennas similar to a screen allowing a multitude of point targets with multiple AoA, which can be utilized for playback of detection lists [Reference Holzinger, Liebau and Metzger17]. A final detection list always needs to be adapted to possible hardware limitations of the RTS, like, e.g., a maximum number of simultaneous targets [Reference Leon, González-Huici and Dallmann3]. For creating a digital twin from the sensor output data, it has to be mentioned that there is no standardized sensor output and recent research is discussing a possible metric to establish the objective validation of point cloud outputs [Reference Schmidt, Widera, Kahlert and Bertram18]. The chosen approach aims to extract range, Doppler, RCS, and AoA information from the sensor for all received frames as these are parameters which can be set by an RTS and approximates sensor information by providing a set of scaled point reflections. An ideal point target has a constant RCS concentrated in a single RD map bin at a single AoA. This can be approximated with a corner reflector or point target being emulated by an RTS [Reference Gowdu, Asghar, Nagel, Rozmann, Stephan and Hein19]. In the mathematical description of a propagation channel, an equivalent representation of a point target is the unit-impulse $\delta$ [Reference Sundararajan20]. The processed output of a radar reporting a point target represents its system response. It includes radar-specific artifacts due to, e.g., antenna, windowing, and additional processing steps and can be used as PSF. A PSF for both datasets was measured, Fig. 6 shows the measured PSF of the TI AWR1642 in normalized power. Windowing (Blackman) was applied to reduce FFT sidelobes below the receiver noise floor. A PSF should be measured as high as possible above the noise floor ensuring scalability to various point target RCS sizes. Based on the radar configuration profile and the processing of ADC values, the PSF will change.

Figure 6.

Measured and normalized point spread function with $-$3 dB contour.

Target list extraction

After FFT processing, the RD map of each radar snapshot is available and compressed to a 1D range plot as described in “ADC samples processing” section. Its dBFS amplitude values are shown in Fig. 7. Peaks correspond to range bins and can include reflections from multiple AoA. The information can be converted in a range vs. RCS plot. Figure 8 shows RCS values according to the antenna gain at the estimated AoA. As only one dominant AoA is estimated per range bin, the curve may have discontinuities. It is used as Dirty Map for searching peaks according to Fig. 2. As iteration stop criterion for the parking garage scenario, an RCS level of 0 dBsm is set. An RCS of 0 dBsm is specified as typical value for a bicyclist as smaller traffic participant on the road [14]. The iteration stop criterion depends on the overall Signal-to-Noise Ratio of the analyzed scenario and is limited by the available dynamic range of the sensor. The iteration stop should be below the average scenario RCS for analyzing clutter statistics. For a sensor with higher dynamic range, a lower limit of $-$8 dBsm is proposed to include typical RCS of pedestrians [14] and thus to include all vulnerable road users. Table 2 summarizes the parameters used in the CLEAN algorithm, as well as those for radar processing. For creating a target list for an RTS with limited point targets per snapshot, a fixed number of targets can be set as stop criterion as well.

Figure 7.

Snapshot range vs. amplitude in dBFS.

Figure 8.

Snapshot range vs. RCS.

Table 2.

CLEAN algorithm and radar processing parameters

Table 2 Long description

Parameters are listed for Dataset 1 and 2. Loop Gain, kappa, is 1. Iteration Stop Criterion is 0 dBsm Residual Map Level. Zero Padding Factor is 4. Window is Blackman. PSF Measurement Method is RTS E8719A point target in RD map.

All snapshots are analyzed, a list of possible targets and the residual are recorded. Figure 9 shows the resulting superposition of point targets in RCS, up to the iteration stop of 0 dBsm, for a single range profile. Figure 10 shows the superposition of corresponding extracted PSF with respect to the iteration stop expressed in the originally received dBFS values. The iteration stop at 0 dBsm is converted into the corresponding dBFS values based on the estimated AoA. However, the list of targets for the digital twin needs to have a fixed AoA for each PSF. For the first dataset, CLEAN operates on reported magnitudes without phase information. A superposition of closely spaced overlapping targets is therefore considered with the same phase. The output of CLEAN shall be used as input for an RTS architecture. However, RTS architectures are not always able to offer the full control of the output phase, and an emulation of targets by RTS is often assumed to happen with random phases. As extension to the previous paper [Reference Sims, Holzinger, Vandeplas, Dallmann and Heberling12], the algorithm operates on reported power for the second dataset and assumes the incoherent superposition of radar targets. A full phase-coherent superposition is expected to yield the most realistic superposition of overlapping targets but requires phase information from the radar and phase control of the RTS architecture, which led to the decision of relying on an incoherent CLEAN implementation.

Figure 9.

Snapshot range vs. RCS as superposition of point targets.

Figure 10.

Snapshot range vs. amplitude recreated by weighted PSF.

Residual analysis after CLEAN

All values below the iteration stop are considered clutter. These values are intended to be added as random fluctuating values in target emulation. The residual of CLEAN is defined by the iteration stop criterion. In this paper, the stop criterion is defined with a fixed RCS value. Thus, the residual contains measured RCS values between the sensor detection threshold and the stop criterion of 0 dBsm. Stronger scattering centers are removed by CLEAN and the residual only shows low fluctuating echoes, which are expected to also originate from multipath. Its values are measured at different time snapshots and different ranges. To correctly interpret the residual statistics, both processing effects and sensor characteristics must be considered:

1) Range- and angle-dependent detectability:

Due to free-space path loss, targets with identical RCS are measured with different receive powers depending on their detected range. This results in a minimum detectable RCS that increases with target distance. Range-dependent detection effects are expected to bias the observed statistics and also affect the true RCS mean values of the environment. With a CLEAN iteration stop at 0 dBsm, the dynamic range of detectable clutter decreases with distance. At larger distances, fewer clutter samples can be observed, and only stronger residual contributions remain, as weaker scatterers fall below the minimum detectable RCS. Therefore, clutter statistics will also vary depending on analyzed range. An angle-dependent antenna gain further affects the minimum RCS by range to be detected at the boresight of the antenna only.

2) PSF-related and signal-processing effects:

Clutter contains the unfiltered residual content of the RD map and thus includes the target RCS roll-off from the PSF. It is therefore influenced by the scenario dynamics in range and time but also contains specific characteristics of the sensor itself. In addition, mathematical artifacts and imperfections of CLEAN can also contribute. The processing, FFT zero-padding, and compression of data influence statistics. Especially the simplified 1D compression can mask weak scattering centers, and a full 2D RCS estimation is expected to reveal a better model of the clutter statistics.

The environment can be considered as homogeneous due to a repeating driving path inside a parking garage. Due to the moving ego-vehicle, recorded reflections change their position each time snapshot. In order to analyze the fluctuation of clutter, a probability density function (PDF) of all received residual RCS can be calculated. The histogram in Fig. 11 shows the residual of the first dataset after CLEAN in RCS. The dataset was recalculated in this paper with improved angle estimation. By applying curve fitting, an exponential PDF can be associated with the histogram. The residual RCS values $\sigma$ are showing an exponential behavior which resembles the fluctuating RCS of a Swerling I target [Reference Gamba13]. This target model corresponds to a multitude of similarly large isotropic reflectors. From snapshot to snapshot, the composition and/or the position of the reflectors changes and so does the reported RCS. While the classical Swerling I model describes the snapshot-to-snapshot variation of an individual target, the observed behavior here results from an ensemble of weak, unresolved scattering contributions whose composition changes due to ego-vehicle movement. The exponential PDF therefore characterizes the aggregate fluctuation behavior of the residual clutter rather than that of a single physical target. In addition, the processed data include interpolated samples obtained via zero-padded FFT processing, which introduces correlations between adjacent bins and should therefore be regarded as oversampled representations of the same measurement rather than independent observations. This is expected to influence statistics and the true temporal scene mean. In Swerling I, the RCS is constant in a single snapshot but varies from snapshot to snapshot according to:

(2)\begin{equation} \ P(\sigma)= \frac{1}{\sigma_{avg}}\exp\left(- \frac{\sigma}{\sigma_{avg}}\right),\quad \sigma \geq 0 \end{equation}

with a mean RCS of $\sigma_{avg}$. A dataset mean RCS is calculated from all received residual clutter samples. In addition to that, extracted target lists from the scenario can also be analyzed for statistics in AoA, RCS, Doppler, and range for model creation of characteristic driving environments.

Figure 11.

Histogram of residual (Dataset 1).

As CLEAN operates without phase information, an incoherent superposition of overlapping targets is further assumed and CLEAN is performed again on both datasets. The histogram in Fig. 12 shows the residual of both datasets using incoherent subtraction in CLEAN. New effects in the PDF due to the changing physical interpretation of overlapping targets can be detected. Compared to Fig. 11, a larger mean residual RCS is observed. This behavior is expected, as incoherent subtraction of overlapping targets yields larger residual values than a subtraction in phase. The exponential fit exhibits a reduced matching, particularly at low RCS values, where the expected PDF peak is less pronounced. Both datasets show similar histogram shapes and mainly deviate at the exponential peak, which could be explained with improved and more sensitive measurement methods and slightly different sensor characteristics and antenna gains. The second dataset is therefore used for further analysis.

Figure 12.

Histogram of residual for incoherent processing (Dataset 1 vs. Dataset 2).

Figure 12 Long description

The vertical axis label is, Probability Density. The vertical axis range is 0 to 4 with tick labels at 0, 1, 2, 3 and 4. The horizontal axis label is, RCS (linear). The horizontal axis range is 0 to 1 with tick labels at 0, 0.2, 0.4, 0.6, 0.8 and 1. A legend lists six entries: Histogram Data Set 2, Fit Data Set 2, Swerling I Data Set 2, Histogram, Fit, Swerling I. Two overlapping histogram bar series span from near 0 on the horizontal axis to near 1. The bar heights are highest near the left side and decrease toward the right side. Two fitted curves are drawn over the histograms. One is a solid line and one is a dashed line. Both start near a probability density around the mid 2 range close to RCS 0 and decrease toward a probability density near 0 by RCS 1. Two marker series are plotted over the bars and curves. One uses asterisk shaped markers and one uses x shaped markers. The markers are densest across the full horizontal range from near 0 to near 1 and follow the same decreasing pattern from higher values near the left to lower values near the right.

Limitations of residual analysis

Effects of a limited minimum detectable RCS mentioned in “Residual analysis after CLEAN” section will be analyzed in the following. By applying the radar equation, a minimum detectable RCS and the available dynamic range of clutter can be estimated. The clutter dynamics are limited by the RCS detection threshold and the iteration stop of 0 dBsm. Although the environment is considered to be homogeneous, this influences its radar perception. The clutter was separated into groups of smaller range patches which feature a change of the minimum detectable RCS with rising range of up to 3 dB. In addition, the angular field of view was limited to 3 dB beamwidth of the combined TX/RX antenna of the sensor. Thus, effects from range- and angle-dependent detectability were reduced, while still a large number of samples is available. The aggregated PDF of measured RCS across time snapshots can be calculated. Depending on range, a different amount of samples is available affecting the PDF. Especially at short ranges, clutter exhibits a larger dynamic range and a high number of measured samples. Figure 13 shows the corresponding PDFs with exponential fits, whose parameters were estimated by fitting the exponential function via nonlinear least-squares using the “fit” function of MATLAB (“Trust-Region” algorithm) with the initial value set to the samples mean. However, a matching of the fitted curve to the Swerling I model does not guarantee an exponential nature of the PDF, the deviation from the histogram can be observed. With rising range cells, migration of the PDF maximum is visible. This can be explained by the rising minimum detectable RCS imposed by the sensor noise floor rather than by the physical environment. Accordingly, the agreement with a Swerling I model evolves from an overemphasized peak at short ranges to a good match at intermediate ranges and degrades clearly again at larger ranges. Interpolated samples are further expected to distort especially the PDF peaks. In addition, larger ranges also show fewer clutter detections; thus, fewer samples and are less suitable for statistical analysis. As example, the first PDF consists of 39,704 samples and features a dynamic range from 0 dBsm down to $-$24.55 dBsm for 6–7.15 m, the last PDF consists of 31,133 samples down to $-$5.74 dBsm for 14.88–17.7 m. In theory, Swerling I describes the fluctuation of a single target at a single range bin. But as scattering centers keep moving, even the analysis of a single range bin contains echoes from multiple migrating targets. For the analysis of the scenario clutter, it therefore makes no difference if adjacent range bins are pooled and thus form a global PDF of multiple equally scattering centers. But in addition, a sufficient and ideally equal dynamic range for clutter samples should be ensured for a statistical analysis. Thus, while the exponential behavior can resemble a Swerling I fluctuation model, the resulting PDF is shaped by the detection capability of the sensor rather than representing the properties of the of the environment alone. If the full scene over all ranges is pooled as a global PDF, the limited sample amount at larger distances contributes less to the statistics.

Figure 13.

Histograms of multiple range patches.

Figure 13 Long description

Six histogram plots arranged in three rows and two columns, titled, Cells 1 to 17, Cells 18 to 37, Cells 38 to 60, Cells 61 to 89, Cells 90 to 123 and Cells 124 to 163. The vertical axis label on the left plots is, Probability Density, with the abbreviation in parentheses. The horizontal axis label on the bottom row plots is, RCS, with the word in parentheses. Cells 1 to 17: The horizontal axis shows tick labels 0, 0.5 and 1. The vertical axis shows tick labels 0, 5 and 10. A histogram has the tallest bars near 0 and then decreases toward 1. A smooth curve overlays the bars and decreases from near the top left toward near 0 at the right. Star markers appear along the curve. Cells 18 to 37: The horizontal axis shows tick labels 0, 0.5 and 1. The vertical axis shows tick labels 0, 2 and 4. A histogram has higher bars near 0 and then decreases toward 1. A smooth curve overlays the bars and decreases toward the right. Star markers appear along the curve. Cells 38 to 60: The horizontal axis shows tick labels 0, 0.5 and 1. The vertical axis shows tick labels 0, 1, 2, 3 and 4. A histogram has higher bars near 0 and then decreases toward 1. A smooth curve overlays the bars and decreases toward the right. Star markers appear along the curve. Cells 61 to 89: The horizontal axis shows tick labels 0, 0.5 and 1. The vertical axis shows tick labels 0, 1, 2 and 3. A histogram has higher bars near 0 and then decreases toward 1. A smooth curve overlays the bars and decreases toward the right. Star markers appear along the curve. Cells 90 to 123: The horizontal axis shows tick labels 0, 0.5 and 1. The vertical axis shows tick labels 0, 1, 2 and 3. A histogram begins with lower bars near 0, rises to a peak at a small x value, then decreases toward 1. A smooth curve overlays the bars and decreases toward the right. Star markers appear along the curve. Cells 124 to 163: The horizontal axis shows tick labels 0, 0.5 and 1. The vertical axis shows tick labels 0, 1 and 2. A histogram begins with low bars near 0, rises to a peak around the left to middle of the axis, then decreases toward 1. A smooth curve overlays the bars and decreases toward the right. Star markers appear along the curve. A legend at the bottom lists three items: Clutter PDF, Fit and Swerling I.

The findings show that sensor characteristics limit the perception of a cluttered scenario. Applying an exponential fluctuating model as Swerling I to resemble clutter appears to be a reasonable global model to emulate a more realistic clutter to target ratio. The sensor characteristics will shape the received clutter values for larger ranges by not displaying low exponential fluctuating values. However, the analysis also shows that sensor characteristics need to be considered while formulating clutter models to be used for emulation. It is recommended to consider this by analyzing smaller range patches with reduced angular coverage of similar antenna gain for local fluctuation models.