Variational Signal Separation for Automotive Radar Interference Mitigation

Scenarios for Interference with FMCW Automotive Radars

 Variational Signal Separation for Automotive Radar Interference Mitigation

Electrical Engineering and Systems Science > Signal Processing

[Submitted on 23 May 2024]

Mate Toth, Erik Leitinger, Klaus Witrisal

Algorithms for joint mutual interference mitigation and object parameter estimation are a key enabler for automotive applications of frequency-modulated continuous wave (FMCW) radar. The underlying signal model poses a challenge for signal separation, since both the coherent radar echo and the non-coherent interference influenced by individual multipath propagation channels must be considered.

In particular, under certain assumptions, the model is described as a superposition of multipath channels weighted by parametric chirp envelopes in the case of interference. In this paper, we introduce a method inspired by sparse Bayesian learning (SBL) to detect and estimate radar object parameters while also estimating and successively canceling the interference signal. An augmented probabilistic model is employed that uses hierarchical Gamma-Gaussian prior model for each multipath channel separately.

Based on this model an iterative inference algorithm is derived using the variational expectation-maximization (EM) methodology. The algorithm is statistically evaluated in terms of object parameter estimation accuracy and robustness, indicating that it is fundamentally capable of achieving the Cramer-Rao lower bound (CRLB) with respect to the accuracy of object estimates and it closely follows the radar performance achieved when no interference is present.

Comments:	18 pages, 8 figures, submitted to IEEE Transactions on Radar Systems on 23rd of May, 2024
Subjects:	Signal Processing (eess.SP)
Cite as:	arXiv:2405.14319 [eess.SP]
	(or arXiv:2405.14319v1 [eess.SP] for this version)
	https://doi.org/10.48550/arXiv.2405.14319

Submission history

From: Mate Toth [view email]
[v1] Thu, 23 May 2024 08:49:55 UTC (897 KB)

Summary

The paper introduces a new method for mitigating automotive radar interference while jointly estimating object parameters, inspired by sparse Bayesian learning (SBL). Key points:

• Develops a probabilistic signal model representing both the coherent radar echo and non-coherent interference influenced by multipath propagation channels. The model is a superposition of multipath channels weighted by parametric chirp envelopes for interference.
• Proposes an iterative variational expectation-maximization (EM) inference algorithm that detects and estimates radar object parameters while estimating and canceling the interference signal.
• The algorithm employs an augmented probabilistic model using a hierarchical Gamma-Gaussian prior for each multipath channel separately.
• Evaluation shows the algorithm can achieve the Cramer-Rao lower bound on object estimation accuracy and closely follows radar performance without interference present. It outperforms established interference mitigation preprocessing methods.
• The method is relatively resilient to model mismatch and poor signal separability scenarios. Preliminary proof-of-concept is provided on measurement data.
• Though presented for automotive radar interference, the proposed algorithm may apply to other similar signal separation problems.

In summary, the paper presents a novel, robust method combining SBL-inspired probabilistic inference with signal modeling to jointly mitigate radar interference and estimate object parameters. 

Algorithm

The proposed algorithm is based on a variational expectation-maximization (EM) approach, inspired by sparse Bayesian learning (SBL), to infer the superposition of sparse line spectra representing the radar object signal and interference.

The main steps of the algorithm are:

• 1. Initialize estimates for the object dictionary, object weights, interference dictionary, interference weights, and noise precision.
• 2. Update the interference estimate:
•    - For each radar ramp, update the chirp parameters and run the interference signal estimation subroutine.
•    - Update the noise precision estimate and object signal residual.
• 3. Update the object signal estimate:
•    - Run the object signal estimation subroutine.
•    - Update the noise precision, interference signal residual, and interference weight estimates.
• 4. Repeat steps 2-3 until convergence or a fixed number of iterations is reached.

The interference and object signal estimation subroutines are similar in structure. They involve:

1. Finding a possible new component to add to the dictionary.
2. Updating all existing components in the dictionary.
3. Updating the respective weight and parameter estimates.

A key aspect is the fast component precision update, which allows for adaptive pruning of dictionary components or addition of new components based on a signal-to-noise ratio threshold.

The algorithm jointly estimates the object range-Doppler parameters, interference multipath channel parameters, and signal parameters. It employs a bottom-up initialization scheme, starting with empty dictionaries for computational efficiency.

Specific design choices, such as grid-less estimation for the object channel and on-grid estimation for the interference channel, as well as the alternating update of object and interference estimates, contribute to the algorithm's robustness and ability to avoid local optima.

The algorithm's derivation involves a structured mean-field variational approximation and optimization of the evidence lower bound (ELBO). Update equations for the channel weights, noise precision, and unknown parameters are obtained by maximizing the ELBO.

Data, Tests/Sim and Artifacts

The paper presents both simulation studies and a real-world measurement example to evaluate the proposed algorithm.

Simulations:

• 1. Single-object estimation: A simplified scenario with a single radar transmit ramp (delay-only estimation) is used to compare the algorithm's performance to the Cramer-Rao lower bound (CRLB). The object channel consists of a single line spectral component, while the interference channel is constant (direct-path interference only). The root mean squared error (RMSE) of the normalized beat frequency estimate is analyzed over varying signal-to-noise ratio (SNR) and signal-to-interference ratio (SIR) levels. Three algorithm variants are considered: neglecting interference, the proposed factorized algorithm, and joint modeling of object and interference channels.
• 2. Multi-object estimation: A more complex scenario with multiple radar ramps, objects, and interference components is simulated. The generalized optimal subpattern assignment (GOSPA) metric is used to evaluate multi-target detection performance, considering localization errors, misdetections, and false alarms. The proposed algorithm is compared to two existing interference mitigation methods: zeroing and morphological component analysis (MCA).

Real-world measurements:

A qualitative example using automotive radar measurement data from a real traffic scenario is presented. As the measurements did not inherently contain interference, a separately measured interference signal was added to the data. The range-Doppler spectrum and detected objects are compared for the interference-free, interfered, and interference-mitigated cases.

Artifacts for independent validation:

The paper does not explicitly provide artifacts for independent validation. However, it includes detailed descriptions of the signal model, inference algorithm, and simulation parameters, which could potentially allow for replication of the results. The measurement example is qualitative and may be more challenging to reproduce without access to the original data.

In summary, the paper presents simulation studies demonstrating the algorithm's performance in terms of object parameter estimation accuracy and multi-target detection metrics, as well as a qualitative real-world measurement example. While detailed, the provided information may not be sufficient for complete independent validation without access to the original code and data.

Authors

The authors of the paper are Mate Toth, Erik Leitinger, and Klaus Witrisal. Their affiliations and previous related work are as follows:

Mate Toth:
- Affiliation: Institute of Communication Networks and Satellite Communications, Graz University of Technology, Graz, Austria.
- Previous related work: Toth has co-authored papers on automotive radar interference analysis and mitigation, including "Analytical Investigation of Non-Coherent Mutual FMCW Radar Interference" (2018) and "Slow-Time Mitigation of Mutual Interference in Chirp Sequence Radar" (2020).

Erik Leitinger:
- Affiliation: Laboratory of Signal Processing and Speech Communication, Graz University of Technology, Graz, Austria; associated with the Christian Doppler Laboratory for Location-aware Electronic Systems.
- Previous related work: Leitinger has co-authored papers on sparse signal processing and estimation, such as "Super-resolution estimation of UWB channels including the dense component – An SBL-inspired approach" (2024, preprint) and "Detection and estimation of a spectral line in MIMO systems" (2020).

Klaus Witrisal:
- Affiliation: Institute of Communication Networks and Satellite Communications, Graz University of Technology, Graz, Austria; associated with the Christian Doppler Laboratory for Location-aware Electronic Systems.
- Previous related work: Witrisal has co-authored papers on radar signal processing, localization, and channel modeling, including "Analysis of automotive radar interference mitigation for real-world environments" (2020) and "Deep interference mitigation and denoising of real-world fmcw radar signals" (2020).

The authors' affiliations suggest a collaboration between researchers from the Institute of Communication Networks and Satellite Communications and the Laboratory of Signal Processing and Speech Communication at Graz University of Technology. Their previous work demonstrates expertise in automotive radar interference, sparse signal processing, and channel modeling, which are relevant to the proposed algorithm in this paper.

The Christian Doppler Laboratory for Location-aware Electronic Systems, with which Leitinger and Witrisal are associated, is mentioned as a supporting institution for this work.

======================================

• Variational Signal Separation for Automotive Radar Interference Mitigation

  Mate Toth, Erik Leitinger, and Klaus Witrisal M. Toth and K. Witrisal are with the Institute of Communication Networks and Satellite Communications, Graz University of Technology, Graz, Austria (e-mail: (mate.a.toth,witrisal)@tugraz.at). E. Leitinger is with the Laboratory of Signal Processing and Speech Communication, Graz University of Technology, Graz, Austria (e-mail: erik.leitinger@tugraz.at). E. Leitinger and K. Witrisal are further associated with the Christian Doppler Laboratory for Location-aware Electronic Systems. This work was supported in part by the Christian Doppler Research Association.

  Abstract

  Algorithms for joint mutual interference mitigation and object parameter estimation are a key enabler for automotive applications of frequency-modulated continuous wave (FMCW) radar. The underlying signal model poses a challenge for signal separation, since both the coherent radar echo and the non-coherent interference influenced by individual multipath propagation channels must be considered. In particular, under certain assumptions, the model is described as a superposition of multipath channels weighted by parametric chirp envelopes in the case of interference.

  In this paper, we introduce a method inspired by sparse Bayesian learning (SBL) to detect and estimate radar object parameters while also estimating and successively canceling the interference signal. An augmented probabilistic model is employed that uses hierarchical Gamma-Gaussian prior model for each multipath channel separately. Based on this model an iterative inference algorithm is derived using the variational expectation-maximization (EM) methodology. The algorithm is statistically evaluated in terms of object parameter estimation accuracy and robustness, indicating that it is fundamentally capable of achieving the Cramer-Rao lower bound (CRLB) with respect to the accuracy of object estimates and it closely follows the radar performance achieved when no interference is present.

  Index Terms:

  line spectral estimation, sparse Bayesian learning, variational EM, signal separation, automotive radar, interference.

  I Introduction

  The growing number of radar-equipped vehicles is expected to cause mutual interference, potentially compromising safety-critical automotive radar applications [1, 2]. As a result, this issue has received significant research attention in the past decade. Large industry research projects have been conducted to analyze interference mitigation methods and strategies for using the automotive radar spectrum [3, 4]. From an academic perspective, the interference problem presents a challenge in signal modeling and processing, making it an interesting area for algorithm design. Early analyses were conducted in [5, 6, 7], and since then, a significant amount of research has been published on the subject.

  I-A State of the Art

  Interference modeling efforts can be categorized into scenario-level and signal-level models. Works in the former category characterize the probability of interference occurring and the received interference power, whereas the latter models the specific form of interference in the received radar signal. Scenario-level modeling approaches include analyses via stochastic geometry [8], graph theory [9] as well as other extensive numerical simulation frameworks [10]. Signal-level models have been developed for different modulation schemes, including the most commonly used FMCW [11, 12] as well as phase-modulated continuous wave (PMCW) [13] and orthogonal frequency division multiplexing (OFDM) [14] radars. Experimental evaluations can be found e.g. in [15, 16]. Recent modeling efforts also consider the interference propagation channel [17, 18]. This paper focuses on the mutual interference of FMCW radars, taking into account the interference channel.

  Approaches to mitigate interference include avoidance strategies using cognitive radar principles [19, 20] as well as exploiting vehicle-to-vehicle communication [21, 22]. In terms of signal processing algorithms to cancel the interference in the receiver, a vast number of algorithms have been proposed. As mutual FMCW interference is commonly of a short duration as compared to the total number of received samples, zeroing, optionally followed by interpolating the zeroed samples, has been discussed first in [6]. This method typically necessitates a separate prior interference detection step, e.g. [23, 24, 25, 26]. Interpolation approaches include autoregressive modeling [27], time-frequency interpolation [28] and sparsity-based algorithms [29, 30, 31]. Alternative approaches to interpolation include spatial beamforming [32, 33], adaptive noise cancellers [34, 35], time-frequency [36, 37] and nonlinear filtering [38, 39]. Denoising based on signal decomposition and thresholding has been proposed in e.g. [40]. Additionally, many works are making use of deep learning [41] to detect [42] and mitigate interference [43, 44, 45].

  If the interference signal is estimated, it may be canceled completely by subtraction. Proposed methods for this include [46, 47, 48]. However, if the useful radar return itself, here termed the object signal, is also explicitly considered in the estimation, the problem can be cast as signal separation. This problem is generally ill-posed, but can be solved by assuming that the object and interference signal components are represented by a sparse basis in different signal bases. In [49], the author applied the framework of morphological component analysis (MCA) [50] modeling the interference and object signals as being sparse in the discrete Fourier and discrete short-time Fourier domains, respectively. This assumption allows for the joint sparse recovery of the signals to be written as a modified dual basis pursuit optimization problem [51], for which an iterative algorithm was proposed. Similar algorithms can be found in [52, 53, 54, 55].

  I-B Contributions and Notations

  Although the existing literature is extensive, we identify the following open topics for which our work provides novel insights. Previously proposed algorithms are applied as preprocessing in a conventional range-Doppler processing chain. This approach is practical, but has inherent limitations. In general, the goal of radar sensing is to provide estimates of object parameters, which is equivalent to estimating the parameters of a multipath propagation channel containing static (such as buildings, traffics signs or trees) and moving objects (such as vehicles, bicycle or pedestrians). Note that parametric channel models typically represent multipath propagation as a linear superposition of weighted Dirac delta distributions - or spectral lines - with distinct supports in the underlying dispersion domain (range, angle of arrival, angle of departure, Doppler frequency, and combinations thereof). Therefore, estimating multipath parameters can be cast as a line spectral estimation (LSE) problem [56].

  To achieve robust inference performance, it is necessary to consider the entire problem of parameter estimation under the influence of interfering signals. Previous methods are also often based on heuristics or require the prior setting of certain algorithm parameters. This makes it necessary to manually fine-tune the range-Doppler processing chain and makes it difficult giving performance guarantees. Our proposed algorithm is derived within a rigorous theoretical framework. Specific contributions are the following:

  • •

    We develop a new probabilistic signal model for the problem of mutual FMCW radar interference. The model incorporates the multipath propagation as line spectra for both the coherently received radar echo and the interference.

  • •

    We propose a novel inference algorithm inspired by SBL [57, 58, 59] for LSE that is able to infer the superposition of sparse line spectra. The proposed algorithm, which is based on the variational EM approach, provides super-resolution by jointly estimating the objects’ range and Doppler parameters, as well as the parameters of the interference multipath channel and signal parameters. This makes it explicitly robust to mutual interference.

  • •

    We quantitatively show near-optimal multi-target detection and parameter estimation performance by statistically evaluating the proposed method in simulation, comparing it to the CRLB as well as to the interference-free case. We also showcase that the algorithm significantly outperforms a few established signal preprocessing methods for interference mitigation.

  • •

    We include investigations on difficult-to-handle scenarios of model mismatch and poor signal separability. We point out the relative resilience of the algorithm to these challenges, and provide an initial proof-of-concept based on measurement data.

  Notation: The complex conjugate of the variable $c$ is denoted by $c^{\ast }$. Bold lowercase letters denote column vectors and bold uppercase letters denote matrices, respectively. For the vector $𝒗$, $\Vert 𝒗\Vert$ is its Euclidean norm; $\mathrm{Diag}(𝒗)$ constructs a diagonal matrix from $𝒗$. The matrix transpose of $𝑴$ is ${𝑴}^{\mathrm{T}}$; ${𝑴}^{\mathrm{H}}$ is the Hermitian; $\mathrm{tr}(𝑴)$ is the trace of the matrix, and $|𝑴|$ denotes the determinant. The notation $𝑴[m,n]$ refers to the element $(m,n)$ of the matrix; $𝑰$ denotes the identity matrix. For the random variable $𝒙$, $p(𝒙)=\text{CN}(𝒙|𝝁,𝑪)$ denotes a complex Gaussian probability density function (PDF) with mean $𝝁$ and covariance matrix $𝑪$; $p(x)=\text{Ga}(x|a,b)$ is a Gamma PDF with shape $a$ and rate parameters $b$. The expectation of the probability distribution $p(x)$ with respect to $q(x)$ (denoted in shorthand as $p$ and $q$) is written as ${⟨p⟩}_q=\int p(x)q(x)dx$. The vertical bar denotes a conditional PDF $p(𝒙|𝒛)$; the semicolon a parametrized PDF $p(𝒙;𝜽)$. $\widetilde{𝜽}$ denotes the vector containing the “true” parameters of the generative model and $\widehat{𝜽}$ denotes an estimate of the parameter vector $𝜽$.

  II Signal Model

  II-A FMCW Radar Under Interference

  The victim radar transmits a transmit signal $x(t)$ in an environment containing $M_{\mathrm{I}}$ interferers with respective transmit signals $x_{\mathrm{I},i}(t),i\in [0,M_{\mathrm{I}}-1]$. The corresponding received intermediate frequency (IF) signal $r(n,p)$, where $n\in [0,N-1]$ denotes the fast-time index with sampling time ${\mathrm{T}}_{\mathrm{s}}$ and $p\in [0,P-1]$ is the slow-time ramp index¹¹We consider a single receive antenna of the radar system. The proposed model and algorithm are extendable to coherent multi-antenna processing, but further investigation is out of scope of the current work., is given by

  $r(t^{\prime }=n{\mathrm{T}}_{\mathrm{s}},p)=r_{\mathrm{O}}(n,p)+{\displaystyle \sum _{i=0}^{M_{\mathrm{I}}-1}}{r}_{\mathrm{I},i}(n,p)+\eta (n,p).$

  (1)

  The first term $r_{\mathrm{O}}(n,p)$ refers to the coherently received radar echo, termed the object signal. The non-coherent received interference signal from the $i$-th interferer is denoted by $r_{\mathrm{I},i}(n,p)$. Finally, $\eta (n,p)$ represents additive white Gaussian noise (AWGN) with precision (inverse variance) $\tilde{\lambda }$. Fig. 1a illustrates an example interference scenario where $M_{\mathrm{I}}=1$, and Fig. 1e shows example signals for a single ramp.

  

  Figure 1: Illustrations of the principles and signal modeling for a mutual automotive FMCW radar interference scenario. (a) A critical automotive radar scenario with a pedestrian to be detected and an interfering vehicle on the opposite lane. (b) Abstracted block diagram representing the applied signal model. (c)-(d) Time-frequency plots of the involved signals. (e) Plot of the received signal components from an interfered ramp.

  The victim radar, a monostatic FMCW radar, transmits a sequence of linearly modulated chirp signals (ramps)

  $x(t)=A{\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})u(t-T_{\mathrm{p}}p;f_0,k)$

  (2)

  with

  	$\overline{x}(t;T_{\mathrm{sw}})$	$=\{\begin{array}{cc}1\hfill & \text{for }0\le t\le T_{\mathrm{sw}}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$		(3)
  	$u(t;f_0,k)$	$=\exp \left(\mathrm{j}(2\pi f_{0}t+\pi k{t}^2)\right)$		(4)

  where we defined the parametrized rectangular function template $\overline{x}(t;T)$ and the chirp $u(t;f_0,k)$. $A$ is some complex transmit amplitude, $P$ is the number of chirps, $f_0$ is the chirp start frequency, $k$ is the chirp slope and $T_{\mathrm{p}}$ is the pulse duration which consists of the actual transmit chirp duration $T_{\mathrm{sw}}$ plus an additional idle time.²²Note that actual transmit signals are real-valued but a complex receiver architecture is assumed here. Thus, we immediately use a complex representation for mathematical simplicity. However, it is possible to model a real receiver architecture as well.

  The transmit signal $x(t)$ in (2) propagates through a multipath channel, yielding a superposition of reflected signals that are coherently demodulated by multiplication with $x^{\ast }(t)$, and filtered by the receiver anti-aliasing filter (AAF) $G(f)$. This is indicated in Fig. 1b starting at the rectangle symbol. The according received object signal $r_{\mathrm{O}}(t^{\prime },p)$ is given by

  	$r_{\mathrm{O}}(t^{\prime },p)$	$\approx H_{\mathrm{O}}(f_0+k{t}^{\prime },T_{\mathrm{p}}p)$
  		$={\displaystyle \sum _{l=0}^{\tilde{L}-1}}{\tilde{\alpha }}_l\exp (-\mathrm{j2}\pi (f_0+k{t}^{\prime }){\tilde{\tau }}_l)\exp (-\mathrm{j2}\pi {\tilde{\nu }}_l{T}_{\mathrm{p}}p)$		(5)

  with $H_{\mathrm{O}}(f,t)$ being the frequency-selective time-variant channel transfer function of the object channel made up of $\tilde{L}$ components³³The object channel describes the full radar measurement scene, where components originate from objects in the scene. Note that a single physical object can potentially generate multiple components. with respective amplitudes ${\tilde{\alpha }}_l$. A detailed derivation of (5) can be found in Appendix A. An illustration of a single object component in the time-frequency plane can also be seen in Fig. 1c. The fast-time dimension with $t^{\prime }=t-T_{\mathrm{p}}p\in (0,T_{\mathrm{sw}})\forall p$ corresponds to the beat frequencies or analogously the delays ${\tilde{\tau }}_l$ proportional to the distances of the object components. The slow-time dimension with $p\in [0,P-1]$ corresponds to the Doppler frequencies ${\tilde{\nu }}_l$ proportional to object velocities. Sampling leads to the discretized fast-time domain with $t^{\prime }=n{\mathrm{T}}_{\mathrm{s}}$.

  The transmit signal of the $i$-th interferer is given by

  	$x_{\mathrm{I},i}(t)=A_{\mathrm{I},i}{\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I},i}-1}}$	$\overline{x}(t-{\overline{T}}_i-T_{\mathrm{pI},i}{p}_{\mathrm{I}};T_{\mathrm{swI},i})$
  		$\times u(t-T_{\mathrm{pI},i}{p}_{\mathrm{I}};f_{\mathrm{I},i},k_{\mathrm{I},i})$		(6)

  where ${\overline{T}}_i$ is some time offset between victim and interferer transmit sequences and other parameters are analogous to (2). Each interferer transmit signal propagates through a multipath channel. Within our formulation, considerations for this interference channel are analogous to the object channel. However, due to non-coherent demodulation the interferer transmit chirp sequence does not vanish, and the effect of the AAF cannot be neglected. The model is represented in Fig. 1b starting at the triangle symbol, and a time-frequency schematic is found in Fig. 1d. The result reads

  	$r_{\mathrm{I},i}(t^{\prime },p)$	$\approx {\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I},i}-1}}{\overline{u}}_{\mathrm{pI},i}(t^{\prime },p)H_{\mathrm{I},i}({\overline{f}}_{\mathrm{I},i}+k_{\mathrm{I},i}{t}^{\prime },T_{\mathrm{pI},i}{p}_{\mathrm{I}})$
  		$\times \overline{G}(\overline{\mathrm{\Delta }}{f}_{0,i}+\mathrm{\Delta }{k}_i{t}^{\prime })$		(7)

  where ${\overline{u}}_{\mathrm{pI},i}(t^{\prime },p)$ is the noncoherently demodulated interferer transmit chirp, and $H_{\mathrm{I},i}(f,t)$ and $\overline{G}(f)$ denote the interference channel and a modified AAF transfer function, respectively. For derivation details, see Appendix A.

  II-B Inference Problem Statement

  The objective is to jointly estimate the object parameters, including the number of objects (see (5)), as well as the interfering signal parameters and the interfering channel (see (II-A)). Using the assumptions given in Appendix A, the generic model in (1) can be rewritten as

  	$r(n,p)$	$={\displaystyle \sum _{l=0}^{\tilde{L}-1}}{\tilde{\alpha }}_l\exp (-\mathrm{j2}\pi {\tilde{\varphi }}_l{\mathrm{T}}_{\mathrm{s}}n)\exp (-\mathrm{j2}\pi {\tilde{\nu }}_l{T}_{\mathrm{p}}p)$
  		$+u({\mathrm{T}}_{\mathrm{s}}n;\mathrm{\Delta }{\tilde{f}}_0^{(p)},\mathrm{\Delta }{\tilde{k}}^{(p)})\overline{G}(\mathrm{\Delta }{\tilde{f}}_0^{(p)}+\mathrm{\Delta }{\tilde{k}}^{(p)}{\mathrm{T}}_{\mathrm{s}}n)$
  		$\times {\displaystyle \sum _{k=0}^{{\tilde{K}}^{(p)}-1}}{\tilde{\beta }}_k^{(p)}\exp (-\mathrm{j2}\pi {\tilde{\vartheta }}_k^{(p)}{\mathrm{T}}_{\mathrm{s}}n)+\eta (n,p).$		(8)

  The object signal described in (5) is simplified, noting that the object channel component beat frequency is equal to $\tilde{\varphi }=k\tilde{\tau }$. The interference signal is recast with separate parameters for each ramp $p$. For each ramp $p$, the delay-only interference channel is assumed to consist of ${\tilde{K}}^{(p)}$ components with respective amplitudes ${\tilde{\beta }}^{(p)}$ and beat frequencies ${\tilde{\vartheta }}^{(p)}$. Note that (8) constitutes a superposition of “source signals”, each represented by a linear combination of a limited number of components, corrupted by AWGN. Each of these source signals is traditionally referred to as a line spectrum. Inferring the underlying model parameters constitutes an instance of a LSE problem. The inference model in accordance to (8) is given as

  $$𝒓=𝚽(𝜻)𝜶+𝑼(𝜽)𝚿(\mathit{\vartheta })𝜷+𝜼$$

  (9)

  where the vector $𝒓\in {ℂ}^{PN\times 1}$ is constructed from stacking the samples of $r(p,n)$ and $p(𝜼)=\text{CN}(𝜼|0,{\lambda }^{-1}𝑰)$. This is an extension of the well-known model from compressed sensing [51] and also described for MCA [50]. $𝚽(𝜻)=[\mathit{\varphi }({𝜻}_0)\mathrm{\cdots }\mathit{\varphi }({𝜻}_{L-1})]$ is an $M\times L$ dictionary matrix with $M\le L$, typically $M\ll L$. Note that the non-linear parameters in $𝜻$, which define the dictionary basis vectors, may be fixed on-grid or considered unknown and adaptively estimated, referred to as grid-less from here on.

  The object signal dictionary $𝚽(𝜻)$ consists of the basis vectors

  $$\begin{array}{c}\mathit{\varphi }({𝜻}_l=[{\phi }_l{\mathrm{T}}_{\mathrm{s}}{\nu }_l{T}_{\mathrm{p}}])=\frac{1}{\sqrt{PN}}(\exp (-\mathrm{j2}\pi {\nu }_k{T}_{\mathrm{p}}𝒑)\hfill \\ \hfill \otimes \exp (-\mathrm{j2}\pi {\phi }_l{\mathrm{T}}_{\mathrm{s}}𝒏))\end{array}$$

  (10)

  where $\otimes$ denotes the Kronecker product. We estimate the normalized beat frequencies ${\phi }_l{\mathrm{T}}_{\mathrm{s}}\in [-1/2,1/2)$ and normalized Doppler frequencies ${\nu }_l{T}_{\mathrm{p}}\in [-1/2,1/2)$.

  As the interference $𝑼(𝜽)𝚿(\mathit{\vartheta })𝜷$ is modeled for every $p$ separately, $𝑼(𝜽)\in {ℂ}^{PN\times PN}$ and $𝚿(\mathit{\vartheta })\in {ℂ}^{PN\times PK}$ with $𝜽=[{𝜽}^0\mathrm{\cdots }{𝜽}^{P-1}]$ and $\mathit{\vartheta }=[{\mathit{\vartheta }}^0\mathrm{\cdots }{\mathit{\vartheta }}^{P-1}]$ are block-diagonal matrices. I.e.,

  $$𝑼(𝜽)=\left[\begin{array}{ccc}𝑼({𝜽}^{(p=0)})& & \mathrm{𝟎}\\ & \mathrm{\ddots }& \\ \mathrm{𝟎}& & 𝑼({𝜽}^{(P-1)})\end{array}\right]$$

  (11)

  with ${𝜽}^{(p)}=[\mathrm{\Delta }{f}_0^{(p)}\mathrm{\Delta }{k}^{(p)}]$ and

  $$\begin{array}{c}𝑼({𝜽}^{(p)})=\mathrm{Diag}(u({\mathrm{T}}_{\mathrm{s}}𝒏;\mathrm{\Delta }{f}_0^{(p)},\mathrm{\Delta }{k}^{(p)})\hfill \\ \hfill \times \overline{G}(\mathrm{\Delta }{f}_0^{(p)}+\mathrm{\Delta }{k}^{(p)}{\mathrm{T}}_{\mathrm{s}}𝒏))\end{array}$$

  (12)

  is an $N\times N$ diagonal matrix that contains the effect of the non-coherent demodulation and the AAF. The channel dictionary $𝚿(\mathit{\vartheta })$ is constructed analogously from blocks $𝚿({\mathit{\vartheta }}^{(p)})$ of size $N\times K$, and on-grid basis vectors $𝝍({\vartheta }_k)=(1/\sqrt{N})\exp (-\mathrm{j2}\pi {\vartheta }_k{\mathrm{T}}_{\mathrm{s}}𝒏)$, where $K\ge N$ is the chosen size of the grid. As the basis vectors are fixed, we do not write out the parameters $\mathit{\vartheta }$ in the sequel.

  The ill-posed problem of estimating the weights $𝜶$ and $𝜷$ and the according parameters $𝜻$ and $𝜽$ may be solvable under the assumption of sparsity, stemming from $\tilde{L}$ and $\tilde{K}$ being small in (8). That is, we have some optimization problem of the form

  $$\begin{array}{c}\widehat{𝜶},\widehat{𝜷},\widehat{𝜻},\widehat{𝜽}=\underset{𝜶,𝜷,𝜻,𝜽}{\arg \min }{\Vert 𝒓-𝚽(𝜻)𝜶-𝑼(𝜽)𝚿𝜷\Vert }^2\hfill \\ \hfill -f(𝜶,𝜷,𝜻,𝜽)\end{array}$$

  (13)

  with $f(𝜶,𝜷,𝜻,𝜽)$ being a sparsity-inducing penalty function. For a single-source model, cost functions of this form and computationally feasible solution algorithms can be derived in a probabilistic framework [60]. In particular, SBL [57] has been developed for LSE. Variational formulations [61] of SBL have been further developed for on-grid [58], as well as grid-less models [59, 62] that enable super-resolution estimation. The SBL-inspired framework is flexible, lending itself to extensions such as structured line spectra [63] and dense multipath channel models [56]. Our proposed algorithm extends the concept to a superposition of line spectra from different source signals.

  II-C Probabilistic Modeling

  Our probabilistic model has the weights $𝜶$ and $𝜷$ and the noise sensitivity $\lambda$ as latent variables. We apply a Gamma-Gaussian hierarchical model introducing the additional variables ${𝜸}_{\alpha }$ and ${𝜸}_{\beta }$, which leads to a sparse estimate [60]. The non-linear parameters $𝜻$ and $𝜽$ are considered as unknown parameters. Hence, the joint PDF reads

  $$\begin{array}{c}p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)\hfill \\ \hfill =p(𝒓|𝜶,𝜷,\lambda ;𝜻,𝜽)p(𝜶|{𝜸}_{\alpha })p({𝜸}_{\alpha })\\ \hfill \times p(𝜷|{𝜸}_{\beta })p({𝜸}_{\beta })p(\lambda )\end{array}$$

  (14)

  The explicit forms of the terms in (14) are

  	$p(𝒓|𝜶,𝜷,\lambda ;𝜻,𝜽)$	$=\text{CN}(𝒓|𝚽(𝜻)𝜶+𝑼(𝜽)𝚿𝜷,{\lambda }^{-1}𝑰)$		(15)
  	$p(𝜶|{𝜸}_{\alpha })p({𝜸}_{\alpha })$	$={\displaystyle \prod _l}\text{CN}({𝜶}_l|0,{\gamma }_{\alpha ,l}^{-1})$
  		$\times \text{Ga}({\gamma }_{\alpha ,l}|a_0,b_0)$		(16)
  	$p(𝜷|{𝜸}_{\beta })p({𝜸}_{\beta })$	$={\displaystyle \prod _p}{\displaystyle \prod _k}\text{CN}({𝜷}_k^{(p)}|0,\gamma _{\beta ,k}^{-1}{}_{}{}^{(p)})$
  		$\times \text{Ga}({\gamma }_{\beta ,k}^{(p)}|c_0,d_0)$		(17)
  	$p(\lambda )$	$=\text{Ga}(\lambda |e_0,f_0).$		(18)

  The likelihood function $p(𝒓|𝜶,𝜷,\lambda ;𝜻,𝜽)$ is described by a Gaussian distribution due to the AWGN assumption in (8). The weights $𝜶$ and $𝜷$ are modeled as conditionally independent zero-mean Gaussian-distributed with individual precisions governed respectively by Gamma-distributed hyper-parameters ${𝜸}_{\alpha }$ and ${𝜸}_{\beta }$ [57]. The Gamma distribution is the conjugate prior for the precision of a Gaussian [64] and is known to promote sparsity [60]. Consequently, the number of components, $\widehat{K}$ and $\widehat{L}$, can be indirectly estimated by inferring the parameters, ${𝜸}_{\alpha }$ and ${𝜸}_{\beta }$. The prior PDF of the noise precision $\lambda$ is also assumed to be Gamma distributed.

  III Inference Algorithm

  III-A Variational Formulation

  Directly solving the high-dimensional non-linear estimation problem of (13) by statistical inference on the probabilistic model of (14) is computationally infeasible. Therefore, we resort to the variational EM solution [61] to iteratively determine all marginal PDFs of the latent variables and maximum likelihood (ML) estimates of the unknown parameters.

  Denoting the inferred proxy posterior distribution of the latent variables as $q(𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda )$, the EM framework iteratively maximizes the functional

  $$\begin{array}{c}ℒ(q(𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ))\hfill \\ \hfill ={⟨\log \frac{p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)}{q(𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda )}⟩}_{q(𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda )}\end{array}$$

  (19)

  termed the evidence lower bound (ELBO).

  For fixed parameter estimates $\widehat{𝜻}$ and $\widehat{𝜽}$, the ELBO is maximized if $q$ is equal to the joint posterior PDF. As this posterior distribution is intractable, we first constrain the form of $q$ and maximize the ELBO under those constraints to obtain a tractable variational approximation. In line with variational SBL [58, 59], we use the structured mean-field approach [61]

  	$q(𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)$	$=q(𝜶){\displaystyle \prod _l}q({\gamma }_{\alpha ,l})$
  		$\times q(𝜷){\displaystyle \prod _k}q({\gamma }_{\beta ,k})q(\lambda ).$		(20)

  I.e., groups of variables are constrained to factorize in the PDF leading to a set of factors $𝒬$, $q={\prod }_{𝒬}{q}_i$, with $q_i$ being the i-th proxy distribution. To be noted is that our model does not factorize elements of the weight vectors $𝜶$ and $𝜷$, hence posterior correlations between these elements are taken into account. However, the proxy PDFs $q(𝜶)q(𝜷)\ne q(𝜶\mathbf{,}𝜷)$, so that different source signals do in fact factorize, which is a simplifying assumption.⁴⁴If this model assumption is not applied, we essentially obtain a single larger concatenated dictionary $[𝚽𝑼(𝜽)𝚿]$ with its appropriate weights and joint covariance matrix. This leads to a similar algorithm, but without distinct object and interference estimation subroutines. It will be evaluated in Section IV-A as a point of comparison to the proposed method, but will not be further discussed in this work.

  It can be shown that the log-distribution of the $i$-th factor $q_i$, when every other variable and parameter is fixed, is computed by

  $\log q_i={⟨\log p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;\widehat{𝜻},\widehat{𝜽})⟩}_{{\overline{q}}_i}+\text{ const.}$

  (21)

  where ${\overline{q}}_i$ is shorthand notation for “proxy PDF of every factor except for the $i$-th”. Hence, we obtain a set of interdependent implicit equations that are solved by iteratively updating the inferred PDFs. As long as the ELBO is increased at every step, the ordering of updates is in principle arbitrary. In practice though, certain update schemes are used to improve the convergence properties of the algorithm. In the proposed algorithm, we only directly use (21) for the noise precision $q_{\lambda }$. The proxies of the weights $q_{\alpha }$ and $q_{\beta }$ and the according unknown parameters $𝜻$ and $𝜽$ are updated jointly [62], which is described in Section III-B. For the weight precision hyper-parameter distributions $q_{\gamma }$ the fast update scheme [58] is used, discussed in Section III-C. The complete resulting scheme is sketched as pseudo-code in Algorithm 1, making use of Algorithms 2-4 as subroutines. The details are discussed in the sequel.

  III-B Estimation of Unknown Parameters

  In the typical variational EM approach, the updates of unknown parameters are derived with the estimates of the proxy PDFs fixed [61]. However, in our application estimates of the object channel parameters $\widehat{𝜻}$ are strongly tied to the inferred weight proxy PDF $q_{\alpha }$. Similarly, the estimation of interference chirp parameters ${\widehat{𝜽}}^{(p)}$ at ramp $p$ strongly depends on the respective channel weights PDF $q_{\beta }^{(p)}$. Term-wise optimization of the ELBO therefore leads to slower convergence as well as the algorithm being more prone to reach local optima. To mitigate this, the ELBO is maximized jointly for $\widehat{𝜻}$ and $q_{\alpha }$ as well as ${\widehat{𝜽}}^{(p)}$ and $q_{\beta }^{(p)}$ for all $p$. Rewriting (19), this yields the same update equations as (21) for the resulting proxy PDFs. For the parameters, however, we obtain

  $$\begin{array}{c}{\widehat{𝜽}}^{(p)}=\underset{{𝜽}^{(p)}}{\arg \max }\hfill \\ \hfill \int \exp {⟨\log p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;\widehat{𝜻},\widehat{𝜽})⟩}_{{\overline{q}}_{\beta }^{(p)}}\mathrm{d}{𝜷}^{(p)}\end{array}$$

  (22)

  for the chirp parameters ${\widehat{𝜽}}^{(p)}$. The result for the object channel dispersion parameter estimates $\widehat{𝜻}$ is analogous, where estimates for the individual components ${\widehat{𝜻}}_l$ can be computed separately [62, 65] for each spectral line $l$. More details on the derivation are found in Appendix B.

  III-C Fast Component Precision Update and Thresholding

  In [58] it has been shown for variational SBL that with the other variables fixed, the estimates of the component precision hyper-parameters “at infinity” can be derived analytically, yielding a test for component convergence. Hence, at every iteration the dictionary may be adaptively pruned of divergent components, or possibly new components may be added. With this, algorithm convergence is much accelerated. Crucially, it can be seen that for our proposed extended model, the update equations for the factors corresponding to each source signal are identical to the single-source case except with $𝒓$ exchanged by the respective “residuals” ${\widehat{𝒓}}_{\alpha }$ and ${\widehat{𝒓}}_{\beta }$. This means that the formulations of [58, 59] can be followed; Appendix B contains more details.

  Using the estimates of the object signal parameters $\widehat{𝚽}$, $\widehat{𝜶}$, and ${\widehat{𝜸}}_{\alpha }$, the update equation, the pruning of existing components, and the addition of new potential components (see Algorithm 4) are based on

  	${\rho }_l$	$={(\widehat{\lambda }{\widehat{\mathit{\varphi }}}_l^{\mathrm{H}}{\widehat{\mathit{\varphi }}}_l-{\widehat{\lambda }}^2{\widehat{\mathit{\varphi }}}_l^{\mathrm{H}}{\widehat{𝚽}}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{\widehat{𝚽}}_{\overline{l}}^{\mathrm{H}}{\widehat{\mathit{\varphi }}}_l)}^{-1}$		(23)
  	${\omega }_l^2$	$={|(\widehat{\lambda }{\rho }_l{\widehat{\mathit{\varphi }}}_l^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }-{\widehat{\lambda }}^2{\rho }_l{\widehat{\mathit{\varphi }}}_l^{\mathrm{H}}{\widehat{𝚽}}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{\widehat{𝚽}}_{\overline{l}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha })|}^2.$		(24)

  The notation of the subscript $l$ means taking the $l$-th element or column corresponding to that component, and $\overline{l}$ denotes the vector or matrix with said element removed or computed as if it were removed. Existing components or new components are kept if and only if ${\rho }_l/{\omega }_l^2>T$ with $T=1$, where $T$ denotes a threshold and ${\omega }_l^2/{\rho }_l-1$ can be interpreted as an estimated component signal-to-noise ratio (SNR). Note that the original convergence test threshold given by $T=1$ leads to positive bias in the number of estimated components [58, 56]. In [66, 56], the authors propose an adaptation of the threshold $T$ for grid-less SBL-type methods based on extreme value analysis.

  Similar update equations are found for the interference parameters ${\widehat{𝜽}}^{(p)}$, ${\widehat{𝚿}}^{(p)}$, ${\widehat{𝜷}}^{(p)}$, and ${\widehat{𝜸}}_{\beta }^{(p)}$ for all $p$ (see Algorithm 2).

  III-D Algorithm Implementation

  We have now derived structured mean-field variational EM update steps for inference of our model. As mentioned, these steps are interdependent and therefore to be applied iteratively. Explicit forms of the update equations are listed in Appendix B. The algorithm converges to an optimum of the ELBO corresponding to an approximate solution of our signal separation problem. However, it does not necessarily converge to the global optimum, nor is any rate of convergence guaranteed. Our signal separation problem is inherently ambiguous, with many possible combinations of estimates reconstructing the received signal. In addition to the specific update schemes presented previously, the overall scheduling of steps as well as the initialization significantly influence the behavior of the algorithm.

  Algorithm 1 Variational Signal Separation (Main)

  Received signal $𝒓$; initial estimates $\widehat{𝚽}$, $\widehat{𝜶}$, ${\widehat{𝑪}}_{\alpha }$, ${\widehat{𝜸}}_{\alpha }$, $\widehat{𝚿}$, $\widehat{𝜷}$, ${\widehat{𝑪}}_{\beta }$, ${\widehat{𝜸}}_{\beta }$, $\widehat{\lambda }$; thresholds $T_{\alpha }$, $T_{\beta }$

  repeat$▷$ Main iteration

  $▷$ Update interference estimate:

       for all $p$ do $▷$ per-ramp update

           ${\widehat{𝜽}}^{(p)}\leftarrow$ from (66)

           run Algorithm 2(${\widehat{𝒓}}_{\beta }^{(p)}$,${𝑼}^{(p)}(\widehat{𝜽})$,${\widehat{𝚿}}^{(p)}$,${\widehat{𝑪}}_{\beta }^{(p)}$,${\widehat{𝜸}}_{\beta }^{(p)}$,$T_{\beta }$)

       end for

       $\widehat{\lambda }\leftarrow$ from (60) $▷$ update noise precision estimate

       ${\widehat{𝒓}}_{\alpha }=𝒓-𝑼(\widehat{𝜽})\widehat{𝚿}\widehat{𝜷}$ $▷$ update object signal residual

       ${\widehat{𝑪}}_{\alpha }\leftarrow$ from (54) $▷$ update with new residual

       $\widehat{𝜶}\leftarrow$ from (53)

  $▷$ Update object signal estimate:

       run Algorithm 3(${\widehat{𝒓}}_{\alpha }$,$\widehat{𝚽}$,${\widehat{𝑪}}_{\alpha }$,${\widehat{𝜸}}_{\alpha }$,$T_{\alpha }$)

       $\widehat{\lambda }\leftarrow$ from (60)

       ${\widehat{𝒓}}_{\beta }=𝒓-\widehat{𝚽}\widehat{𝜶}$

       ${\widehat{𝑪}}_{\beta }\leftarrow$ from (57)

       $\widehat{𝜷}\leftarrow$ from (56)

  until convergence $▷$ criterion or fixed number

  For notational simplicity we have consistently omitted iteration counters from Algorithm 1. Within the sequence of update steps, it uses the most recent values of other estimates. The main iteration consists of two separated subroutines for the estimations of the respective object and interference signals, similar in structure and connected by computations of ${\widehat{𝒓}}_{\alpha }$, ${\widehat{𝒓}}_{\beta }$, as well as $\widehat{\lambda }$. I.e., the estimated residuals are used at every iteration as inputs to the respective interference and object signal estimation subroutines. This scheme is obtained as the interference and object signals are additive and their components factorize according to (14) and (20). Furthermore, as interference is modeled independently over the ramps, separate steps over the index $p$ are applied.

  Algorithm 2 Interference Signal Estimation (Subroutine)

  Estimates $\widehat{𝒓}$, $\widehat{𝜽}$, $\widehat{𝚿}$, $\widehat{𝑪}$, $\widehat{\lambda }$, $\widehat{𝜸}$; threshold $T$

  $\widehat{k}\leftarrow {\arg \max }_{\overline{k}}(\omega _{}^2{}_{\overline{k}}{}^{}/{\rho }_{\overline{k}})$ where $\overline{k}\in \overline{𝒦}$

  $▷$ find possible new component

  run Algorithm 4($\widehat{𝒓}$,$𝑼(\widehat{𝜽})\widehat{𝚿}$,$\widehat{𝑪}$,$\widehat{𝜸}$,$\widehat{k}$,$𝑼(\widehat{𝜽}){\widehat{𝝍}}_{\widehat{k}}$,$T$)

  $▷$ add component if above threshold

  for all $k\in 𝒦$ do $▷$ update all existing components

       run Algorithm 4($\widehat{𝒓}$,$𝑼(\widehat{𝜽})\widehat{𝚿}$,$\widehat{𝑪}$,$\widehat{𝜸}$,$k$,${\widehat{𝝍}}_k$,$T$)

  end for

  $\widehat{𝜷}\leftarrow$ from (56)

  $\widehat{𝜽}\leftarrow$ from (66) $▷$ update chirp parameter estimates

  Algorithm 3 Object Signal Estimation (Subroutine)

  Estimates $\widehat{𝒓}$, $\widehat{𝚽}$, $\widehat{𝑪}$, $\widehat{\lambda }$, $\widehat{𝜸}$; threshold $T$

  $\widehat{l}\leftarrow L$ $▷$ possible new component index

  ${\widehat{𝜻}}_{\widehat{l}}\leftarrow$ from (65)

  ${\widehat{\mathit{\varphi }}}_{\widehat{l}}\leftarrow$ generate according to signal model

  run Algorithm 4($\widehat{𝒓}$,$\widehat{𝚽}$,$\widehat{𝑪}$,$\widehat{𝜸}$,$\widehat{l}$,${\widehat{\mathit{\varphi }}}_{\widehat{l}}$,$T$)

  $▷$ add component if above threshold

  $\widehat{𝜶}\leftarrow$ from (53)

  for all $l$ do $▷$ update all existing components

       ${\widehat{𝜻}}_l\leftarrow$ from (65)

       ${\widehat{\mathit{\varphi }}}_l\leftarrow$ update with new ${\widehat{𝜻}}_l$

       run Algorithm 4($\widehat{𝒓}$,$\widehat{𝚽}$,$\widehat{𝑪}$,$\widehat{𝜸}$,$l$,${\widehat{\mathit{\varphi }}}_l$,$T$)

       $\widehat{𝜶}\leftarrow$ from (53)

  end for

  Algorithm 4 Fast Component Precision Update (Subroutine)

  Estimates $\widehat{𝒓}$, $\widehat{𝚽}$, $\widehat{𝑪}$, $\widehat{\lambda }$, $\widehat{𝜸}$; existing or new component index $l$, basis vector ${\widehat{\mathit{\varphi }}}_l$, threshold $T$

  ${\rho }_l\leftarrow$ from (23)

  ${\omega }_l^2\leftarrow$ from (24)

  if ${\omega }_l^2/{\rho }_l>T$ then

       ${\widehat{\gamma }}_l\leftarrow (1/({\omega }_l^2-{\zeta }_l))$

       if $l\notin ℒ$ then $▷$ add new component

           $\widehat{𝜸}\leftarrow (\widehat{𝜸}\cup {\widehat{\gamma }}_l)$

           $\widehat{𝚽}\leftarrow (\widehat{𝚽}\cup {\widehat{\mathit{\varphi }}}_l)$

           $L\leftarrow L+1$

       end if

  else

       if $l\in ℒ$ then $▷$ prune existing component

           $\widehat{𝜸}\leftarrow (\widehat{𝜸}\setminus {\widehat{\gamma }}_l)$

           $\widehat{𝚽}\leftarrow (\widehat{𝚽}\setminus {\widehat{\mathit{\varphi }}}_l)$

           $L\leftarrow L-1$

       end if

  end if

  $\widehat{𝑪}\leftarrow$ from (54) using updated $\widehat{𝜸}$ and $\widehat{𝚽}$

  The estimation subroutines of Algorithms 2 and 3 are based on the repeated addition and then update of components by applying (22) and Algorithm 4. The routine adds components it can “find” given the selected threshold, which are then refined or pruned using the now updated estimate of the weight covariance matrix. Hence, the proposed algorithm structure envisions a bottom-up initialization scheme, i.e., that the dictionaries $\widehat{𝚽}$ and analogously $\widehat{𝚿}$ are initialized as empty, which is computationally advantageous. The noise precision is initialized as $\widehat{\lambda }=(2(NP-1))/{\Vert 𝒓\Vert }^2$, i.e., half the received signal power is assumed to be from noise initially.

  Algorithm 1 starts by estimating the interference as described by Algorithm 2. The interference channel is estimated on a grid. Hence, the possible “complete” set of basis vectors of the dictionary matrix is fixed, denoted by the index sets $𝒦$ and $\overline{𝒦}=\{[0,K-1]\setminus 𝒦\}$ for the respective active and passive sets of basis vectors. At every iteration, we select the passive component with the largest estimated component SNR to potentially add to the active dictionary. Irrespective whether or not a new component was added, we proceed with updating all active components. Only a single component is potentially added per main iteration to avoid $\widehat{𝚿}$ converging before the chirp parameter estimates $\widehat{𝜽}$. $\widehat{𝜽}$ is updated once before and after this dictionary update process. The initial chirp parameter estimates are provided by evaluating the cost function on a coarse grid. Further estimates are obtained by applying any established constrained optimization algorithm, with possible values of $\widehat{𝜽}$ constrained to a reasonable range.

  For the grid-less object signal estimation of Algorithm 3, the basis vectors are generated and refined adaptively according to the respective estimates ${\widehat{𝜻}}_l$. Also here, only at most a single component is added per iteration. In this way, the chance of converging to local optima by wrongly assigning part of the interference to the object estimate is minimized. Finally, note that the value of the component threshold for the coherently received object signal $T_{\alpha }$ can be systematically set using the analysis in [66] to achieve a certain false alarm rate. This is not the case for its interference counterpart $T_{\beta }$ due to the unknown $𝜽$. As the object signal is coherently processed over all ramps, $T_{\alpha }>T_{\beta }$ is reasonable. Empirical testing has shown that a low threshold between $0\mathrm{d}\mathrm{B}$ and $6\mathrm{d}\mathrm{B}$ for the interference is generally desirable.

  IV Analysis and Results

  IV-A Fundamental Object Estimation Performance

  The main task of the radar sensor is to accurately determine the distances and velocities of surrounding objects. Therefore, we evaluate the proposed algorithm based on the quality of object parameter estimates. First, a comparison to the CRLB [67] on the estimation error variance is carried out. In order to simplify the analysis and numerical evaluations, a scenario with a single transmit ramp, and hence delay-only estimation, is set up. The object channel consists of a single line spectral component whereas the interference channel is constant, i.e., there is only direct-path interference corresponding to a single component at excess delay zero. The interference signal template is as described in Section II and Appendix A. The AAF frequency transfer function is a raised cosine with Nyquist bandwidth $f_{\mathrm{s}}/4$ and roll-off factor 0.25. Relevant signal parameters for this experiment are in Table I (column denoted Simulation I).

  TABLE I: Summary of relevant signal parameters for the analysis. ${\mathrm{c}}_0$ denotes the speed of light in a vacuum.

  Parameter	Simulation I	Simulation II	Measurement I	Measurement II
  Victim Radar and Object Channel
  Ramp start frequency $f_0(\mathrm{GHz})$	$79$	$79$	$79$	$76.2$
  Ramp slope $k(\mathrm{GHz}/\mathrm{s})$	${10}^4$	${10}^4$	$2.08\times {10}^4$	$2.08\times {10}^4$
  Ramp duration $T_{\mathrm{sw}}(\mu \mathrm{s})$	$25$	$25$	$6.39$	$24$
  Pulse duration $T_{\mathrm{p}}(\mu \mathrm{s})$	$25$	$25$	$44$	$44$
  No. of ramps $P$	1	16	32	32
  No. of samples $N$	256	128	256	256
  Sampling frequency ${\mathrm{f}}_{\mathrm{s}}=1/{\mathrm{T}}_{\mathrm{s}}(\mathrm{MHz})$	10.2	5.1	40	40
  No. of components $\tilde{L}$	1	10	unknown	unknown
  Component weights $\tilde{\alpha }$	$1\exp (\mathrm{j}{\phi }_{\alpha })$, ${\phi }_{\alpha }\sim 𝒰(0,2\pi )$	$\log {|\alpha |}_{[\mathrm{dB}]}^2$ $=-40\log (\tilde{\tau }{\mathrm{c}}_0+1)+x$, $x\sim 𝒰(-3,3)$	unknown	unknown
  Component delays $\tilde{\tau }(\mathrm{ns})$	$80.06$	$𝒰(3.97,127)$	unknown	unknown
  Component Doppler frequencies $\tilde{\nu }(\mathrm{kHz})$	$0$	$𝒰(-5,5)$	unknown	unknown
  Interferer and Interference Channel
  Ramp start frequency $f_{\mathrm{I}}(\mathrm{GHz})$	$79.01$	$\{79.002,79.004,79.008\}$	-	$76.1$
  Ramp slope $k_{\mathrm{I}}(\mathrm{GHz}/\mathrm{s})$	$9.2\times {10}^3$	$\{9.8321,9.7122,9.3925\}\times {10}^3$	-	$2.6923\times {10}^4$
  Ramp duration $T_{\mathrm{swI}}(\mu \mathrm{s})$	$25$	$25.02$	-	$26$
  Pulse duration $T_{\mathrm{pI}}(\mu \mathrm{s})$	$25$	$\{75.01,50.01,25.02\}$	-	$25$
  No. of ramps $P_{\mathrm{I}}$	1	$\{2,4,8\}$	-	$8$
  No. of components $\tilde{K}$	1	10	unknown	unknown
  Component weights $\tilde{\beta }$	$|\beta |\exp (\mathrm{j}{\phi }_{\beta })$, $|\beta |$ according to $\mathrm{SIR}$, ${\phi }_{\beta }\sim 𝒰(0,2\pi )$	$\log {|\beta |}_{[\mathrm{dB}]}^2$ $=-20\log ({\tilde{\tau }}_{\mathrm{I}}{\mathrm{c}}_0+1)+x$, $x\sim 𝒰(-10,0)$	-	unknown
  Component delays ${\tilde{\tau }}_{\mathrm{I}}(\mathrm{ns})$	$0$	0 for first component, $𝒰(3.97,127)$ for others	-	unknown

  We study the expected root mean squared error (RMSE) of the normalized beat frequency estimate $\sqrt{⟨{\Vert \widehat{\phi }-\tilde{\phi }\Vert }^2⟩}$ over $500$ simulated realizations. The expectation is replaced by the sample average, where only realizations with correct model order estimates $\widehat{L}=\tilde{L}=1$ are taken into account. Different realizations consist of varying the noise as well as the phases of the complex scalar weights $\tilde{\alpha }$ and $\tilde{\beta }$. The RMSE is analyzed for different SNRs and signal-to-interference ratios, defined as

  	$\mathrm{SNR}$	$=\lambda {\Vert 𝚽(\widetilde{𝜻})\widetilde{𝜶}\Vert }^2$		(25)
  	$\mathrm{SIR}$	$={\Vert 𝚽(\widetilde{𝜻})\widetilde{𝜶}\Vert }^2/{\Vert 𝑼(\widetilde{𝜽})\widetilde{𝚿}\widetilde{𝜷}\Vert }^2.$		(26)

  Three algorithm variants are considered. First, mutual interference is completely neglected in the inference model, i.e., Algorithm 2 and associated steps are removed. Second, the factorized algorithm as proposed in this work can be applied. Third, as mentioned in Section III, it is possible to jointly model the interference and object channels within a larger concatenated dictionary by modeling the joint weight proxy $q(𝜶,𝜷)$ in (20). For all variants, component acceptance thresholds are set to $T_{\alpha }=9\mathrm{d}\mathrm{B}$ and $T_{\beta }=3\mathrm{d}\mathrm{B}$. Results are compared to the root mean CRLB whose values are computed from the given signal model of (9). Fig. 2 shows an example result of the proposed algorithm.

  

  Figure 2: Example signal and object estimate at $\mathrm{SNR}=30\mathrm{d}\mathrm{B}$ and $\mathrm{SIR}=-15\mathrm{d}\mathrm{B}$. To reduce visual clutter, only certain signals are shown. The upper plot contains the fast-time view of the complete received signal $𝒓$, the object source signal ${𝒓}_{\mathrm{O}}$ and its estimate ${\widehat{𝒓}}_{\mathrm{O}}$. The lower plot contains the same signals in the delay domain, as well as the true object channel components ($\widetilde{𝜻}$ and $\widetilde{𝜶}$), their estimates ($\widehat{𝜻}$ and $\widehat{𝜶}$), and the estimated interference channel components ($\widehat{\mathit{\vartheta }}$ and $\widehat{𝜷}$).

  

  Figure 3: Statistical results of single-object delay estimation RMSE. (a) Evaluated over varying SNR at a fixed SIR of $0\mathrm{d}\mathrm{B}$. (b) Evaluated over varying SIR at a fixed SNR of $30\mathrm{d}\mathrm{B}$. (c) Evaluated for the AAF model error scenario.

  In Fig. 3a it can be seen for a fixed $\mathrm{SIR}=0\mathrm{d}\mathrm{B}$ that the algorithm jointly modeling object and interference channels achieves the CRLB once the SNR is high enough for the object channel component to be detected. The factorized algorithm only experiences a slight performance loss visible at high SNRs, whereas the “object only” algorithm, as expected, cannot achieve an estimation error below a certain value determined by the disturbing interference.

  We also observe the performance over varying SIR values, plotted in Fig. 3b for a fixed $\mathrm{SNR}=30\mathrm{d}\mathrm{B}$. The performance of object-only channel estimation strongly deteriorates with increasing interference power. Although the RMSE seemingly recovers at lower SIR, this is misleading as the estimate gets dominated by the interference signal and the number of missed detections is growing. Below around $-30\mathrm{d}\mathrm{B}$, the component cannot be detected at all. This is not the case for the signal separation-based methods. The proposed factorized algorithm shows its largest evaluated deviation from the theoretical optimum at $-15\mathrm{d}\mathrm{B}$ (highlighted by a close-up). We observe that the worst case for the factorized algorithm is related to the condition

  $$\Vert {\widetilde{𝚽}}^{\mathrm{H}}\widetilde{𝚽}\widetilde{𝜶}\Vert \approx \Vert {\widetilde{𝚽}}^{\mathrm{H}}𝑼(\widetilde{𝜽})\widetilde{𝚿}\widetilde{𝜷}\Vert$$

  (27)

  measuring the similarity between object and interference signals when projected onto the object signal base $\widetilde{𝚽}$, with equality reached around $-20<\mathrm{SIR}<-15\mathrm{d}\mathrm{B}$. We hypothesize that such a condition has theoretical significance for the signal separation problem, but further rigorous analysis is necessary. The algorithm again tends to the CRLB at lower SIR. From a practical point of view, it is certainly most important in a radar application to achieve some guaranteed performance at critical low-SIR and low-SNR scenarios.

  IV-B Robustness to Interference Model Error

  Any model-based algorithm breaks down if the inference model is not sufficiently accurate to reality. In our application, it cannot for example be assumed that we have perfect knowledge of the AAF transfer function $G(f)$ for specific radar hardware. Hence, we investigate the behavior of the algorithm when assuming, as previously, the raised cosine envelope for inference but in fact generating the tested interference signal using convolution with a Butterworth lowpass filter with equivalent filter parameters. This constitutes a small, but not negligible error in the filter’s assumed transfer function. An example result is shown in Fig. 4.

  

  Figure 4: Example signal and estimation results at $\mathrm{SNR}=30\mathrm{d}\mathrm{B}$ and $\mathrm{SIR}=-15\mathrm{d}\mathrm{B}$ with model error. The included signals are the same as for Fig. 2. To note in particular is that the estimate for the interference channel is biased, meaning that it estimates more than a single component. This compensates for the model error, resulting in a more robust estimate of the object component.

  Fig. 3c shows the same analysis as for Fig. 3b, but with the introduced model error. As expected, an increasing error in the object channel parameter estimate becomes visible below a certain $\mathrm{SIR}$. Nevertheless, the algorithm does not abruptly break down and behaves robustly, reaching the CRLB for $\mathrm{SIR}\ge -30\mathrm{d}\mathrm{B}$. The reason for this adaptability is illustrated in Fig. 4. Note in Fig. 4 that several line spectral components are estimated for the interference channel, even though the signal is simulated with line-of-sight interference only. These additional components compensate for the AAF model error.

  IV-C Statistical Performance Comparison

  

  Figure 5: Statistical results of multi-object simulation scenario. (a)-(c) Evaluated over varying SNR at a fixed SIR of $-20\mathrm{d}\mathrm{B}$. (d)-(f) Evaluated over varying SIR at a fixed SNR of $60\mathrm{d}\mathrm{B}$.

  

  Figure 6: Example signals and estimates of a single interfered ramp at $\mathrm{SNR}=60\mathrm{d}\mathrm{B}$ and $\mathrm{SIR}=-15\mathrm{d}\mathrm{B}$. The included signals are the same as for Fig. 2.

  We study the performance of the proposed algorithm in broader, more complex scenarios. Table I shows the relevant simulation parameters in the column denoted Simulation II. The radar system receives $N=128$ fast-time samples from $P=16$ ramps each, with the task to estimate objects’ range and velocity parameters. The interferer transmit ramp parameters are chosen randomly for each realization from three different sets representing three scenarios. The three sets represent increasing interference burst durations while decreasing the number of interfered ramps, which arises naturally due to the form of FMCW transmit signals.

  Both object and interference channels consist of $\tilde{L}=\tilde{K}=10$ line spectral components. Dispersion parameters are chosen randomly on a uniform distribution, while component magnitudes are set proportionately according to the radar equation [68] with some random variation to account for the radar cross section of different physical objects. Note also the different exponents in the radar equation for victim and interference channels [69]. For the interference, the first channel component is fixed to a delay of zero and relative magnitude of one. This is the direct path interference between sensors, which has been observed to dominate in power [70].

  As these are multi-target scenarios, performance metrics appropriate for multi-target detection are selected. The generalized optimal subpattern assignment (GOSPA) metric [71] is suited to such cases. Its evaluation starts with an assignment of the estimated components to the true ones based on an upper-bounded estimation error metric. GOSPA can be split into a so-called localization and a cardinality error term. The former is the sum of estimation errors for the assigned components, while the latter penalizes non-assigned detections. The cardinality error can be further split into two terms proportional to the number of misdetections and false alarms, respectively. We use the assignment method of GOSPA and analyze variants of these three error terms, leading to an insightful set of performance metrics. For the localization error, we change the sum of errors of [71] to the mean of Euclidean distances. This yields a metric comparable to the single-component RMSE, here termed the mean assigned RMSE. For the two cardinality error terms, we evaluate the absolute number of misdetections and false alarms directly.

  Our proposed algorithm is compared with two other methods for interference mitigation found in literature. Note that these are preprocessing methods to be employed prior to a separate object detection step in the automotive radar signal processing chain. In order to ensure a fair comparison, the “object only” algorithm of Section IV-A is applied after preprocessing, as well as in the No Mitigation and No Interference test cases. The two methods, both introduced in Section I, are:

  • •

    Zeroing - zeroing relies on the explicit prior detection of interfered samples, which is assumed perfect here. Comparison of different algorithms for the detection of interfered samples is out of the scope of this work.

  • •
    \Acf

    mca - using the concept and model of [49]. Our implementation follows that of [72, Algorithm 30].

  Fig. 5 illustrates the mean of the resulting metrics computed from 500 realizations. Over SNR, interference leads to a lower limit on the achievable estimation error and the number of misdetections and causes false alarms. With preprocessing, the estimation error and number of misdetections improve, but the performance is still limited and improvements come at the cost of false alarms. This behavior can be explained by considering the effects of preprocessing on the form of the resulting signal. Zeroing merely exploits the time-limited nature of the interference, and hence leads to a signal with missing clusters of samples. MCA considers a signal separation problem, but does not employ a parametric model for the interference bursts and only applies single-ramp on-grid LSE for the object signal. Preprocessing therefore typically causes a distortion of the object signal that is a model error for the subsequent detector, primarily leading to false alarms. Our proposed method mitigates this behavior and reaches performance very close to the No Interference case. The estimation error and the number of misdetections decrease with increasing SNR, and the number of false alarms is constant according to the set detection threshold.

  The bottom row of Fig. 5 shows the results at a high SNR, over SIR levels ranging from highly interfered to weakly interfered object signals. Ideally, perfect interference cancellation leads to values identical to the No Interference case at $\mathrm{SNR}=60\mathrm{d}\mathrm{B}$. The proposed method approximately achieves this up until $\mathrm{SIR}\approx 0\mathrm{d}\mathrm{B}$. Its performance then slightly decreases with decreasing interference power up until $\mathrm{SIR}\approx 15\mathrm{d}\mathrm{B}$, after which it slowly recovers. This is the same effect as discussed in Section IV-A. Note that unlike for the simulation scenario of Fig. 3, in the scenario of Fig. 5 the condition of (27) is met roughly at $5\mathrm{d}\mathrm{B}<\mathrm{SIR}<10\mathrm{d}\mathrm{B}$. This strengthens our hypothesis on the significance of this condition. For the other evaluated methods, we see a bounded improvement in RMSE and misdetection rate at the cost of further false alarms. Naturally, perfect zeroing is independent of interference power while MCA performance is generally improved with increasing SIR.

  We also include in Fig. 6 an example plot analogous to Figs. 2 and 4. We typically observe $\widehat{K}>\tilde{K}$, i.e., the amount of interference channel line spectral components is overestimated in our interference estimation subroutine. This necessarily occurs when there are several relatively closely-spaced components, as we employ on-grid estimation on a small grid $(K=2N)$ with a relatively low detection threshold $T_{\beta }$. Nevertheless, it can be seen that the object signal is successfully estimated.

  In summary, our algorithm statistically outperforms the investigated preprocessing methods across the whole investigated SNR/SIR space and only slightly diverges from optimal interference cancellation in a less critical high-SIR region of this exploration space.

  IV-D Measurement Example

  

  Figure 7: Example of algorithm applied to real automotive radar measurement data. (a) Photo of the measurement scene. (b)-(c) Range-Doppler spectrum of the measured signal without/with interference and corresponding object detections. (d) Range-Doppler spectrum of the interference-cancelled signal and corresponding object detections as a result of the proposed algorithm.

  In order to illustrate the practical application of the proposed method, we apply it to radar measurement data from a real traffic scenario. Details on the measurement campaign can be found in previous works [73, 74]. Measurement parameters are summarized in Table I (Measurement I). Note that there is no ground truth of the channel parameters available, as the measurements were done with a single radar sensor in completely uncontrolled scenarios. We select a single measurement and showcase the results qualitatively. Furthermore, the measurements do not inherently contain any interference. In prior work, simulated interference was added to the measured signals. In this work however, measured interference is added, as described in Appendix C, with parameters listed in Table I (Measurement II). The interference was obtained from a separate measurement, due to the difficulty of generating and measuring interference separately in the same inner city traffic environment. Nevertheless, the applied scheme enables a qualitative comparison of our mitigation result to an interference-free one, with both object and interference components obtained using actual radar sensors.

  Fig. 7a shows a photo of the selected Measurement I scene. Particularly, visible on the photo is the car in front and to the left of the ego-radar, as well as the facade of a row of buildings to the right. The interference-free spectrum in Fig. 7b hence contains a few close object components with a positive relative velocity, i.e., the car moving away from the sensor. The buildings lead to a large number of components distributed over distances, with negative relative velocities with respect to the measurement vehicle. Other detected components are reflections not easily associated with the camera image or false alarms due to noise. When strong interference is introduced, the range-Doppler spectrum becomes dominated by it, see Fig. 7c. Directly applying object detection onto this signal leads to no discernibly useful result. Finally, the spectrum of the measurement when subtracting the interference estimated by the proposed processing algorithm, and the estimated object components are shown in Fig. 7d. The resulting spectrum is visually very close to the interference-free one and the detected components correspond well to the significant parts of the described scene. This very promising result is nevertheless only preliminary, and a larger scale analysis with better known hardware and measurement environment is a necessary part of future work.

  We also qualitatively illustrate the fast-time signal of a single ramp in Fig. 8 as an example. The top graph shows a large interference burst not too dissimilar to, but distinct from the simulated ones of Figs. 2 and 4. The bottom plot provides a close-up view of the signal after subtraction of the interference estimate as well as shows the object signal estimate.

  

  Figure 8: Qualitative example of the processed signal at a single interfered ramp and object signal estimation results. It shows a ramp of the measurement $𝒓$, of the mitigated signal being the object residual ${\widehat{𝒓}}_{\alpha }$, and of the object signal estimate ${\widehat{𝒓}}_{\mathrm{O}}=\widehat{𝚽}\widehat{𝜶}$. The top and bottom graphs are the same result with different y-axis limits, the latter showing the much lower amplitude object signal.

  V Conclusion

  This paper presents an algorithm to mitigate the mutual interference of automotive FMCW radars. We fully consider the underlying signal model and systematically design a model-based inference algorithm. Our description of the interference considers both sensor-based effects and the multipath propagation environment. We propose to infer the delay-Doppler object channel and ramp-wise delay-only interference channels. The result is a superposition of LSE problems, with the non-coherent interference chirp envelopes modeled by additional unknown parameters. We extend the state-of-the-art sparse probabilistic inference approach to such a superposition within the variational EM framework. Our proposed inference model leads to an iterative algorithm consisting of subsequent object and interference estimation subroutines. This results in robust object detection performance that is often comparable to the interference-free case. Conventional interference mitigation preprocessing is shown to often lead to increasing false alarms in the object detection step, which the proposed algorithm is not susceptible to, in comparison. It is also shown that considering the estimation of interference as an LSE problem can offset minor model errors.

  Specific choices in the proposed algorithm design, initialization and scheduling are essential for its performance. This includes employing grid-less estimation for the object channel, yielding super-resolution accuracy, while using an on-grid algorithm for the interference channels. Parameter optimization and inference of the respective channel weights are derived from a joint update, which decreases the probability of converging to a local optimum of the inherently highly ambiguous signal separation problem.

  The proposed method has some limitations due to assumptions made in its derivation, as elaborated in the paper. Most prominently, at most a single interference chirp is assumed per ramp, and interference chirps are present over the full passband of the receiver’s AAF. While ways to adapt the algorithm are briefly discussed, more work is needed to verify their practical viability in different scenarios. Furthermore, we posit that signal separability is fundamentally limited by how distinct the involved signal bases are. Investigating such limits is also an interesting prospect for further work. Finally, note that although here presented for a specific application, the proposed algorithm is likely applicable to other signal separation problems of this kind.

  Appendix A Derivation of Signal Model

  A-A Object Signal

  To derive (5), we start with the transmit signal of (2). It propagates through the environment and reflects off of (moving) objects, i.e., a weighted sum of delayed transmit signals reach the receiver. The radar echo as a function of time $t$ can be derived for ideal point reflections starting with $y_{\mathrm{O}}(t)={\sum }_{l=0}^{\tilde{L}-1}{\tilde{\alpha }}_{l}x(t-{\tilde{\tau }}_l(t))$, where ${\tilde{\tau }}_l(t)\approx 2(d_l+v_{l}t)/{\mathrm{c}}_0$ with $d_l$ and $v_l$ being radial distance and velocity [75], and ${\mathrm{c}}_0$ denotes the speed of light. In our derivations, we use an analogous formulation based on channel modeling [76]. I.e., the transmit signal propagates through the radar channel which can be described by the operation

  	$y_{\mathrm{O}}(t)$	$={\displaystyle {\int }_{\nu }}{\displaystyle {\int }_{\tau }}{h}_{\mathrm{O}}(\tau ,\nu )x(t-\tau )\exp (\mathrm{j2}\pi \nu t)d\tau d\nu$
  		$={\displaystyle \sum _l}{\tilde{\alpha }}_{l}x(t-{\tilde{\tau }}_l)\exp (\mathrm{j2}\pi {\tilde{\nu }}_{l}t)$		(28)

  with

  $h_{\mathrm{O}}(\tau ,\nu )={\displaystyle \sum _{l=0}^{\tilde{L}-1}}{\tilde{\alpha }}_l\delta (\tau -{\tilde{\tau }}_l)\delta (\nu -{\tilde{\nu }}_l)$

  (29)

  where the channel is described by its spreading function $h_{\mathrm{O}}(\tau ,\nu )$ in delay $\tau$ and Doppler frequency $\nu$. The model assumes negligible coupling between delay and Doppler dispersion domains. Expanding on (28), we have

  	$y_{\mathrm{O}}(t)$	$=A{\displaystyle \sum _{l=0}^{\tilde{L}-1}}{\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p-{\tilde{\tau }}_l;T_{\mathrm{sw}}){\tilde{\alpha }}_l\exp (-\mathrm{j2}\pi {\tilde{\nu }}_{l}t)$
  		$\times \exp (\mathrm{j}(2\pi f_0(t-T_{\mathrm{p}}p-{\tilde{\tau }}_l)+\pi k{(t-T_{\mathrm{p}}p-{\tilde{\tau }}_l)}^2))$		(30)

  where we will assume $A=1$ for simplicity, without loss of generality.

  To simplify this, we first note that in a continuous wave radar system, the maximum effective delay $\max ({\tilde{\tau }}_l)$ of the channel is much shorter than the duration of a transmit ramp $T_{\mathrm{sw}}$. We can hence approximately neglect the spreading of the ramp envelope itself. I.e., we swap the order of summations and apply $\overline{x}(t-T_{\mathrm{p}}p-{\tilde{\tau }}_l;T_{\mathrm{sw}})\approx \overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})$. Furthermore, the stop-and-go approximation ${\tilde{\nu }}_{l}t\approx {\tilde{\nu }}_l{T}_{\mathrm{p}}p$ as commonly known in automotive radar [68] is applied. We obtain

  	$y_{\mathrm{O}}(t)$	$\approx {\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}}){\displaystyle \sum _{l=0}^{\tilde{L}-1}}{\tilde{\alpha }}_l\exp (-\mathrm{j2}\pi {\tilde{\nu }}_l{T}_{\mathrm{p}}p)$
  		$\times \exp (\mathrm{j}(2\pi f_0(t-T_{\mathrm{p}}p-{\tilde{\tau }}_l)+\pi k{(t-T_{\mathrm{p}}p-{\tilde{\tau }}_l)}^2)).$		(31)

  Multiplying out and rearranging the terms in the chirp exponential straightforwardly leads to

  	$y_{\mathrm{O}}(t)$	$={\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})$
  		$\times \exp (\mathrm{j}(2\pi f_0(t-T_{\mathrm{p}}p)+\pi k{(t-T_{\mathrm{p}}p)}^2)))$
  		$\times {\displaystyle \sum _{l=0}^{\tilde{L}-1}}{\tilde{\alpha }}_l\exp (\mathrm{j}(-2\pi (f_0+k(t-T_{\mathrm{p}}p)){\tilde{\tau }}_l)$
  		$\times \exp (\mathrm{j}(\pi k{\tilde{\tau }}_l^2))\exp (-\mathrm{j2}\pi {\tilde{\nu }}_l{T}_{\mathrm{p}}p).$		(32)

  Reasoning again based on the short duration of the channel, we neglect the term $k{\tilde{\tau }}_l^2\approx 0\forall l$. The rest of the expression in the last brackets is then the Fourier transform of the channel spreading function. Defining this as the object channel transfer function $H_{\mathrm{O}}(f,t)$, we write ⁵⁵Somewhat uncommonly for automotive radar literature, in our formulation the fast-time signal measures the frequency domain of the channel whose inverse Fourier spectrum represents the delay dispersion. Note that delay and frequency are always dual domains of the radio channel, whereas delay being proportional to a beat frequency is a specific “quirk” of FMCW radar.

  	$y_{\mathrm{O}}(t)$	$\approx {\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})\exp (\mathrm{j}(2\pi f_0(t-T_{\mathrm{p}}p)$
  		$+\pi k{(t-T_{\mathrm{p}}p)}^2))H{}_{\mathrm{O}}(f_0+k(t-T_{\mathrm{p}}p),T_{\mathrm{p}}p).$		(33)

  Next, demodulation (demixing) according to stretch processing [77] is applied. Demixing is coherent for the object signal, so the transmit ramp envelope $\overline{x}(t)$ over the ramp duration stays constant and the chirp term vanishes. I.e.,

  	$y_{\mathrm{O}}^{\prime }(t)$	$=(x^{\ast }(t){\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})\exp (\mathrm{j}(2\pi f_0(t-T_{\mathrm{p}}p)$
  		$+\pi k{(t-T_{\mathrm{p}}p)}^2))H_{\mathrm{O}}(f_0+k(t-T_{\mathrm{p}}p),T_{\mathrm{p}}p))$
  		$={\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})H_{\mathrm{O}}(f_0+k(t-T_{\mathrm{p}}p),T_{\mathrm{p}}p).$		(34)

  Finally, we assume that the object channel contains no significant amount of energy outside of the constant passband of the AAF, hence the convolution with $g(t)$ can be neglected. The signal is projected onto the fast-time $t^{\prime }\in (0,T_{\mathrm{sw}})=t-T_{\mathrm{p}}p\forall p$ and slow-time $p$ so that the explicit sum over $p$ vanishes as well. This leads to the result in (5).⁶⁶This is a standard result that is well established. Nevertheless, its detailed derivation within our channel modeling-based formulation in this work is warranted to set the derivation of the interference signal into the proper context.

  A-B Interference Signal

  We derive the received interference signal from a single interferer, $r_{\mathrm{I},i}(t^{\prime },p)$ of (II-A), starting from the interferer transmit signal $x_{\mathrm{I},i}(t)$ of (6). Note that for notational simplicity, we drop the explicit interferer subscript $i$ in the sequel. First, propagation through a multipath channel leads to the expression

  	$y_{\mathrm{I}}(t)$	$\approx x_{\mathrm{I}}(t){\displaystyle \sum _k^{\tilde{K}-1}}{\tilde{\beta }}_k\exp (-\mathrm{j2}\pi (f_{\mathrm{I}}+k_{\mathrm{I}}(t-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})){\tilde{\tau }}_{\mathrm{I},k})$
  		$\times \exp (-\mathrm{j2}\pi T_{\mathrm{pI}}{p}_{\mathrm{I}}{\tilde{\nu }}_{\mathrm{I},k})$
  		$=x_{\mathrm{I}}(t)H_{\mathrm{I}}(f_{\mathrm{I}}+k_{\mathrm{I}}(t-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}}),T_{\mathrm{pI}}{p}_{\mathrm{I}})$		(35)

  analogously as for the object signal $y_{\mathrm{O}}(t)$. However, this signal is non-coherently demodulated by multiplication with $x^{\ast }(t)$. To elaborate, we explicitly write out the mixing, i.e.,

  	$y_{\mathrm{I}}^{\prime }(t)$	$=\overline{A}{\displaystyle \sum _{p=0}^{P-1}}\overline{x}(t-T_{\mathrm{p}}p;T_{\mathrm{sw}})$
  		$\times \exp (-\mathrm{j}(2\pi f_0(t-T_{\mathrm{p}}p)+\pi k{(t-T_{\mathrm{p}}p)}^2))$
  		$\times ({\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I}}-1}}\overline{x}(t-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}};T_{\mathrm{swI}})$
  		$\times \exp (\mathrm{j}(2\pi f_{\mathrm{I}}(t-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})$
  		$+\pi k_{\mathrm{I}}{(t-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})}^2))$
  		$\times H_{\mathrm{I}}(f_{\mathrm{I}}+k_{\mathrm{I}}(t-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}}),T_{\mathrm{pI}}{p}_{\mathrm{I}}))$		(36)

  setting $\overline{A}=A_{\mathrm{I}}{A}^{\ast }=1$ without loss of generality in the sequel.

  For convenience, we rewrite this in the fast- and slow-time domain description by substituting $t=t^{\prime }+T_{\mathrm{p}}p$, leading to

  	$y_{\mathrm{I}}^{\prime }(t^{\prime },p)$	$=\exp (-\mathrm{j}(2\pi f_0{t}^{\prime }+\pi kt_{}^{\prime }{}_{}{}^2))$
  		$\times {\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I}}-1}}\overline{x}(t^{\prime }+T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}};T_{\mathrm{swI}})$
  		$\times \exp (\mathrm{j}(2\pi f_{\mathrm{I}}(t^{\prime }+T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})$
  		$+\pi k_{\mathrm{I}}{(t^{\prime }+T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})}^2))$
  		$\times H_{\mathrm{I}}(f_{\mathrm{I}}+k_{\mathrm{I}}(t^{\prime }+T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}}),T_{\mathrm{pI}}{p}_{\mathrm{I}})$		(37)

  noting again that the sum over $p$ vanishes.

  Carrying out the multiplication of the victim and interferer transmit chirps and rearranging the resulting phase terms with respect to $t^{\prime }$-dependency yields

  	$y_{\mathrm{I}}^{\prime }(t^{\prime },p)$	$={\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I}}-1}}\overline{x}(t^{\prime }+T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}};T_{\mathrm{swI}})$
  		$\times \exp (\mathrm{j}(2\pi (f_{\mathrm{I}}+k_{\mathrm{I}}(T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})-f_0)t^{\prime }$
  		$+\pi (k_{\mathrm{I}}-k)t_{}^{\prime }{}_{}{}^2))$
  		$\times \exp (\mathrm{j}(2\pi f_{\mathrm{I}}(T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})$
  		$+\pi k_{\mathrm{I}}{(T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}})}^2))$
  		$\times H_{\mathrm{I}}(f_{\mathrm{I}}+k_{\mathrm{I}}(t^{\prime }+T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}}),T_{\mathrm{pI}}{p}_{\mathrm{I}}).$		(38)

  Due to non-coherent demodulation, the effect of the AAF cannot be neglected for the interference. We can approximate the AAF for every slow-time index $p$ separately as the convolution with an impulse response $g(t^{\prime })$. We define ${\overline{T}}_0=T_{\mathrm{p}}p-\overline{T}-T_{\mathrm{pI}}{p}_{\mathrm{I}}$ to de-clutter notation, as well as the non-coherent demixed chirp

  	${\overline{u}}_{\mathrm{pI}}(t^{\prime },p)$	$=\exp (\mathrm{j}(2\pi \mathrm{\Delta }{f}_0{t}^{\prime }+\pi \mathrm{\Delta }\tilde{k}t_{}^{\prime }{}_{}{}^2))$
  		$\times \exp \left(\mathrm{j}(2\pi f_{\mathrm{I}}{\overline{T}}_0+\pi k_{\mathrm{I}}{\overline{T}}_0^2)\right)$		(39)

  where $\mathrm{\Delta }{f}_0=f_{\mathrm{I}}-f_0$ and $\mathrm{\Delta }\tilde{k}=k_{\mathrm{I}}-k$. This leads to the convolution integral expression

  	$r_{\mathrm{I}}(t^{\prime },p)$	$={\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I}}-1}}{\overline{u}}_{\mathrm{pI}}(t^{\prime },p)({\displaystyle {\int }_{\varsigma }}\overline{x}(t^{\prime }-\varsigma +{\overline{T}}_0;T_{\mathrm{swI}})$
  		$\times \exp \left(-\mathrm{j}(2\pi (\mathrm{\Delta }{f}_0+k_{\mathrm{I}}{\overline{T}}_0)\varsigma +2\pi \mathrm{\Delta }\tilde{k}{t}^{\prime }\varsigma -\pi \mathrm{\Delta }\tilde{k}{\varsigma }^2)\right)$
  		$\times H_{\mathrm{I}}(f_{\mathrm{I}}+k_{\mathrm{I}}((t^{\prime }-\varsigma )+{\overline{T}}_0),T_{\mathrm{pI}}{p}_{\mathrm{I}})g(\varsigma )\mathrm{d}\varsigma ).$		(40)

  Incorporating (40) directly into the inference model is possible, but leads to significantly increased complexity. In order to derive a simplified form used in this work, certain further approximations are considered. We introduce here the first two out of four main assumptions employed for inference:

  Assumption 1

  The envelope term $\overline{x}(t^{\prime }-\varsigma +{\overline{T}}_0;T_{\mathrm{swI}})$ can be neglected if it is assumed to be long enough in comparison to the filter impulse response. This assumption is correct in cases where the interferer’s transmit ramp is active during the whole time its frequency course would be inside the victim receiver’s IF bandwidth.

  Assumption 2

  We take the channel term $H_{\mathrm{I}}$ outside the convolution integral, neglecting its dependence on $\varsigma$. This is based on the reasoning that the chirp term ${\overline{u}}_{p\mathrm{I}}$ dominates the result of the convolution.

  We then rearrange the remaining chirp terms and obtain

  	$r_{\mathrm{I}}(t^{\prime },p)$	$\approx {\displaystyle \sum _{p_{\mathrm{I}}=0}^{P_{\mathrm{I}}-1}}{\overline{u}}_{\mathrm{pI}}(t^{\prime },p)H_{\mathrm{I}}({\overline{f}}_{\mathrm{I}}+k_{\mathrm{I}}{t}^{\prime },T_{\mathrm{pI}}{p}_{\mathrm{I}}){\displaystyle {\int }_{\varsigma }}\exp (\mathrm{j}(\pi \mathrm{\Delta }\tilde{k}{\varsigma }^2))$
  		$\times g(\varsigma )\exp (-\mathrm{j}(2\pi (\mathrm{\Delta }{\tilde{f}}_0\varsigma +\mathrm{\Delta }\tilde{k}{t}^{\prime }\varsigma )))\mathrm{d}\varsigma$		(41)

  where we have further defined ${\overline{f}}_{\mathrm{I}}=f_{\mathrm{I}}+k_{\mathrm{I}}{\overline{T}}_0$ and $\mathrm{\Delta }{\tilde{f}}_0=\mathrm{\Delta }{f}_0+k_{\mathrm{I}}{\overline{T}}_0$.

  The non-coherently demixed signal with which the AAF impulse response is convolved is a chirp. Thus, similarly to the derivation for the multipath channel, we may describe the convolution result as approximately a multiplication with the AAF frequency transfer function evaluated over the frequency course of this chirp. However, the term $\exp (\mathrm{j}(\pi \mathrm{\Delta }\tilde{k}{\varsigma }^2))$ cannot in general be neglected as the length of the AAF impulse response is long enough so that $\mathrm{\Delta }\tilde{k}{\varsigma }^2\not\approx 0$. Hence, we define the modified AAF impulse response as $\overline{g}(t^{\prime })=g(t^{\prime })\exp (\mathrm{j}(\pi \mathrm{\Delta }\tilde{k}t_{}^{\prime }{}_{}{}^2))$ and its frequency transfer function $\overline{G}(f)$. In this way, the remaining expression is of an analogous form to (32), finally yielding (II-A).

  A-C Inference Model

  While our derivations yield parametric models for all terms of (1), the full model contains an unknown number of interferers. These cause interference bursts distributed over the received ramps depending on the victim and interferer transmit parameters as well as the propagation channel. To relax the resulting inference problem, we make the following further model assumptions in this work:

  Assumption 3

  We limit the maximum amount of interference bursts at a certain slow-time index $p$ to one, i.e., the sum over $p_{\mathrm{I}}$ in (II-A) vanishes.
  If the number of bursts is known or estimated a-priori, our algorithm can be trivially adapted, although signal separation becomes increasingly ill-posed with a growing number of bursts. If it is unknown, it may be included in the inference model by considering a structured dictionary [78]. Algorithms for problems of this kind within the variational framework have been proposed in recent works [63, 79].

  Assumption 4

  We consider every interference burst over the slow-time independently. The channel term in (38) therefore becomes for every ramp $p$ a different delay-only line spectrum with frequency transfer function $H_{\mathrm{I}}^{(p)}(n)$.
  For a single interferer, it is possible to consider the interference signal as the sequence of all received bursts propagating through a delay-Doppler channel. We have not investigated such an approach in this work, and note that more complex model assumptions can increase the risk of major model errors, as well as lead to increased algorithm complexity.

  The above assumptions together mean that the sum ${\sum }_{i=0}^{M_{\mathrm{I}}-1}{r}_{\mathrm{I},i}(n,p)$ can be replaced by a single equivalent $r_{\mathrm{I}}(n,p)$ irrespective of $M_{\mathrm{I}}$. Furthermore, terms of this interference signal $r_{\mathrm{I}}(n,p)$ that do not explicitly depend on the fast-time index $n$ become constants for inference. These terms include the second complex phasor factor of (39), as well as certain terms of the transfer functions $H_{\mathrm{I}}(f)$ and $\overline{G}(f)$. These constants are simply inferred within the complex channel weights ${𝜷}^{(p)}$. Hence, the model is finally rewritten as (8).

  Appendix B Mean-field Variational EM Updates

  For the structured mean-field proxy PDF $q={\prod }_{𝒬}{q}_i$, the ELBO (19) can be recast to express its dependence of any one factor $q_i$ as [64, 61]

  	$ℒ(q)$	$={⟨\log p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)⟩}_q-{⟨\log q⟩}_q$
  		$={⟨{⟨\log p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)⟩}_{{\overline{q}}_i}⟩}_{q_i}$
  		$-{⟨\log q_i⟩}_{q_i}-E_i$
  		$=Z_i+{⟨\log {\overline{p}}_i⟩}_{q_i}-{⟨\log q_i⟩}_{q_i}-E_i$		(42)

  with

  	$E_i$	$={\displaystyle \sum _{𝒬\setminus q_i}}{⟨\log q_j⟩}_{q_j}$		(43)
  	${\overline{p}}_i$	$=Z_i^{}{}_{}{}^{-1}\exp \left({⟨\log p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)⟩}_{{\overline{q}}_i}\right)$		(44)
  	$Z_i$	$={\displaystyle \int \exp \left({⟨\log p(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)⟩}_{{\overline{q}}_i}\right)d{𝒛}_i}$		(45)

  where ${𝒛}_i$ denotes the latent variables corresponding to the factor $q_i$. With all other latent variables ${\overline{q}}_i$ as well as parameters fixed, the term $E_i$ and the normalization factor $Z_i$ are both independent of $q_i$. The remaining terms form a negative Kullback-Leibler divergence (KLD), which is maximized for $q_i\propto {\overline{p}}_i$, yielding (21). When only ${\overline{q}}_i$ are fixed (i.e., parameters are to be optimized jointly with $q_i$), $E_i$ is still a constant, but $Z_i$ is a function of the parameters. The properties of the negative KLD term remain, still yielding the solution $q_i\propto {\overline{p}}_i$ for the factor proxy PDF. Hence, the ELBO is maximized when $Z_i$ is maximum, which directly leads to (22) for parameter estimation.

  We proceed to summarize the solution steps needed to obtain the update equations for our model. Writing out (14), the joint PDF reads

  	$\log p$	$(𝒓,𝜶,{𝜸}_{\alpha },𝜷,{𝜸}_{\beta },\lambda ;𝜻,𝜽)$
  		$=-PN\log \pi +PN\log \lambda$
  		$-\lambda {(𝒓-𝚽𝜶-𝑼𝚿𝜷)}^{\mathrm{H}}(𝒓-𝚽𝜶-𝑼𝚿𝜷)$
  		$+(e_0-1)\log \lambda -f_0\lambda +e_0\log f_0-\log (\mathrm{\Gamma }(e_0))$
  		$+{\displaystyle \sum _l}(-\log \pi +\log {\gamma }_{\alpha ,l}-|{\alpha }_l{|}^2{\gamma }_{\alpha ,l}$
  		$+(a_0-1)\log {\gamma }_{\alpha ,l}-b_0{\gamma }_{\alpha ,l}$
  		$+a_0\log b_0-\log (\mathrm{\Gamma }(a_0)))$
  		$+{\displaystyle \sum _p}{\displaystyle \sum _k}(-\log \pi +\log {\gamma }_{\beta ,k}^{(p)}-|{\beta }_k^{(p)}{|}^2{\gamma }_{\beta ,k}^{(p)}$
  		$+(c_0-1)\log {\gamma }_{\beta ,k}^{(p)}-d_0{\gamma }_{\beta ,k}^{(p)}$
  		$+c_0\log d_0-\log (\mathrm{\Gamma }(c_0))).$		(46)

  B-A Updates of Channel Weights

  Using (21) to solve for the object channel weights $q_{\alpha }$, first note that the expectation of (46) with respect to $q_{\lambda }$ and the $q_{\gamma }$’s is obtained simply by applying $\lambda \to \widehat{\lambda }$ and $\gamma \to \widehat{\gamma }$. Neglecting constant terms independent of $𝜶$, we are left with

  $$\begin{array}{c}\log q_{\alpha }{\propto }^{\mathrm{e}}-\widehat{\lambda }{⟨{(𝒓-\widehat{𝚽}𝜶-\widehat{𝑼}\widehat{𝚿}𝜷)}^{\mathrm{H}}(𝒓-\widehat{𝚽}𝜶-\widehat{𝑼}\widehat{𝚿}𝜷)⟩}_{q_{\beta }}\hfill \\ \hfill -{𝜶}^{\mathrm{H}}{\widehat{𝚪}}_{\alpha }𝜶\end{array}$$

  (47)

  where ${\propto }^{\mathrm{e}}$ denotes proportionality after taking the exponential, and $𝚪=\mathrm{Diag}(𝜸)$. For the remaining expectation we have, using the shorthand ${𝒓}_{\beta }=𝒓-\widehat{𝚽}𝜶$,

  	${⟨{({𝒓}_{\beta }-\widehat{𝑼}\widehat{𝚿}𝜷)}^{\mathrm{H}}({𝒓}_{\beta }-\widehat{𝑼}\widehat{𝚿}𝜷)⟩}_{q_{\beta }}$
  	$={𝒓}_{\beta }^{\mathrm{H}}{𝒓}_{\beta }-{𝒓}_{\beta }^{\mathrm{H}}\widehat{𝑼}\widehat{𝚿}\widehat{𝜷}-{\widehat{𝜷}}^{\mathrm{H}}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}{𝒓}_{\beta }+{⟨{𝜷}^{\mathrm{H}}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\widehat{𝑼}\widehat{𝚿}𝜷⟩}_{q_{\beta }}.$		(48)

  We can further rewrite this, as

  	${⟨{𝜷}^{\mathrm{H}}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\widehat{𝑼}\widehat{𝚿}𝜷⟩}_{q_{\beta }}$	$={⟨{(𝜷-\widehat{𝜷})}^{\mathrm{H}}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\widehat{𝑼}\widehat{𝚿}(\beta -\widehat{𝜷})⟩}_{q_{\beta }}$
  		$+{\widehat{𝜷}}^{\mathrm{H}}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\widehat{𝜷}$		(49)

  and

  	${⟨{(𝜷-\widehat{𝜷})}^{\mathrm{H}}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\widehat{𝑼}\widehat{𝚿}(\beta -\widehat{𝜷})⟩}_{q_{\beta }}$
  	$=\mathrm{tr}\left(\widehat{𝑼}\widehat{𝚿}{⟨(𝜷-\widehat{𝜷}){(𝜷-\widehat{𝜷})}^{\mathrm{H}}⟩}_{q_{\beta }}{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\right)$
  	$=\mathrm{tr}\left(\widehat{𝑼}\widehat{𝚿}{\widehat{𝑪}}_{\beta }{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\right).$		(50)

  The trace term of (50) is independent of $𝜶$ and can be neglected, hence defining ${\widehat{𝒓}}_{\alpha }=𝒓-\widehat{𝑼}\widehat{𝚿}\widehat{𝜷}$ and recasting the rest of the terms in the form of a complex normal log-PDF yields

  	$\log q_{\alpha }$	${\propto }^{\mathrm{e}}-\widehat{\lambda }{({\widehat{𝒓}}_{\alpha }-\widehat{𝚽\alpha })}^{\mathrm{H}}({\widehat{𝒓}}_{\alpha }-\widehat{𝚽\alpha })-{𝜶}^{\mathrm{H}}{\widehat{𝚪}}_{\alpha }𝜶$
  		$=-{𝜶}^{\mathrm{H}}\left(\widehat{\lambda }{\widehat{𝚽}}^{\mathrm{H}}\widehat{𝚽}+{\widehat{𝚪}}_{\alpha }\right)𝜶+2\mathrm{\Re }({\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}𝜶)$
  		$-{\widehat{\lambda }}^2{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}{\left(\widehat{\lambda }{\widehat{𝚽}}^{\mathrm{H}}\widehat{𝚽}+{\widehat{𝚪}}_{\alpha }\right)}^{-1}{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }$
  		$+{\widehat{\lambda }}^2{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}{\left(\widehat{\lambda }{\widehat{𝚽}}^{\mathrm{H}}\widehat{𝚽}+{\widehat{𝚪}}_{\alpha }\right)}^{-1}{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }$
  		$=\log \text{CN}(𝜶|\widehat{𝜶},{\widehat{𝑪}}_{\alpha })+\log |{\widehat{𝑪}}_{\alpha }|$
  		$+{\widehat{\lambda }}^2{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }$		(51)

  where by coefficient comparison we obtain the result

  	$q(𝜶)$	$=\text{CN}(𝜶|\widehat{𝜶},{\widehat{𝑪}}_{\alpha })$		(52)
  	$\widehat{𝜶}$	$=\widehat{\lambda }{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }$		(53)
  	${\widehat{𝑪}}_{\alpha }$	$={\left(\widehat{\lambda }{\widehat{𝚽}}^{\mathrm{H}}\widehat{𝚽}+{\widehat{𝚪}}_{\alpha }\right)}^{-1}.$		(54)

  The derivations for the interference channel weights $𝜷$ are analogous, leading to

  	$q(𝜷)$	$=\text{CN}(𝜷|\widehat{𝜷},{\widehat{𝑪}}_{\beta })$		(55)
  	$\widehat{𝜷}$	$=\widehat{\lambda }{\widehat{𝑪}}_{\beta }{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}{\widehat{𝒓}}_{\beta }$		(56)
  	${\widehat{𝑪}}_{\beta }$	$={\left(\widehat{\lambda }{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\widehat{𝑼}\widehat{𝚿}+{\widehat{𝚪}}_{\beta }\right)}^{-1}.$		(57)

  B-B Update of Noise Precision

  Similarly for the noise precision update, noting only $\lambda$-dependent terms we have

  $$\begin{array}{cc}\hfill \log q_{\lambda }& {\propto }^{\mathrm{e}}(e_0+PN-1)\log \lambda -f_0\lambda \hfill \\ & -\lambda {⟨{⟨{(𝒓-\widehat{𝚽}𝜶-\widehat{𝑼}\widehat{𝚿}𝜷)}^{\mathrm{H}}(𝒓-\widehat{𝚽}𝜶-\widehat{𝑼}\widehat{𝚿}𝜷)⟩}_{q_{\beta }}⟩}_{q_{\alpha }}.\hfill \end{array}$$

  The expectations are analogous to the one in (47), denoting ${\widehat{𝒓}}_{\lambda }=𝒓-\widehat{𝚽}\widehat{𝜶}-\widehat{𝑼}\widehat{𝚿}\widehat{𝜷}$, leading to

  $$\begin{array}{cc}& \log q_{\lambda }\hfill \\ & {\propto }^{\mathrm{e}}(e_0+PN-1)\log \lambda -f_0\lambda \hfill \\ & -\lambda \left({\widehat{𝒓}}_{\lambda }^{\mathrm{H}}{\widehat{𝒓}}_{\lambda }+\mathrm{tr}\left(\widehat{𝚽}{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}}\right)+\mathrm{tr}\left(\widehat{𝑼}\widehat{𝚿}{\widehat{𝑪}}_{\beta }{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}}\right)\right)\hfill \end{array}$$

  (58)

  and therefore

  	$q(\lambda )$	$=\text{Ga}(\lambda |e,f)$		(59)
  	$\widehat{\lambda }={\displaystyle \frac{e}{f}}$	$={\displaystyle \frac{e_0+PN}{f_0+{\widehat{𝒓}}_{\lambda }^{\mathrm{H}}{\widehat{𝒓}}_{\lambda }+\mathrm{tr}(\widehat{𝚽}{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}})+\mathrm{tr}(\widehat{𝑼}\widehat{𝚿}{\widehat{𝑪}}_{\beta }{\widehat{𝚿}}^{\mathrm{H}}{\widehat{𝑼}}^{\mathrm{H}})}}.$		(60)

  For the hyper-parameters of the Gamma distribution we let $e_0=f_0=0$ that is the non-informative Jeffery’s improper hyper-prior.

  B-C Estimations of Parameters

  The cost function for the object channel dispersion parameters $𝜻$ is $Z_{\alpha }$ as in (45). For optimization, the equivalent $\log Z_{\alpha }$ will be used. The expectation term is identical to the solution for $q_{\alpha }$, except that $𝚽\mathbf{(}𝜻\mathbf{)}$-dependent terms are to be explicitly considered. The result in (51) contains all these non-constant terms. Then, marginalization of $𝜶$ integrates out $\text{CN}(𝜶|\widehat{𝜶},{\widehat{𝑪}}_{\alpha })$ to one as it is a valid PDF. I.e.,

  	$\log Z_{\alpha }$	${\propto }^{\mathrm{e}}\log {\displaystyle \int }\exp (\log \text{CN}(𝜶|\widehat{𝜶},{\widehat{𝑪}}_{\alpha })+\log |{\widehat{𝑪}}_{\alpha }|$
  		$+{\widehat{\lambda }}^2{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha })\mathrm{d}𝜶$
  		$=\log |{\widehat{𝑪}}_{\alpha }|+{\widehat{\lambda }}^2{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }$
  		$+\log {\displaystyle \int \text{CN}(𝜶|\widehat{𝜶},{\widehat{𝑪}}_{\alpha })d𝜶}$		(61)

  and hence the result

  $$\widehat{𝜻}=\underset{𝜻}{\arg \max }\left(\log |{\widehat{𝑪}}_{\alpha }(𝜻)|+{\widehat{\lambda }}^2{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}𝚽(𝜻){\widehat{𝑪}}_{\alpha }(𝜻){𝚽}^{\mathrm{H}}(𝜻){\widehat{𝒓}}_{\alpha }\right)$$

  (62)

  where we explicitly noted all $𝜻$-dependencies including inside the covariance matrix estimate ${\widehat{𝑪}}_{\alpha }(𝜻)$. This result can be recast to express the cost function for the parameters of a single component ${𝜻}_l$ with all other parameters fixed. We order the components so that the index $l$ is last and denote all others by $\overline{l}$, i.e., $𝚽(𝜻)=𝚽({[{𝜻}_{\overline{l}}{𝜻}_l]}^{\mathrm{T}})=[{𝚽}_{\overline{l}}{\mathit{\varphi }}_l]$. Then, ${\widehat{𝑪}}_{\alpha }$ can be written as the block matrix

  	${\widehat{𝑪}}_{\alpha }$	$={\left[\begin{array}{cc}{\widehat{𝑪}}_{\alpha ,\overline{l}}^{-1}& \widehat{\lambda }{𝚽}_{\overline{l}}^{\mathrm{H}}{\mathit{\varphi }}_l\\ \widehat{\lambda }{\mathit{\varphi }}_l^{\mathrm{H}}{𝚽}_{\overline{l}}& {\widehat{c}}_{\alpha ,l}^{-1}\end{array}\right]}^{-1}$		(63)
  		$=\left[\begin{array}{cc}{\widehat{𝑪}}_{\alpha ,\overline{l}}+& -\widehat{\lambda }{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}{\mathit{\varphi }}_l{\widehat{c}}_{\alpha ,l}^{\prime }\\ {\widehat{\lambda }}^2{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}{\mathit{\varphi }}_l{\widehat{c}}_{\alpha ,l}^{\prime }{\mathit{\varphi }}_l^{\mathrm{H}}{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}& \\ -\widehat{\lambda }{\widehat{c}}_{\alpha ,l}^{\prime }{\mathit{\varphi }}_l^{\mathrm{H}}{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}& {\widehat{c}}_{\alpha ,l}^{\prime }\end{array}\right]$

  where we have applied the formula for block matrix inversion and defined ${\widehat{𝑪}}_{\alpha ,\overline{l}}={(\widehat{\lambda }{𝚽}_{\overline{l}}^{\mathrm{H}}{𝚽}_{\overline{l}}+{\widehat{𝚪}}_{\alpha ,\overline{l}})}^{-1}$, ${\widehat{c}}_{\alpha ,l}={(\widehat{\lambda }{\mathit{\varphi }}_l^{\mathrm{H}}{\mathit{\varphi }}_l+{\widehat{\gamma }}_{\alpha ,l})}^{-1}$ and ${\widehat{c}}_{\alpha ,l}^{\prime }={({\widehat{c}}_{\alpha ,l}^{-1}-{\widehat{\lambda }}^2{\mathit{\varphi }}_l^{\mathrm{H}}{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}{\mathit{\varphi }}_l)}^{-1}$. In (62), for the log-determinant of the block matrix we simply obtain $\log |{\widehat{𝑪}}_{\alpha }|=\log |{\widehat{c}}_{\alpha ,l}^{\prime }|+\log |{\widehat{𝑪}}_{\alpha ,\overline{l}}|$. ${\widehat{𝑪}}_{\alpha ,\overline{l}}$ is independent of ${𝜻}_l$ and hence can be neglected in the resulting cost function. Concerning the other term of (62), noting that

  	$𝚽{\widehat{𝑪}}_{\alpha }{𝚽}^{\mathrm{H}}$	$={𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}+{\widehat{\lambda }}^2{\widehat{c}}_{\alpha ,l}^{\prime }{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}{\mathit{\varphi }}_l{\mathit{\varphi }}_l^{\mathrm{H}}{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}$
  		$-2\mathrm{\Re }\left(\widehat{\lambda }{\widehat{c}}_{\alpha ,l}^{\prime }{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}{\mathit{\varphi }}_l{\mathit{\varphi }}_l^{\mathrm{H}}\right)+{\mathit{\varphi }}_l{\widehat{c}}_{\alpha ,l}^{\prime }{\mathit{\varphi }}_l^{\mathrm{H}}$
  		$={𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}$		(64)
  		$+{(\widehat{\lambda }{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}-𝐈)}^{\mathrm{H}}{\mathit{\varphi }}_l{\widehat{c}}_{\alpha ,l}^{\prime }{\mathit{\varphi }}_l^{\mathrm{H}}(\widehat{\lambda }{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}-𝐈)$

  with the first term being a constant, we finally arrive at the result

  $$\begin{array}{c}{\widehat{𝜻}}_l=\underset{{𝜻}_l}{\arg \max }(\log |{\widehat{c}}_{\alpha ,l}^{\prime }({𝜻}_l)|\hfill \\ \hfill +{\widehat{\lambda }}^2{\widehat{c}}_{\alpha ,l}^{\prime }({𝜻}_l)\parallel {\mathit{\varphi }}_l^{\mathrm{H}}({𝜻}_l)(\widehat{\lambda }{𝚽}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{𝚽}_{\overline{l}}^{\mathrm{H}}-𝐈){\widehat{𝒓}}_{\alpha }{\parallel }^2).\end{array}$$

  (65)

  For the interference chirp parameter estimates the derivations are analogous to (62) for every ramp independently, i.e.,

  $$\begin{array}{c}\widehat{𝜽}=\underset{𝜽}{\arg \max }(\log |{\widehat{𝑪}}_{\beta }(𝜽)|\hfill \\ \hfill +{\widehat{\lambda }}^2{\widehat{𝒓}}_{\beta }^{\mathrm{H}}𝑼(𝜽)\widehat{𝚿}{\widehat{𝑪}}_{\beta }(𝜽){\widehat{𝚿}}^{\mathrm{H}}𝑼{(𝜽)}^{\mathrm{H}}{\widehat{𝒓}}_{\beta })\end{array}$$

  (66)

  where the dependence of variables on the ramp index $p$ has been left out to reduce notational clutter.

  B-D Fast Updates of Component Precisions

  As discussed in Section III-C, for our proposed algorithm we also apply the “fast update rule” for components as derived in [58]. Applying (21) yields

  	$\log q_{{\gamma }_{\alpha ,l}}$	${\propto }^{\mathrm{e}}{a}_0\log {\gamma }_{\alpha ,l}-(b_0+{⟨{|{\alpha }_l|}^2⟩}_{q_{\alpha }}){\gamma }_{\alpha ,l}$		(67)
  		${\propto }^{\mathrm{e}}{a}_0\log {\gamma }_{\alpha ,l}-(b_0+{|{\widehat{\alpha }}_l|}^2+{\widehat{𝑪}}_{\alpha }[l,l]){\gamma }_{\alpha ,l}$		(68)

  which leads to the solution

  	$q({\gamma }_{\alpha ,l})$	$=\text{Ga}({\gamma }_{\alpha ,l}|a_l,b_l)$		(69)
  	${\widehat{\gamma }}_{\alpha ,l}={\displaystyle \frac{a_l}{b_l}}$	$={\displaystyle \frac{a_0+1}{b_0+\left({|{\widehat{\alpha }}_l|}^2+{\widehat{𝑪}}_{\alpha }[l,l]\right)}}$		(70)

  where we take $a_0=b_0=0$. However, we can further investigate the behavior of this implicit equation when repeatedly updating $q_{\alpha }$ and $q_{{\gamma }_{\alpha ,l}}$. We insert (53) into (70) and recast the resulting terms to “isolate” the $l$-th component as already shown for the derivation of (65). I.e.,

  	${\widehat{\gamma }}_{\alpha ,l}^{-1}$	$=({\widehat{\lambda }}^2{\widehat{𝑪}}_{\alpha }{\widehat{𝚽}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }{\widehat{𝒓}}_{\alpha }^{\mathrm{H}}\widehat{𝚽}{\widehat{𝑪}}_{\alpha }^{\mathrm{H}}+{\widehat{𝑪}}_{\alpha })[l,l]$
  		$={\widehat{\lambda }}^2{|c_{\alpha ,l}^{\prime }{\widehat{\mathit{\varphi }}}_l^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }-\lambda c_{\alpha ,l}^{\prime }{\widehat{\mathit{\varphi }}}_l^{\mathrm{H}}{\widehat{𝚽}}_{\overline{l}}{\widehat{𝑪}}_{\alpha ,\overline{l}}{\widehat{𝚽}}_{\overline{l}}^{\mathrm{H}}{\widehat{𝒓}}_{\alpha }|}^2+c_{\alpha ,l}^{\prime }$		(71)

  and by defining (23) and (24), and hence $c_{\alpha ,l}^{\prime }={\widehat{\gamma }}_{\alpha ,l}^{-1}{\rho }_l/({\widehat{\gamma }}_{\alpha ,l}^{-1}+{\rho }_l)$, this can be further rewritten as

  	${\widehat{\gamma }}_{\alpha ,l}^{-1}$	$={|c_{\alpha ,l}^{\prime }|}^2{\rho }_l^{-2}{\omega }_l^2+c_{\alpha ,l}^{\prime }$
  		$={\displaystyle \frac{{\widehat{\gamma }}_{\alpha ,l}^{-2}({\omega }_l^2+{\rho }_l)+{\widehat{\gamma }}_{\alpha ,l}^{-1}{\rho }_l^2}{{|{\widehat{\gamma }}_{\alpha ,l}^{-1}+{\rho }_l|}^2}}.$		(72)

  By analyzing (72) in the context of iterative updates as a rational map, it can be found that its stationary points are located at $\{0,{\omega }_l^2-{\rho }_l\}$. It can furthermore be proven that the non-zero stationary point is reached (i.e., the estimate ${\widehat{\gamma }}_{\alpha ,l}$ converges) if and only if ${\rho }_l/{\omega }_l^2>T$ with $T=1$, which leads to the thresholding scheme with threshold $T$ (see Algorithm 4). For a detailed discussion see [80, 58]. The component precisions of the interference channel ${\widehat{\gamma }}_{\beta ,k}$ can be found in a similar manner.

  Appendix C Measurement of Interference

  As introduced in Section IV-D, a measurement scheme whereby a separate interference signal is made available was devised. The parameters of the interference measurement can be found as Measurement II in Table I. The measurement took place in an indoor office with the victim and interfering sensors on a desk opposite each other. Crucially, the scenario was therefore completely static with no movement of objects. First, the location and duration of the interference bursts received by the victim radar were identified manually. The static object signal could then be estimated by taking the sample mean of the non-interfered data. This estimate was subtracted from the measured signal, approximately removing the object signal, leaving interference and measurement noise. Then, the time-limited interference bursts were cut out to exclude irrelevant noise-only sections.

  In Section II, the effect of the AAF in the sensor hardware on the resulting interference burst was discussed. This is typically assumed to be known to some precision in the application by the algorithm designers of an automotive radar manufacturer. Section IV-B investigated the resilience of the proposed algorithm to a minor inaccuracy of this kind in the inference model. However, because the used hardware was a radar evaluation board from an external supplier, no details on the transfer function of the system was accessible to the authors. This can lead to a major model error. To mitigate this, the signal was additionally filtered digitally with a filter whose magnitude response dominates that of the sensor hardware. The concatenation of a high- and a low-pass was applied, both being 8-th order Butterworth filters with respective cut-off frequencies of $0.5\mathrm{MHz}$ and $5\mathrm{MHz}$. The known response of this digital filter chain was then used in our inference model, noting that this still constitutes a rather significant model error. Finally, the interference signal prepared in this way was added onto the measurement from Measurement I, yielding our test radar signal.

  References

  • [1]
• M. Kunert, H. Meinel, C. Fischer, and M. Ahrholdt, “D16.1 - Report on interference density increase by market penetration forecast,” the MOSARIM Consortium, Tech. Rep., 2010.
• [2]
• F. Roos, J. Bechter, C. Knill, B. Schweizer, and C. Waldschmidt, “Radar Sensors for Autonomous Driving: Modulation Schemes and Interference Mitigation,” IEEE Microwave Magazine, vol. 20, no. 9, pp. 58–72, Sep. 2019.
• [3]
• M. Kunert, “The EU project MOSARIM: A general overview of project objectives and conducted work,” in 2012 9th European Radar Conference, Oct. 2012, pp. 1–5. [Online]. Available: https://ieeexplore.ieee.org/document/6450766
• [4]
• the IMIKO Consortium, “Gesamtbewertung und Empfehlung für Standardisierung,” Deliverable AP09, May 2022, in German. [Online]. Available: https://www.elektronikforschung.de/dateien/projekte/imiko-radar-ap09-gesamtbewertungen-der-rim-verfahren.pdf
• [5]
• M. Hischke, “Collision warning radar interference,” in IEEE Intell. Veh. Symp., 1995, pp. 13–18.
• [6]
• B. E. Tullsson, “Topics in FMCW radar disturbance suppression,” in Radar 97 (Conf. Publ. No. 449), 1997, pp. 1–5.
• [7]
• G. M. Brooker, “Mutual interference of millimeter-wave radar systems,” IEEE Trans. Electromagn. Compat., vol. 49, no. 1, pp. 170–181, 2007.
• [8]
• A. Al-Hourani, R. J. Evans, S. Kandeepan, B. Moran, and H. Eltom, “Stochastic Geometry Methods for Modeling Automotive Radar Interference,” IEEE Trans. Intell. Transp. Syst., pp. 1–11, 2017.
• [9]
• L. L. T. Torres and C. Waldschmidt, “Analysis of Automotive Radar Interference in Complex Traffic Scenarios using Graph Theory,” in 2022 23rd International Radar Symposium (IRS), Sep. 2022, pp. 269–274, iSSN: 2155-5753.
• [10]
• T. Schipper, S. Prophet, M. Harter, L. Zwirello, and T. Zwick, “Simulative Prediction of the Interference Potential Between Radars in Common Road Scenarios,” IEEE Trans. Electromagn. Compat., vol. 57, no. 3, pp. 322–328, 2015.
• [11]
• M. Toth, P. Meissner, A. Melzer, and K. Witrisal, “Analytical Investigation of Non-Coherent Mutual FMCW Radar Interference,” 2018.
• [12]
• G. Kim, J. Mun, and J. Lee, “A Peer-to-peer Interference Analysis for Automotive Chirp Sequence Radars,” IEEE Trans. Veh. Technol., vol. PP, no. 1, pp. 1–1, 2018. [Online]. Available: https://ieeexplore.ieee.org/document/8388252/
• [13]
• A. Bourdoux and M. Bauduin, “PMCW Waveform Cross-correlation Characterization and Interference Mitigation,” in 2020 17th European Radar Conference (EuRAD), Jan. 2021, pp. 164–167.
• [14]
• B. Schweizer, C. Knill, D. Werbunat, S. Stephany, and C. Waldschmidt, “Mutual Interference of Automotive OFDM Radars—Analysis and Countermeasures,” IEEE Journal of Microwaves, vol. 1, no. 4, pp. 950–961, 2021. [Online]. Available: https://ieeexplore.ieee.org/document/9552961/
• [15]
• M. Rameez, M. Dahl, and M. I. Pettersson, “Experimental Evaluation of Adaptive Beamforming for Automotive Radar Interference Suppression,” in 2020 IEEE Radio and Wireless Symposium (RWS), Jan. 2020, pp. 183–186, iSSN: 2164-2974.
• [16]
• A. Ossowska, L. Sit, S. Manchala, T. Vogler, K. Krupinski, and U. Luebbert, “IMIKO-Radar Project: Laboratory Interference Measurements of Automotive Radar Sensors,” in 2020 21st International Radar Symposium (IRS), Oct. 2020, pp. 334–338, iSSN: 2155-5753.
• [17]
• L. L. T. Torres, M. Steiner, and C. Waldschmidt, “Channel Influence for the Analysis of Interferences Between Automotive Radars,” 2020, p. 4.
• [18]
• N. Cardona, J. S. Romero, W. Yang, and J. Li, “Integrating the Sensing and Radio Communications Channel Modelling From Radar Mutual Interference,” in ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Jun. 2023, pp. 1–5.
• [19]
• G. Hakobyan, K. Armanious, and B. Yang, “Interference-Aware Cognitive Radar: A Remedy to the Automotive Interference Problem,” IEEE Transactions on Aerospace and Electronic Systems, pp. 1–1, 2019.
• [20]
• L. L. Tovar Torres, T. Grebner, and C. Waldschmidt, “Automotive Radar Interference Avoidance Strategies for Complex Traffic Scenarios,” in 2023 IEEE Radar Conference (RadarConf23), May 2023, pp. 1–6.
• [21]
• J. Khoury, R. Ramanathan, D. McCloskey, R. Smith, and T. Campbell, “RadarMAC: Mitigating radar interference in self-driving cars,” in 2016 13th Annu. IEEE Int. Conf. Sensing, Commun. Networking, SECON 2016.   IEEE, 2016, pp. 1–9.
• [22]
• C. Aydogdu, M. F. Keskin, N. Garcia, H. Wymeersch, and D. W. Bliss, “Radchat: Spectrum sharing for automotive radar interference mitigation,” IEEE Transactions on Intelligent Transportation Systems, vol. 22, no. 1, pp. 416–429, 2021.
• [23]
• C. Fischer, H. L. Blöcher, J. Dickmann, and W. Menzel, “Robust detection and mitigation of mutual interference in automotive radar,” in 2015 16th International Radar Symposium (IRS), Jun. 2015, pp. 143–148.
• [24]
• Z. Liu, W. Lu, J. Wu, S. Yang, and G. Li, “A PELT-KCN Algorithm for FMCW Radar Interference Suppression Based on Signal Reconstruction,” IEEE Access, pp. 1–1, 2020.
• [25]
• T. Pernstål, J. Degerman, H. Broström, V. T. Vu, and M. I. Pettersson, “GIP test for Automotive FMCW interference Detection and Suppression,” in 2020 IEEE Radar Conference (RadarConf20), Sep. 2020, pp. 1–6, iSSN: 2375-5318.
• [26]
• T. Shimura, M. Umehira, Y. Watanabe, X. Wang, and S. Takeda, “An Advanced Wideband Interference Suppression Technique using Envelope Detection and Sorting for Automotive FMCW Radar,” in 2022 IEEE Radar Conference (RadarConf22), Mar. 2022, pp. 1–6.
• [27]
• M. Rameez, M. Dahl, and M. I. Pettersson, “Autoregressive Model-Based Signal Reconstruction for Automotive Radar Interference Mitigation,” IEEE Sensors Journal, pp. 1–1, 2020.
• [28]
• S. Neemat, O. Krasnov, and A. Yarovoy, “An Interference Mitigation Technique for FMCW Radar Using Beat-Frequencies Interpolation in the STFT Domain,” IEEE Transactions on Microwave Theory and Techniques, vol. 67, no. 3, pp. 1207–1220, Mar. 2019. [Online]. Available: https://ieeexplore.ieee.org/document/8551288/
• [29]
• J. Bechter, F. Roos, M. Rahman, and C. Waldschmidt, “Automotive Radar Interference Mitigation Using a Sparse Sampling Approach,” in Proc. 14th Eur. Radar Conf., 2017, pp. 90–93.
• [30]
• S. Chen, J. Taghia, U. Kühnau, T. Fei, F. Grünhaupt, and R. Martin, “Automotive Radar Interference Reduction Based on Sparse Bayesian Learning,” in 2020 IEEE Radar Conference (RadarConf20), Sep. 2020, pp. 1–6, iSSN: 2375-5318.
• [31]
• J. Wang, M. Ding, and A. Yarovoy, “Matrix-Pencil Approach-Based Interference Mitigation for FMCW Radar Systems,” Feb. 2021.
• [32]
• J. Bechter, A. Demirlika, P. Hügler, F. Roos, and C. Waldschmidt, “Blind Adaptive Beamforming for Automotive Radar Interference Suppression,” in 2018 19th International Radar Symposium (IRS), Jun. 2018, pp. 1–10, iSSN: 2155-5753.
• [33]
• M. Rameez, M. Dahl, and M. I. Pettersson, “Adaptive Digital Beamforming for Interference Suppression in Automotive FMCW Radars,” in 2018 IEEE Radar Conf., 2018, pp. 252–256.
• [34]
• F. Jin and S. Cao, “Automotive Radar Interference Mitigation Using Adaptive Noise Canceller,” IEEE Transactions on Vehicular Technology, vol. 68, no. 4, pp. 3747–3754, Apr. 2019. [Online]. Available: https://ieeexplore.ieee.org/document/8651538/
• [35]
• P. Wang, X. Yin, J. Rodríguez-Piñeiro, Z. Chen, P. Zhu, and G. Li, “A Dual-Recursive-Least-Squares Algorithm for Automotive Radar Interference Suppression,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–15, 2023, conference Name: IEEE Transactions on Intelligent Transportation Systems.
• [36]
• M. Rameez, M. I. Pettersson, and M. Dahl, “Interference Compression and Mitigation for Automotive FMCW Radar Systems,” IEEE Sensors Journal, vol. 22, no. 20, pp. 19 739–19 749, Oct. 2022, conference Name: IEEE Sensors Journal.
• [37]
• M. Toth, P. Meissner, A. Melzer, and K. Witrisal, “Slow-Time Mitigation of Mutual Interference in Chirp Sequence Radar,” in 2020 IEEE MTT-S International Conference on Microwaves for Intelligent Mobility (ICMIM), Nov. 2020, pp. 1–4.
• [38]
• M. Wagner, F. Sulejmani, A. Melzer, P. Meissner, and M. Huemer, “Threshold-Free Interference Cancellation Method for Automotive FMCW Radar Systems,” 2018.
• [39]
• R. Muja, A. Anghel, R. Cacoveanu, and S. Ciochina, “Interference Mitigation in FMCW Automotive Radars using the Short-Time Fourier Transform and L-Statistics,” in 2022 IEEE Radar Conference (RadarConf22), Mar. 2022, pp. 1–6.
• [40]
• J. Wu, S. Yang, W. Lu, and Z. Liu, “Iterative modified threshold method based on EMD for interference suppression in FMCW radars,” Sonar Navigation IET Radar, vol. 14, no. 8, pp. 1219–1228, 2020.
• [41]
• T. Oyedare, V. K. Shah, D. J. Jakubisin, and J. H. Reed, “Interference Suppression Using Deep Learning: Current Approaches and Open Challenges,” IEEE Access, vol. 10, pp. 66 238–66 266, 2022, conference Name: IEEE Access.
• [42]
• J. Hille, D. Auge, C. Grassmann, and A. Knoll, “FMCW radar2radar Interference Detection with a Recurrent Neural Network,” in 2022 IEEE Radar Conference (RadarConf22), Mar. 2022, pp. 1–6.
• [43]
• J. Rock, W. Roth, M. Toth, P. Meissner, and F. Pernkopf, “Resource-Efficient Deep Neural Networks for Automotive Radar Interference Mitigation,” IEEE Journal of Selected Topics in Signal Processing, vol. 15, no. 4, pp. 927–940, Jun. 2021, conference Name: IEEE Journal of Selected Topics in Signal Processing.
• [44]
• J. Fuchs, A. Dubey, M. Lübke, R. Weigel, and F. Lurz, “Automotive Radar Interference Mitigation using a Convolutional Autoencoder,” in 2020 IEEE International Radar Conference (RADAR), Apr. 2020, pp. 315–320, iSSN: 2640-7736.
• [45]
• C. Oswald, M. Toth, P. Meissner, and F. Pernkopf, “End-to-End Training of Neural Networks for Automotive Radar Interference Mitigation,” in 2023 IEEE International Radar Conference (RADAR).   Sydney, Australia: IEEE, Nov. 2023, pp. 1–6. [Online]. Available: https://ieeexplore.ieee.org/document/10371003/
• [46]
• J. Bechter, K. D. Biswas, and C. Waldschmidt, “Estimation and cancellation of interferences in automotive radar signals,” in Proc. Int. Radar Symp., 2017, pp. 1–10.
• [47]
• A. Correas-Serrano and M. A. Gonzalez-Huici, “Sparse Reconstruction of Chirplets for Automotive FMCW Radar Interference Mitigation,” in 2019 IEEE MTT-S International Conference on Microwaves for Intelligent Mobility (ICMIM).   Detroit, MI, USA: IEEE, Apr. 2019, pp. 1–4. [Online]. Available: https://ieeexplore.ieee.org/document/8726758/
• [48]
• L. L. T. Torres, T. Grebner, D. Werbunat, and C. Waldschmidt, “Automotive Radar Interference Mitigation by Subtraction of the Interference Component,” IEEE Microwave and Wireless Technology Letters, pp. 1–0, 2023, conference Name: IEEE Microwave and Wireless Technology Letters.
• [49]
• F. Uysal, “Synchronous and Asynchronous Radar Interference Mitigation,” IEEE Access, vol. PP, no. 1, p. 1, 2018.
• [50]
• J.-L. Starck, M. Elad, and D. Donoho, “Redundant Multiscale Transforms and Their Application for Morphological Component Separation,” Advances in Imaging and Electron Physics, vol. 132, pp. 287–348, Jan. 2004. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1076567004320069
• [51]
• S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Rev., vol. 43, no. 1, pp. 129–159, 2001. [Online]. Available: https://doi.org/10.1137/S003614450037906X
• [52]
• Z. Xu and M. Yuan, “An Interference Mitigation Technique for Automotive Millimeter Wave Radars in the Tunable Q-Factor Wavelet Transform Domain,” IEEE Transactions on Microwave Theory and Techniques, pp. 1–1, 2021, conference Name: IEEE Transactions on Microwave Theory and Techniques.
• [53]
• Z. Xu, “Bi-Level $l_1$ Optimization-Based Interference Reduction for Millimeter Wave Radars,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–11, 2022, conference Name: IEEE Transactions on Intelligent Transportation Systems.
• [54]
• J. Wang, M. Ding, and A. Yarovoy, “Interference Mitigation for FMCW Radar With Sparse and Low-Rank Hankel Matrix Decomposition,” arXiv:2106.06748 [eess], Jun. 2021, arXiv: 2106.06748. [Online]. Available: http://arxiv.org/abs/2106.06748
• [55]
• Y. Wang, Y. Huang, C. Wen, X. Zhou, J. Liu, and W. Hong, “Mutual Interference Mitigation for Automotive FMCW Radar With Time and Frequency Domain Decomposition,” IEEE Transactions on Microwave Theory and Techniques, pp. 1–17, 2023, conference Name: IEEE Transactions on Microwave Theory and Techniques.
• [56]
• S. Grebien, E. Leitinger, K. Witrisal, and B. H. Fleury, “Super-resolution estimation of UWB channels including the dense component – An SBL-inspired approach,” IEEE Trans. Wireless Commun., Feb. 2024.
• [57]
• M. E. Tipping and A. C. Faul, “Fast marginal likelihood maximisation for sparse bayesian models,” in Proc. 9th Int. Workshop Artif. Intell. and Statist., vol. R4, Jan. 03–06, 2003, pp. 276–283. [Online]. Available: https://proceedings.mlr.press/r4/tipping03a.html
• [58]
• D. Shutin, T. Buchgraber, S. R. Kulkarni, and H. V. Poor, “Fast variational sparse bayesian learning with automatic relevance determination for superimposed signals,” IEEE Trans. Signal Process., vol. 59, no. 12, pp. 6257–6261, Dec. 2011.
• [59]
• D. Shutin, W. Wand, and T. Jost, “Incremental sparse bayesian learning for parameter estimation of superimposed signals,” in 10th Int. Conf. Sampling Theory and Appl., Bremen, Germany, Jul. 1–5, 2013, pp. 513–516. [Online]. Available: https://eurasip.org/Proceedings/Ext/SampTA2013/papers/p513-shutin.pdf
• [60]
• D. P. Wipf, B. D. Rao, and S. Nagarajan, “Latent variable bayesian models for promoting sparsity,” IEEE Trans. Image Process., vol. 57, no. 9, pp. 6236–6255, Sep. 2011.
• [61]
• D. G. Tzikas, A. C. Likas, and N. P. Galatsanos, “The variational approximation for Bayesian inference,” IEEE Signal Processing Magazine, vol. 25, no. 6, pp. 131–146, Nov. 2008. [Online]. Available: https://ieeexplore.ieee.org/document/4644060/
• [62]
• M.-A. Badiu, T. L. Hansen, and B. H. Fleury, “Variational bayesian inference of line spectra,” IEEE Trans. Signal Process., vol. 65, no. 9, pp. 2247–2261, May 2017.
• [63]
• J. Möderl, F. Pernkopf, K. Witrisal, and E. Leitinger, “Variational Inference of Structured Line Spectra Exploiting Group-Sparsity,” May 2023, arXiv:2303.03017 [eess]. [Online]. Available: http://arxiv.org/abs/2303.03017
• [64]
• C. M. Bishop, Pattern recognition and machine learning, 8th ed., M. Jordan, J. Kleinberg, and B. Schölkopf, Eds.   New York, NY, USA: Springer-Verlag, 2009.
• [65]
• T. L. Hansen, B. H. Fleury, and B. D. Rao, “Superfast Line Spectral Estimation,” IEEE Transactions on Signal Processing, vol. 66, no. 10, pp. 2511–2526, May 2018, conference Name: IEEE Transactions on Signal Processing.
• [66]
• E. Leitinger, S. Grebien, B. Fleury, and K. Witrisal, “Detection and estimation of a spectral line in MIMO systems,” in 2020 54th Asilomar Conf. Signals, Syst. and Computers, Pacific Grove, CA, USA, Nov. 01–04, 2020, pp. 1090–1095.
• [67]
• S. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory.   Prentice Hall, 1993.
• [68]
• M. A. Richards, Fundamentals of radar signal processing, 2nd ed.   New York, NY: McGraw-Hill Education, 2014, oCLC: 931498102.
• [69]
• T. Schipper, T. Mahler, M. Harter, L. Reichardt, and T. Zwick, “An Estimation of the Operating Range for Frequency Modulated Radars in the Presence of Interference,” Proc. 10th Eur. Radar Conf., pp. 227–230, 2013.
• [70]
• T. Schipper, “SIMULATION OF EFFECTS AND IMPACT OF ENVIRONMENT , TRAFFIC PARTICIPANTS AND INFRASTRUCTURE,” the MOSARIM Consortium, Tech. Rep., 2012.
• [71]
• A. S. Rahmathullah, A. F. Garcia-Fernandez, and L. Svensson, “Generalized optimal sub-pattern assignment metric,” in 2017 20th International Conference on Information Fusion (Fusion), Jul. 2017, pp. 1–8, arXiv:1601.05585 [cs]. [Online]. Available: http://arxiv.org/abs/1601.05585
• [72]
• J.-L. Starck, F. Murtagh, and J. M. Fadili, Sparse image and signal processing: wavelets, curvelets, morphological diversity.   Cambridge ; New York: Cambridge University Press, 2010, oCLC: ocn471237449.
• [73]
• M. Toth, J. Rock, P. Meissner, A. Melzer, and K. Witrisal, “Analysis of automotive radar interference mitigation for real-world environments,” in 2020 17th European Radar Conference (EuRAD), 2021, pp. 176–179.
• [74]
• J. Rock, M. Toth, P. Meissner, and F. Pernkopf, “Deep interference mitigation and denoising of real-world fmcw radar signals,” in 2020 IEEE International Radar Conference (RADAR), 2020, pp. 624–629.
• [75]
• A. Stove, “Linear FMCW radar techniques,” IEE Proc. F Radar Signal Process., vol. 139, no. 5, p. 343, 1992. [Online]. Available: http://digital-library.theiet.org/content/journals/10.1049/ip-f-2.1992.0048
• [76]
• F. Hlawatsch and G. Matz, Eds., Wireless Communications Over Rapidly Time-Varying Channels.   Elsevier, 2011. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/C20090185069
• [77]
• W. J. Caputi, “Stretch: A Time-Transformation Technique,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-7, no. 2, pp. 269–278, Mar. 1971.
• [78]
• D. Wipf, “Sparse Estimation with Structured Dictionaries,” in Advances in Neural Information Processing Systems, vol. 24.   Curran Associates, Inc., 2011. [Online]. Available: https://papers.nips.cc/paper_files/paper/2011/hash/fd2c5e4680d9a01dba3aada5ece22270-Abstract.html
• [79]
• J. Möderl, F. Pernkopf, K. Witrisal, and E. Leitinger, “Fast Variational Block-Sparse Bayesian Learning,” Dec. 2023, arXiv:2306.00442 [eess]. [Online]. Available: http://arxiv.org/abs/2306.00442
• [80]
T. Buchgraber, “Variational Sparse Bayesian Learning: Centralized and Distributed Processing,” Ph.D. dissertation, Graz University of Technology, 2013.