Iterative 2D sparse signal reconstruction with masked residual updates for automotive radar interference mitigation

Compressive sensing has attracted considerable attention in automotive radar interference mitigation. However, these algorithms usually cannot be applied directly to commercial automotive radar as most of them are computationally intense. In this paper, we therefore introduce a computationally efficient two-dimensional masked residual updates (2D MRU) compressive sensing framework. By utilizing the sparsity of the beat signal in the frequency domain, the range-Doppler (RD) spectrum can be reconstructed with the help of undistorted samples in the beat signal. Unlike the other schemes, where a 2D signal measurement is vectorized into a 1D signal, the proposed 2D MRU can directly take a 2D signal measurement and reconstruct the corresponding RD spectrum. Furthermore, the 2D MRU framework can be easily integrated into well-known optimization schemes such as basis pursuit, iterative hard thresholding, iterative soft thresholding, orthogonal matching pursuit, and approximate message-passing algorithm. In addition to the standard iterative thresholding algorithms, we propose a novel prior-model-based iterative thresholding method to further reduce the computation time and reconstruction error. Theoretical analysis shows that the proposed framework can successfully reconstruct the RD spectrum with high probability. Moreover, numerical experiments demonstrate the superiority of the proposed framework in terms of computational complexity.


Prior work
Different signal reconstruction algorithms, such as the zeroing method, autoregressive (AR) model-based interpolation and adaptive noise cancellation have been proposed as countermeasures for automotive radar interference mitigation [1,[4][5][6]. The zeroing method [1] is a basic technique in which the signal segment disturbed by interference is simply set to zero. In [4], a complex band receiver and an adaptive noise canceller are utilized for interference mitigation, it cancels the interference in the positive half of the frequency spectrum ( [0, π ] normalized frequency), while it uses the correlated interference in the negative half of the frequency spectrum ( [−π , 0] normalized frequency) as a reference. Solving the mutual interference issue in the frequency domain for frequency-modulated continuous wave (FMCW) radar is also considered in [5], where the signal disturbed by interference is recovered by linear predictive coding (LPC) in the short-time Fourier transform (STFT) domain. However, the interference in real-world scenarios can contain much diversity and randomness, e.g., the chirp rates of different radar sensors can be similar to each other, and there can be multiple sources of interference. Hence, the samples disturbed by interference in some chirps may spread across the entire measurement. Thus, the number of available interference-free samples in these chirps is not sufficient for LPC to provide adequate recovery in the STFT domain. Different to [5], the AR model is used in [6] for the discrete beat signal interpolation in the positions which are disturbed by interference. Recent research results have shown that compressive sensing (CS) approaches can be used to tackle this problem [7][8][9][10]. In [7], the samples disturbed by interference are first substituted with zeros and these zeroed samples are further interpolated using the iterative method with an adaptive thresholding (IMAT) algorithm. In [8], the interference is removed by using the orthogonal matching pursuit (OMP) algorithm to project the interference-contaminated signals on a reduced chirplet basis which contains possible hypotheses for slopes and time-shifts of chirps. The interference mitigation problem is defined as a dual-basis pursuit problem and the morphological component analysis is used in [9] for separating the interference from the discrete beat signal. In [10], the sparse Bayesian learning algorithm is adopted for automotive radar interference mitigation. The range-Doppler (RD) spectrum can be acquired from the mean of the maximum a posteriori (MAP) estimate under the given remaining undistorted samples in the discrete beat signal. Most recently, the wavelet denoising technique was proposed to suppress the mutual interference in the time domain [11]. However, to apply this method in practical automotive radar systems, the threshold values related to the wavelet coefficients need to be optimized on the basis of large amounts of real world data and the interference with similar amplitude as the echo signal of the target may still remain after reconstruction.

Motivation and contribution
As exemplified above, CS algorithms usually provide better signal recovery than conventional algorithms but might require more computational resources. Therefore, we propose a novel solution for the interference mitigation that can efficiently reconstruct the RD spectrum by directly processing the two-dimensional (2D) discrete beat signal. The proposed 2D masked residual updates (2D MRU) CS framework can be easily integrated into most existing CS schemes as well as the proposed novel prior-model-based iterative thresholding methods. Specifically, the main contributions of this paper are as follows: Firstly, we propose a novel prior-model-based iterative thresholding method that achieves smaller reconstruction error and computation time than corresponding baseline algorithms. The update rate (typically 20 Hz) [12] in automotive radar systems indicates that rich prior information can be provided for interference mitigation, since the perceived environment may change only slightly during the short measurement cycles.
Secondly, a computationally efficient 2D MRU CS framework for RD spectrum recovery is introduced. 2D MRU avoids large vectors and high-dimensional matrix operations which typically come along with an increasing size of the measurement matrix.
Thirdly, the measurement transform matrix used in this paper is proven to satisfy the restricted isometry property (RIP) with high probability. In other words, the proposed framework can successfully reconstruct the RD spectrum with a high probability.
Fourthly, it is demonstrated that the proposed 2D MRU can be easily integrated to the established CS schemes e.g., basis pursuit, iterative soft thresholding (IST), iterative hard thresholding (IHT), OMP, approximate message-passing (AMP) and the proposed priormodel-based iterative thresholding algorithms. The performance of different algorithms that incorporate 2D MRU is evaluated and the superiority of the proposed framework in terms of computational complexity is shown.
Finally, a detection method with a higher true-positive rate is proposed for the detection of interference-contaminated samples. The detection method can be further combined with the state-of-the-art CS frameworks as well as with 2D MRU.

Organization and notations
The organization of this paper is as follows. In Sect. 2, the frequency-modulated continuous wave radar signal model is introduced. Section 3 describes the details of the proposed prior-model-based iterative thresholding algorithms and 2D framework and we prove that the proposed framework can successfully reconstruct the RD spectrum with high probability. In Sect. 4, an algorithm for the detection of interference-contaminated samples is proposed. The performance of the proposed method is evaluated through real measurements in Sect. 5. Finally, Sect. 6 concludes this paper.
Notations: We denote vectors by boldface lower-case symbols and matrices by boldface upper-case symbols. The support set of a sparse signal is denoted as supp(·). The empirical mean of a signal is denoted as �·� and the standard deviation is represented as std(·) . || · || F is the Frobenius norm and vec(·) denotes the vectorization of a matrix by stacking the columns. The lower bound and upper bound on the complexity of an algorithm are represented by and O , respectively.

Automotive FMCW radar signal model
The chirp sequence modulation [13] which is a currently widespread variant of FMCW in automotive radar system, contains M consecutive chirps, and its transmit waveform can be represented as The individual transmit chirp signal is given by where f c is the carrier frequency, α = B/T denotes the slope of the transmit signal with B and T denoting the sweep bandwidth and the chirp duration, respectively, and rect T (t) is the square pulse of duration T. The mth chirp of the target's echo signal is delayed by τ relative to the transmit signal with a normalized relative Doppler shift d [2]: where A m is the received amplitude, τ = 2R/c , d = 2v/c , and v m denotes complex Gaussian noise. Here, R and v denote, respectively, the distance and relative radial velocity between the radar and the target, and c represents the speed of light. The beat signal at the intermediate frequency can be obtained after stretch processing, namely mixing r m with the complex conjugate of the transmitted signal: ŷ m (t) = r m (t)x * (t) . After filtering and sampling with a period of T s and collecting N samples per chirp, the discrete beat signal can be approximated as [2] With the assumption that y m,n represents the additive in-band interference, which can be induced by various sources of interference [14,15], the discrete beat signal y = {y (m·N +n) } for 0 ≤ m ≤ M − 1 and 0 ≤ n ≤ N − 1 can be therefore summarized as where H is the set of sample indices without interference and m · N + n / ∈ H denotes the indices of the samples containing interference. Figure 1 shows examples of interferencefree and interference-contaminated discrete beat signals. We now explain the two different signal models which will play an essential role in this paper.

Summary of signal model in matrix form
To obtain target range and velocity, the 2D discrete Fourier transform (DFT) can be applied on the discrete beat signal matrix Y = {y n,m } ∈ C N×M , where W N ∈ C N×N , W M ∈ C M×M are DFT matrices (as defined in Appendix A), and X ∈ C N×M denotes the 2D RD spectrum. Figure 4a shows an example RD spectrum without interference and Fig. 4b shows an example RD spectrum with interference.

Summary of signal model in vector form
Based on the 2D signal model in matrix form, the one-dimensional (1D) vector form of the signal model can be formulated. x ∈ C MN×1 denotes the vector form of the RD spectrum and can be obtained by where ⊗ represents the Kronecker product and y = (y 0 , ..., y MN −1 ) T ∈ C MN×1 .

Compressive sensing
Compressive sensing has shown its strength in reconstructing sparse signals using far fewer samples than required by the Nyquist sampling theorem [16]. It requires a transform domain that provides a sparse representation of the observed signal. Its sensing structure, i.e. the measurement transform matrix, has to satisfy the restricted isometry property [17]. The transform domain with the low-dimensional representation is called sparse domain. A signal having k nonzero coefficients in the sparse domain is called k-sparse. Generally, the sparsity of the signal y is measured by the The representation of the beat signal in most automotive radar application scenarios is sparse in the RD spectrum. Hence, the interference-pruned discrete beat signal can be seen as the beat signal with a reduced sampling rate and its sparse representation can be restored by the CS algorithm. The 1D signal model in vector form can be easily connected to the CS framework. The measurement transform matrix can be written as The dimensions of the inverse DFT (IDFT) matrices W N and W M are N × N and M × M (see Appendix A), respectively. Assume that the number of all undisturbed samples across all chirps in (5) is q. Thus, the resulting beat signal vector is given by ỹ = (y i 0 , ..., y i q−1 ) T with 0 < q < MN where ψ i q−1 denotes the i q -th row vector in matrix W . The radar interference mitigation problem can be rephrased as the reconstruction of the sparse vector x from the noisy measurement ỹ =˜ x . This problem is equivalent to solving an underdetermined set of linear equa- tions. An illustration of how to utilize the 1D CS framework for the interference mitigation of the automotive radar signal is presented in Fig. 2. Under the condition that x is sparse, the problem can be reduced to a minimization problem: with the Lagrange multiplier ν [18]. However, due to the discrete and discontinuous nature of the l 0 pseudo norm, the l 0 -minimization is NP-hard in general [18]. Thus, P 0 is computationally intractable. The l 1 -minimization or basis pursuit [19] can be interpreted as the convex relaxation of the l 0 -minimization, and the l 1 -minimization is as follows Since P 1 is convex, efficient solvers can be used, such as iterative shrinkage-thresholding pursuit [20]. Alternative reconstruction algorithms include greedy-type methods such as OMP [21,22], as well as thresholding-based methods [23][24][25] and the AMP algorithm [26]. These algorithms can be easily integrated in the framework shown in Fig. 2. However, the efficiency of this framework is limited as it vectorizes the 2D signal measurement in automotive radar system to a 1D vector of dimension MN.

Restricted isometry property
In order to successfully recover a good estimate of signal x , the measurement transform matrix ˜ in (10) should satisfy the restricted isometry property [17]. The selection of the measurement transform matrix has been analyzed in [18,27]. It is shown that a random partial Fourier matrix satisfies a near-optimal RIP with high probability [27,28]. In this work, the theoretical analysis on RIP of the measurement transform matrix ˜ is further conducted in Lemma 1 and Theorem 1 given in Appendix B. Theorem 1 shows that ˜ satisfies the RIP with high probability.

Structure of 1D algorithms for RD spectrum recovery
In this subsection we introduce the novel prior-model-based iterative sparsity-promoting algorithms employed in conjunction with the 1D formulation of our interference mitigation approach.

Prior-model-based iterative thresholding
When it comes to solving large-scale systems of linear equations often iterative gradient descent methods are employed rather than the Gaussian elimination [29]. The basic form of the update step of such an iterative algorithm for solving underdetermined systems like (9) and (10) is given by Here ˜ T r t =˜ T (ỹ −˜ x t ) represents the gradient of the approximation error (residual r t ) in the t-th iteration of the algorithm. T (·) denotes a nonlinear function that promotes the sparsity of the solution. The parameter 0 < ϑ ≤ 1 influences the convergence speed. Hard thresholding and soft thresholding have been shown to be two possibilities for the nonlinear function T (·) , where denotes the threshold value.
For the algorithms to run as efficiently as possible, the choice of the threshold value is crucial. It can remain constant for all iterations, decrease by a fixed multiplicative factor in each iteration, or be adaptively adjusted to the signal properties in each iteration. An adaptive adjustment can be derived by considering the gradient term ˜ T r t . Assuming that the values in this term correspond to a Gaussian distribution with zero mean and standard deviation std (˜ T r t ) , the threshold is given by where β is the threshold control parameter, typically in the range 2 < β < 4 [30]. This effectively reduces the noise in the representation by assuming that small signal values are part of the noise.
Since the typical measurement cycle of an automotive radar sensor is around 50 milliseconds [12], the observed movements of the targets in the RD spectra of successive measurement cycles are rather small in most application scenarios. Therefore, the prior information about the positions of the targets can be utilized to expedite the update steps of the iterative thresholding algorithms described in (11).
Instead of using the standard IST [16,31] and IHT [32], we incorporate prior information into the thresholding process for solving (9) and (10). The prior-model-based soft thresholding function is defined as: where i denotes the index of the elements of vector θ . The prior probability p i ensures that if a sparse maximum is likely at the i-th position, the threshold is scaled down accordingly by (1 − ζ(p i )) . This helps to reduce the number of iterations for searching an optimal estimate of x and facilitates the detection of local maxima, thereby reducing the reconstruction error. The mapping function is introduced for regulating the prior model, where a, b, e ∈ R are the control parameters.
Similarly, the prior-model-based hard thresholding function is defined as: where i denotes the index of the elements of vector θ.

Determination of the prior probability
The prior of the presence of the target at the i-th position of x is assumed to follow a normal distribution whose expected value µ i and variance σ 2 i are equal to the empirical mean and variance of the peak values at this position in the latest Q-measurements: where ξ = ξ 0 ∪ ξ 1 ... ∪ ξ Q−1 represents the set of positions of target peaks detected by a cell averaging constant false alarm rate (CA-CFAR) algorithm in the (original or recovered) RD spectra 1 of the latest Q measurements. ξ 0 and ξ Q−1 denote the sets of detected positions of target peaks in the RD spectrum of the current measurement x η and the measurement at time η − Q + 1 , respectively. The prior probability for the presence of the target at the i-th position ( i ∈ ξ ) of the next measurement x η+1 is determined by x i in x η . Because the presence of target peaks at other positions ( i / ∈ ξ ) was not observed in the latest Q-measurements, the prior probability of the these positions is initially set to zero. Then, a prior probability matrix P ∈ R N×M can be constructed. However, since the targets may move slightly in the next measurement cycle, the prior probability should optimally be propagated from the historical target positions to the neighboring positions surrounding them. The new prior probability matrix is then recalculated as P = G ⊛P , where G represents a 2D window function and ⊛ denotes the convolution operator. p i is then the i-th element of vec(P).

Integration of 2D masked residual updates
Since the multiplication of a vector with the Fourier matrix can utilize the fast Fourier transform (FFT), it can significantly improve the efficiency of recovery algorithms. The microcontrollers of most automotive radar sensors have an accelerator for FFT processing with reduced computational latency. However, the previously discussed radar interference mitigation framework that vectorizes the radar measurement to match the general CS framework cannot take advantage of this benefit. More precisely, recalling the framework illustrated in Fig. 2, by removing the interference-contaminated measurement samples in y , the corresponding rows of the measurement transform matrix W are pruned. Then, the remaining measurement signal ỹ and the pruned measurement transform matrix ˜ are used with different CS solvers to compute a sparse solution of x . Therefore, the FFT operator cannot be directly incorporated. In order to utilize the computational advantage of the FFT, a 2D masked residual updates framework is proposed that can be easily incorporated into existing CS solvers.
Recalling the residual updates in (11), the choice of ˜ ∈ C q×MN and r t ∈ C q have a clear dependency on the q-rows of interference-free samples. As the positions of interference-contaminated samples can be random [33], the size of ˜ can vary with interference scenarios. However, the residual updates always correspond to the remaining interference-free samples. Therefore, it is possible to control the residual updates using (16)  p·M ← − − R t . We refer to this residual updates process as 2D masked residual updates. From Theorem 1 (see Appendix B), it is known that a random row sub-matrix of W satisfies a near-optimal RIP with high probability. For the analysis of the RIP condition of the 2D MRU framework, we consider the equation Y rec =W N · X t ·W T M instead of the 1D formulation y rec =W · x t . The theoretical guarantee of successful RD spectrum recovery discussed in Theorem 1 also applies to 2D algorithms, where q is now substituted by p · M . More concretely, with the help of the mask matrix B , Y t obtains values from Y rec only at the selected p positions in each chirp, meaning that the valid updates of Y rec are preserved in these p · M entries (the same for y rec ). In this way, the rows with the same indices as the indices of these p · M entries are "subsampled" from W in this particular manner for the 2D case. The advantage of incorporating 2D MRU is that the matrix-vector multiplications for Fourier transforms can be solved quickly with hardware acceleration in the automotive radar system. Algorithm 1 describes the incorporation of 2D MRU with prior-model-based IST/IHT (PM-IHT/PM-IST) for RD spectrum recovery, where ǫ is the threshold parameter of the relative residual update. The proposed 2D MRU framework can be easily integrated with other well-known CS solvers, e.g., fast iterative shrinkage-thresholding algorithm (FISTA) [20], OMP [21], compressive sampling matching pursuit (CoSaMP) [22], YALL1 [34], AMP [26], and generalized AMP (GAMP) [35], since these solvers make use of residual updates as well.

Computational complexity
The largest part of the computational effort of the 2D algorithms is taken by the nested FFT in each case. Since FFT in column direction is performed for all rows of the matrix and FFT in row direction is performed for all columns of the matrix, the total computational effort for 2D MRU in each iteration loop can be described by The computational effort of 1D solvers in each iteration loop can be described by  O(qMN ) , where q denotes the number of all undisturbed samples across all chirps. With the proposed approach the complexity is reduced considerably by a factor of O log(NM)/q . Considering a discrete beat signal Y ∈ C 128×256 where 10% samples are disturbed by interference (i.e., q equals 29491), log(NM)/q ≈ 1/6531.

Detection of disturbed samples
In [36], a method was presented in which the interference is detected by iterative adaptive thresholding. Samples from a signal vector y of length L are considered to be interference-contaminated samples when their absolute values exceed a threshold and the detected   where γ denotes the control parameter. If the threshold T D changes more than a predefined value T D , the newly calculated threshold is used to detect further interferencecontaminated samples. As soon as T D changes by less than T D between two iterations, the algorithm terminates. In [10], a classical edge detector, namely, a Laplacian filter [37] is used to identify the anomalies and the interference-contaminated samples. The Laplacian filter based method is fast and hardly deletes interference-free samples. Still, it misses some interference-contaminated samples, since the Laplacian filter assigns a larger weight to edges. Thus, if the interference-contaminated samples are grouped, the distorted samples in the middle often cannot be detected completely. However, this flaw can be compensated by a combination with the iterative adaptive thresholding detection method. To quantitatively analyze the detection performance of different algorithms, the recall and F-measure are introduced, which are calculated as: Recall = TP TP+FN and F = 2 TP 2 TP + FP + FN where TP, FP, and FN denote the number of true-positive, false-positive, and false-negative estimates.
Since the interference-contaminated samples are the minority in the discrete beat signal, the F-measure is theoretically more important for evaluating the detection performance. However, in radar interference mitigation with CS, the amount of correctly detected interference positions has a greater impact on the recovery results. If a small amount of interference-free samples is accidentally discarded, this results in a small change in the compression ratio as defined in Sect. 5 and does not have a large impact on the recovery results. The evaluation results in Table 1 show that the combination of the two methods can correctly detect about 96% of the interference-contaminated samples. The combined method is therefore more suitable for interference mitigation with CS approaches. The undetected 4% of the interference-contaminated samples should have low amplitudes, as they would otherwise have been detected by the iterative adaptive thresholding method. Thus, these distorted samples may also generate little additional noise in the frequency domain and therefore do not affect the weak targets. Figure 3 shows an example of using this combined detection method to determine the position of the anomalies. The combined method due to its superior performance in terms of the true-positive rate (recall) is used in this work for the detection of interference-contaminated samples.

Evaluation methods and results
In this section, the performance of algorithms incorporating 2D MRU is analyzed and the evaluation results are further discussed. Firstly, the reconstruction performance of RD spectrum of real radar measurements is evaluated in terms of the run time and signal-to-interference-plus-noise ratio (SINR) in comparison with the state-of-theart algorithms. Secondly, the time consumption and the relative reconstruction error are evaluated under different metrics, namely, the compression ratio, the amount of (19)  remaining interference-contaminated samples used in reconstruction, and the sparsity level ratio. The compression ratio δ and the sparsity level ratio ρ are defined as where p is the number of samples selected for reconstruction in each chirp and k denotes the number of target peaks. For the prior-model-based iterative thresholding algorithms, the prior probability is determined based on the latest five measurement cycles (i.e., Q = 5 in (16) for evaluations in this work).

Performance evaluation on real radar measurements
Recent research results have shown that the AR model [6] and CS algorithms like the iterative method with an adaptive thresholding for compressed sensing (IMATCS) [7] or the YALL1 algorithm [38] provide satisfactory recovery results of the discrete beat signal for interference mitigation. The performance of these methods along with the simple zeroing method [1], and the prior-model-based iterative thresholding algorithms (PM-IHT and PM-IST) are evaluated in terms of the run time and SINR on a real radar measurement. The real radar measurement is recorded in a test chamber with radiationabsorbent materials and the targets (at distances of 27m, 50m, 73m, 95m and 120m) are created by a target generator. The interference is produced by a radar that emits signals in the same frequency band with a larger slope than the victim radar. Figure 4 shows the recovered RD spectra after the elimination of the severely distorted samples (ca. 10% of the total measurement) for four reference algorithms, AR, IMATCS, YALL1 2 and the zeroing method, as well as for the proposed prior-model-based iterative thresholding algorithms. The recovery performance of IMATCS may vary slightly when different thresholds are selected. It should be clarified that the parameters of IMATCS for implementation are chosen according to [7], where the value of the highest peak is used to initialize the linear threshold. The targets in the distorted RD spectrum are hard to detect, however, the targets in the recovered RD spectra can be found easily. Here, the threshold parameter of relative residual update ǫ of the CS algorithms is set to 10 −6 . The AR, IMATCS, and YALL1 algorithms can properly restore the RD spectra, while the zeroing method produces many artifact peaks. The prior-model-based iterative thresholding algorithms, however, produce superior recoveries. Since the proposed 1D and 2D prior-model-based iterative thresholding algorithms differ only in their computational load, the recovered RD spectra of the prior-model-based iterative thresholding algorithms in Fig. 4 are solely presented for the 2D case.
Since the amount of disturbed samples in the discrete beat signal may vary under different types of interference, the robustness of the algorithms will be further evaluated by eliminating more samples. Table 2 shows the performance comparison of different algorithms in terms of the run time and SINR for different amounts of interference and correspondingly discarded interference-contaminated samples. For interference detection, the combined method discussed in Sect. 4 is used in this evaluation.  CPU@1.70GHz. It is found that the run time of the algorithms PM-IST and PM-IHT incorporating 2D MRU is less than 50 milliseconds for a dimension of Ỹ up to 128 × 256 , while the run time of the 1D implementation requires seconds to complete the recovery. The code implementation of IMATCS 3 is modified based on [39] in which every iteration step needs the eigenvalue decomposition of the measurement transform matrix. Thus, the IMATCS algorithm takes more time for recovery as the size of the measurement transform matrix is large. The size of the measurement transform matrix decreases accordingly when the size of the remaining discrete beat signal decreases. This explains why the run time of IMATCS decreases as more samples are discarded. The run time of the zeroing method is instantaneously fast and therefore difficult to measure accurately. The other advantage of the prior-model-based iterative thresholding algorithms is that they provide superior SINR improvement, since the proposed PM-IST and PM-IHT can strongly suppress the measurement noise. This can be verified in Fig. 5, where the velocity cut of the recovered RD spectrum of the target at 73 meter is shown. The signal peak of the target at 73 m is much easier to detect in the recovered RD spectra provided by the PM-IHT and PM-IST methods.

Performance evaluation of 2D algorithms on different metrics
Since different metrics such as the amount of interference-free samples used for reconstruction, the sparsity level of RD spectrum and the amount of remaining interferencecontaminated samples all have an impact on the performance of different CS solvers [40], the performance of 2D algorithms are therefore further evaluated under these metrics. In this evaluation part, the representative greedy algorithm CoSaMP [22] and AMP with the proposed 2D extension are also included. Firstly, the computation time and the mean relative absolute error (MRAE) of 2D algorithms for different compression ratio (20) is examined. The different amount of interference-free samples is used to recover the RD spectrum. The value of δ is set between 0.1 and 1. For each δ , the simulation is repeated fifty times and the threshold parameter ǫ is set to 10 −3 . The MRAE is defined as the mean of the relative absolute error between the target peaks in the interference-free RD spectrum X clean and in the reconstructed RD spectrum X rec : where ñ and m are row and column indices of the RD matrix, Z is the set of target peaks. # Z denotes the cardinality of set Z. Figure 6 shows that as more interference-free samples p of the signal are used for reconstruction, δ becomes larger and the time needed for the execution of the algorithms decreases. This is a logical consequence of the fact that the number of candidate solutions of an underdetermined system increases and thus more iterations are required in order to find an optimal solution when fewer interference-free samples remain. Figure 7 shows that MRAE of all target peaks (except for the target peaks recovered by YALL1, which aims to achieve a high accuracy and therefore also tries to reconstruct some noise [34]) in the restored RD spectra becomes smaller as the value of δ increases. This shows that the more interference-free samples are left, the more accurately the target peaks will be restored. The proposed PM-IHT and PM-IST further improve the MRAE of the standard IHT and IST, although the improvements in run time are not significant. It should be noted that MRAE provides a measure of the restoration quality of target peaks, namely the sparse representations. There is therefore no direct connection between MRAE and SINR, since it is possible that an algorithm recovers the target peak with a small MRAE but a high background noise level, so that the SINR is small. Since it is difficult to accurately identify all interference-contaminated samples, in some scenarios few interference-contaminated samples may remain after the interference detection. Therefore, we secondly evaluate the robustness of the algorithms in the case that some of the remaining samples contain interference. To simulate the reconstruction on these scenarios, different interference-free beat signals with an average sparsity level ratio ( ρ = 0.005 ) are generated and superimposed with interference. The simulations are performed with a δ that is always adjusted to use the maximum number of interference-free samples for reconstruction. Amounts of interference-contaminated samples between 0.1% to 15% of total remaining samples are tested. Since the distorted sample usually contains larger amplitudes, it changes the SINR of the RD spectrum of the input data. Empirically, it is shown that when the percentage of distorted samples in Ỹ ranges from 0.1 to 15%, the SINR (originally ca. 32 dB) decreases accordingly by a value ranging from 7 dB to 25 dB. Furthermore, as shown in Fig. 8, the number of iterations increases slightly for all algorithms, as the increasing interference quantity can influence the search direction of the sparse solution of the algorithm. Although CoSaMP requires the smallest number of iterations, its run time is still longer than that of iterative thresholding algorithms due to the relatively long computation time in each iteration. In general, the run time of the iterative thresholding algorithms is less than 20 milliseconds. Figure 9 shows that the MRAE of all target peaks also increases slightly with the amount of remaining interference-contaminated samples. This is due to the fact that the remaining interference is also used for reconstruction, thus the recovery contains more errors. Overall, the relative errors are nevertheless low with values below 6% and CoSaMP has the smallest MRAE. Thirdly, the concordance of the local maxima and the MRAE of 2D algorithms for different sparsity level (20) is examined. The value of δ is set to 0.5. The sparsity level ρ is set between 0.001 to 0.1. For each ρ the simulation is repeated fifty times and the threshold parameter ǫ is set to 10 −3 . Figure 10 shows that the MRAE of PM-IHT, PM-IST and AMP increases slightly as ρ increases. The value of ρ has more effect on CoSaMP, and interestingly, the MRAE of YALL1 decreases while ρ increases. Another metric worth evaluating is the concordance of the local maxima. It is important to check whether all local maxima are reconstructed, since the local maxima represent the individual targets. Moreover, for radar target tracking, the position of targets in the RD spectrum is decisive. To generate an evaluation criterion from this context, all local maxima are determined in the RD spectrum of the interference-free signal as well as in the RD spectrum of the reconstructed signal. Then, the percentage of the maxima from the interferencefree signal that are present in the reconstruction is calculated. As shown in Fig. 11, the concordance of local maxima decreases as the sparsity level ratio ρ increases, i.e., more valid targets are present. The concordance of the local maxima of PM-IHT can achieve around 86% when ρ equals 0.1 while concordance of the local maxima of the other algorithms are still better than 80%. This means that more than 80% of the targets are still recognizable after recovery. This shows the robustness of these 2D CS solvers. Note that ρ is determined by the reference signal without interference, and ρ = 0.1 actually corresponds to a large number of target peaks, since the denominator in ρ (see (20)) is large.

Discussion
In comparison with the 2D framework, the 1D framework is limited not only in computation time but also in hardware resources in terms of memory, since the measurement transform matrix ˜ typically has a large size. For example, for a radar measurement of dimension Y ∈ C 512×128 [4], where 10% samples of the total measurement are distorted, the corresponding matrix ˜ ∈ C 58982×65536 becomes excessively large. In contrast to the 1D case, the 2D framework avoids additional memory for the measurement transform matrix by utilizing the built-in FFT operator. Although 2D MRU can also be integrated into the baseline method YALL1, its performance in terms of SINR and the run time is still worse than the greedy algorithms (CoSaMP) and the proposed prior-model-based iterative thresholding algorithms in general. The prior-model-based IHT shows its robustness both in the evaluation of real radar measurements and in the evaluation of CS metrics. Its run time can be even faster when the threshold parameter ǫ is relaxed to 10 −2 , for example. Besides the benefits in terms of run time and reconstruction error, the prior-modelbased iterative thresholding algorithm also has the potential to help in the reconstruction of weak targets. As shown in Fig. 12, a weak target (within the red dashed square with a power of ca. 5 dB) is slightly shifted from its position in Fig. 12c to the new position in Fig. 12a. Since its power is relatively close to the background noise compared to the other target peaks, which range from about 7 dB (bottom left) to 37 dB (top left), it cannot be properly reconstructed with the standard IHT algorithm. However, in the prior-model-based IHT, due to the presence of the target peak (ca. 9 dB) in the previous measurement, the probability of the presence of the target in the red dotted square area is increased for the next measurement cycle (namely, the probability of the presence of the target is propagated to the whole red dotted square area by convolution with a window function) and the thresholds for this area are decreased accordingly. This weak target is therefore properly reconstructed, and in the meantime, the false-positive prior information does not affect the reconstruction result, as shown in Fig. 12e.
Thus, the proposed 2D framework in conjunction with prior-model-based IHT can be a good candidate for automotive radar interference mitigation.

Conclusion
In this paper, an effective 2D masked residual updates compressive sensing framework as well as the prior-model-based iterative thresholding algorithms are proposed for automotive radar interference mitigation. 2D masked residual updates can be easily incorporated into the state-of-the-art compressive sensing algorithms and reduces the computational complexity of these algorithms. In particular, the proposed 2D masked where ω N := exp − 2π j N is the N th root of unity. The IDFT matrix W N ∈ C N ×N can be obtained by replacing the ω N with ω * N := exp 2π j N .

Lemma 1
Given two unitary IDFT (or DFT) matrices W N ∈ C N ×N and W M ∈ C M×M , the Kronecker product of these two IDFT (or DFT) matrices W =W M ⊗W N is also a unitary matrix.
Proof Since the IDFT (or DFT) matrices W M ∈ C M×M and W N ∈ C N ×N are both invertible, also their Kronecker product W =W M ⊗W N is invertible [41] and W −1 = (W M ⊗W N ) −1 =W −1 M ⊗W −1 N . As shown in Definition 1, W M and W N are both unitary matrices, thus W * M =W −1 M and W * N =W −1 N . As conjugate transposition is distributive over the Kronecker product, it follows that W −1 =W * M ⊗W * N = (W M ⊗W N ) * =W * . This means that W is a unitary matrix.

Proof of Theorem 1
Definition 2 A complex matrix A ∈ C M×N satisfies the restricted isometry property of order k with restricted isometry constant δ k > 0 , if For some q = O(log 2 (1/δ k )δ −2 k · k · log 2 (k/δ k ) · log(MN )) , let ˜ ∈ C q×MN be a matrix whose q rows are chosen uniformly at random, and independently from the rows of √ MN /q ·W . Then, matrix ˜ satisfies the RIP of order k with probability 1 − 2 −�(log(MN )·log(k/δ k )) . Proof From Lemma 1, it is known that the Kronecker product of two IDFT or DFT matrices results in a unitary matrix. [28] gives the proof that random sub-matrices of unitary matrices satisfy a near-optimal RIP with high probability. Combining Lemma 1 and the proof in [28], we can draw the conclusion of Theorem 1.