\thesection Introduction
###### Abstract

Most existing approaches to co-existing communication/radar systems assume that the radar and communication systems are coordinated, i.e., they share information, such as relative position, transmitted waveforms and channel state. In this paper, we consider an un-coordinated scenario where a communication receiver is to operate in the presence of a number of radars, of which only a sub-set may be active, which poses the problem of estimating the active waveforms and the relevant parameters thereof, so as to cancel them prior to demodulation. Two algorithms are proposed for such a joint waveform estimation/data demodulation problem, both exploiting sparsity of a proper representation of the interference and of the vector containing the errors of the data block, so as to implement an iterative joint interference removal/data demodulation process. The former algorithm is based on classical on-grid compressed sensing (CS), while the latter forces an atomic norm (AN) constraint: in both cases the radar parameters and the communication demodulation errors can be estimated by solving a convex problem. We also propose a way to improve the efficiency of the AN-based algorithm. The performance of these algorithms are demonstrated through extensive simulations, taking into account a variety of conditions concerning both the interferers and the respective channel states.

\title

Adaptive Interference Removal for Un-coordinated Radar/Communication Co-existence \authorLe Zheng, \emphMember, \emphIEEE, Marco Lops, \emphSenior Member and Xiaodong Wang, \emphFellow, \emphIEEE \thanksLe Zheng and Xiaodong Wang are with Electrical Engineering Department, Columbia University, New York, USA, 10027, e-mail: le.zheng.cn@gmail.com, wangx@ee.columbia.edu.

Marco Lops is with the DIEI, Universita degli Studi di Cassino e del Lazio Meridionale, Cassino 03043, Italy (e-mail: lops@unicas.it, e.grossi@unicas.it).

\maketitle{IEEEkeywords}

Radar/communication co-existance, atomic norm, compressed sensing, off-grid, sparsity.

## \thesection Problem Formulation

We consider a situation with one communication system and active radars. Suppose the -th radar transmits the coded waveform

 sj(t)=N−1∑n=0gj(n)ξ(t−nT), (\theequation)

where is the -th code, is the code length, and satisfies the Nyquist criterion with respect to . From now on we assume that such a basic pulse is a Square Root Rased Cosine (SRRC) with excess bandwidth , so that the transmit bandwidth is . We also set with when is odd and when is even. Let be the Discrete Fourier Transform (DFT) of , i.e., with denoting the DFT matrix. In practice, the radar waveforms usually have some characteristics to achieve certain performance, so we assume lives in a low-dimensional subspace of , spanned by the columns of a known matrix with and , i.e., for some unknown such that for . An example of such a situation is waveform diversity, wherein radar systems may vary the transmit waveforms by selecting them in a set - a dictionary - so as to cope with interference, clutter, co-existence with competing systems [hassanien2016dual]. Obviously, if , the number of unknown variables exceeds , and the problem is for sure infeasible: in general, since the radar waveforms may contain a number of unknown parameters, the condition is mathematically necessary, on top of being per se reasonable. \colorblack For simplicity, we assume that there is one path between the radar TX and the communication RX. This assumption is true for narrow-band radar systems [li2016mimo] or as the interference is dominated by the direct path between the radar TX and the communication RX. The interference produced by active radars - with possibly unknown - onto the communication RX can be expressed as

 yI(t)=J∑j=1N−1∑n=0cjgj(n)ξ(t−nT−τj), (\theequation)

where and denote the delay and complex coupling coefficient of the -th radar, respectively. We assume the communication TX is not moving and its position is known by the radars, so its effect on the radar can be dealt with via beamforming. The communication TX transmits data symbols with . Let where denotes the set of possible values. Defining the data received at the RX side as , we have

 x=HAb, (\theequation)

where , and is the channel matrix. This model subsumes a number of communication systems. For example, in a Code-Division Multiple Access (CDMA) system [guo2005randomly], the elements in are the symbols transmitted by active users: The columns of are the signatures of the users, and is a diagonal matrix representing the channel gains. Another example is an OFDM system [zhang2015maximum], in which is the IDFT matrix with , is the channel matrix and is the received data in time domain. In practice, is known and the channel matrix can be obtained through the transmission of pilot signals. Let be the DFT of , i.e., . The received communication signal is given by

 yC(t)=N−1∑n=0x(n)ξ(t−nT−τC), (\theequation)

where denotes the overall delay of the communication transmission. As the communication TX and RX are synchronized, we let for the convenience of derivation. Notice also that the full bandwidth overlap scenario considered in this paper makes it reasonable to assume that the communication system occupies the full available bandwidth and uses the same SRRC basic pulse as the radar systems. At the communication RX, the signal contains both the communication signal and the radar interference, i.e.,

 y(t) = yI(t)+yC(t)+~w(t), (\theequation) = J∑j=1N−1∑n=0cjgj(n)ξ(t−nT−τj) +N−1∑n=0x(n)ξ(t−nT)+~w(t),

where is the measurement noise. Projecting onto results in

 r(t′) = ⟨y(t),ξ(t−t′)⟩ (\theequation) = N−1∑n=0x(n)Rξ(t′−nT) +J∑j=1N−1∑n=0cjgj(n)Rξ(t′−nT−τj)+w(t′),

where is the auto-correlation function of , i.e., with denoting inner product, and . The auto-correlation function is considered substantially time-limited in a finite interval, say: This is a common assumption which is justified by the fact that is always vanishingly small for large , and in particular goes to zero as in the considered scenario (see Appendix A). Letting and be the minimum and the maximum delays, respectively, tied to the corresponding minimal and maximum distances of all of the potential radar systems from the receiver, we define in \eqrefeq:truebarr

for , where . For simplicity, we define . It is assumed that is a complex Gaussian variable with zero mean and variance , i.e., . In communication receivers, is typically available and can be estimated off-line. Suppose that satisfies for . The above condition is rigorously true for small excess bandwith (in particular, ), which includes the relevant case that the communication system employs an OFDM format, corresponding to , while being only approximately true for larger values of [glover2010digital]. Since in an efficient spectrum exploitation context it is mandatory to choose small excess bandwiths, we henceforth assume

 ∫∞−∞Rξ(t−nT)e−i2πktNTdt≃e−i2πnkN. (\theequation)

Plugging \eqrefeq:integ into \eqrefeq:barr, we have

 ¯r(k) ≈ N−1∑n=0x(n)e−i2πnkN+¯w(k) (\theequation) +J∑j=1e−i2πkτ′jN−1∑n=0cjgj(n)e−i2πnkN, = ¯x(k)+¯w(k)+J∑j=1cj¯gj(k)e−i2πkτ′j,

where and are the -th element of and , respectively. We define . As outlined in the introduction, the communication RX has to remove the radar interference from the measurement with no knowledge of the delays and the waveforms of the active radars, i.e., under uncertainty concerning and . In principle, demodulation may be undertaken simply ignoring the presence of interference, i.e., through the operation , with the decoding function operating on the received signal, which would obviously lead to an uncontrolled symbol-error-rate (SER). The approach we take here instead relies on a joint interference-estimation symbol-demodulation process. In particular, defining as the “estimated communication signal”, the presence of errors in the decision process results into a non-zero difference vector where . Obviously, the -th element of is given by

 z(k) = ⟨HAv,fk⟩+J∑j=1cj¯gj(k)e−i2πkτ′j+¯w(k), (\theequation) = ⟨HAv,fk⟩+J∑j=1cja(τ′j)Hek¯dHkhj+¯w(k), = ⟨HAv,fk⟩+⟨X,¯dkeHk⟩+¯w(k),

where we have defined , with , and the -th column of the identity matrix . Once estimates of and are available, the radar interference can be obtained and canceled from the measurements and the symbols re-demodulated. Hence, the main problem is to estimate and from the noisy measurements . Notice that contains both the radar interference and the residual of communication signal caused by the mis-demodulations. The mixing of both signals causes great difficulties for the estimation, which inspires us to exploit some structural information about the desired solution, and in particular sparsity, as detailed in the next section.

## \thesection Proposed Algorithms

Equation \eqrefeq:z highlights that data demodulation and interference mitigation are coupled, in the sense that they should be accomplished jointly and that poor performance in estimating either one has detrimental effects on the estimate of the other. In fact, in order to remove the radar interference, we need to estimate the matrix from the observations \eqrefeq:z, but this would obviously require also estimating the error vector , which boils down to correctly demodulating the data block. To this end, we design iterative algorithms, exploiting structural information on the desired solution. In the -th iteration, the demodulated symbols are denoted as , and where is the estimate of . \colorblack The following two types of sparsity are exploited in the problem:

1. The signal is a combination of complex exponentials with unknown modulation , and the number of complex exponentials is much smaller than the dimension of , i.e., .

2. Ideally, the vector should be an all-zero vector. As a consequence, denoting the result of the -th iteration, we want to be as small as possible, and in any case we want to force the condition .

### \thesubsection Joint Waveform Estimation and Demodulation Based on On-grid CS Algorithm

In the first iteration, so . According to \eqrefeq:z, jointly estimating , and is a non-linear problem. However, it can be linerized by using an overcomplete dictionary matrix

 ~A=[a(~τ′1),a(~τ′2),...,a(~τ′~J)]∈CN×~J, (\theequation)

with denoting sets of uniformly spaced points of the radar delays. Define as the number of columns of where . For sufficiently large , the delay is densely sampled. Let be the sparse vector whose non-zero elements correspond to in \eqrefeq:z. As usual, forcing a constraint onto the -norm is impractical, since it results in an NP-hard non-convex optimization problem, and -norm regularization is used instead, i.e., and . Define as the set of all possible differences when the two vectors both belong to . We notice that the constraint results in a non-convex problem, and would cause much difficulty in solving the optimization problem. This constraint is relaxed and the non-linear joint estimation problem is thus reduced to a linear parameter estimation problem, i.e., the estimation of the linear amplitude vectors and , under a sparsity constraint:

 (^α(1),^v(1)) = argminα∈C~JK×1,v(1)∈CN×112∥∥z−Φv(1)−Υα∥∥22 (\theequation) +~λ∥α∥1+~γ∥v(1)∥1,

where

 Φ=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣fH−~N1HAfH−~N1+1HA⋮fH~N2HA⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (\theequation)

and is given by \eqrefeq:Upsilon.

The parameters and are weights determining the sparsity of the reconstruction. In practice, we set . As \eqrefeq:cs is convex, it can be solved with standard convex solvers. The computational complexity of problem solving depends on the dimension of and . Specifically, if \eqrefeq:cs is solved by using the interior point method [boyd2004convex], the complexity in each iteration is . By solving \eqrefeq:cs, we obtain the estimates and , whereby the symbols can be corrected and re-demodulated as:

 ^b(1)=argminb∈B∥∥∥b−^b(0)−^v(1)∥∥∥2. (\theequation)

Notice that the re-demodulation process makes use of the structural information of the communication symbols, i.e., . After the re-demodulation, the demodulated symbols belong to the constellation alphabet, whereby . When the estimates of the symbols in the first iteration is accurate and the measurement noise is small, i.e., and are both small, then all the mistakenly demodulated symbols can be corrected by applying \eqrefeq:estb. In many cases, however, the interference from the radars is strong, and contains many demodulation errors. Thus, we need to iterate the joint interference removal/data demodulation process. Specifically, in the -th iteration (), and are estimated as

 (^α(l),^v(l)) = argminα∈C~JK×1,v(l)∈CN×1~λ∥α∥1+~γ∥v(l)∥1 (\theequation) +12∥∥z(l)−Φv(l)−Υα∥∥22.

Then re-demodulation is undertaken by \eqrefeq:estb with replaced by and replaced by . As the iteration goes, some wrong symbols are corrected, and the demodulation error becomes sparser as increases. The proposed algorithm iterates until or the maximum number of iterations is reached. Let be the estimate of when the algorithm terminates. The radar delays can be identified by locating the non-zero entries of . If the solution is either non-zero or has elements larger than a pre-set threshold, i.e., , then a radar interference exists at delay . Notice that one cannot resolve the inherent scaling ambiguity between each and the corresponding , which is in any case not essential since it is the product the measure of interest for radar interference removal purposes. The estimated time domain radar waveform is then given by

 ^cj^gj=FH¯D~cj~hj. (\theequation)

The CS algorithm based on -minimization (CS-L1) is capable of super-resolving the spectrum of the sparse signal under certain conditions of the matrices and [eldar2009robust]. For comparison purposes, and to quantify, at the performance assessment stage, the loss induced by the uncertainty on the active radars positions, it is worth exploring also the case that the delays of such active radar systems are known at the communication receiver. In such a simplified scenario, which assumes that training pilot signals are periodically transmitted for channel sensing purposes, a modified version of the proposed algorithm can be easily derived, by using arguments similar to those employed above. Indeed, since the radar delays are known precisely, we let , and \eqrefeq:z can be re-written as

 z(k)=⟨HAv,fk⟩+ϕk¯α+¯w(k), (\theequation)

where . The algorithm with known timing and radar delay has the same procedure as that of CS-L1. Notice that is not sparse. In the -th iteration, and the demodulation error can be estimated by solving \eqrefeq:cs2 with , removed and replaced by .

### \thesubsection Joint Waveform Estimation and Demodulation Based on Off-grid CS Algorithm

As anticipated, \eqrefeq:z reveals that is a linear combination of modulated complex exponentials with arbitrary phases, which do not in general correspond to the point of a discrete grid: the off-grid radar position can lead to mismatches in the model and deteriorate the performance. In this subsection, we use the atomic norm to build a sparse representation which does not suffer from the off-grid problem. We define the atomic norm [yang2016super] associated to as

 ∥X∥A = inf{μ>0:X∈μconv(A)} = infcj,τ′j,∥hj∥2=1{∑j|cj|:X=∑jcjhja(τ′j)H},

where denotes the convex hull of the input atom set, and the set of atoms is defined as

 A={ha(τ′)H:τ′∈[0,1),∥h∥2=1,h∈CK×1}. (\theequation)

For future developments, we introduce the following equivalent form of the atomic norm for the atom set [yang2016super]:

 ∥X∥A=infu,T⎧⎪ ⎪⎨⎪ ⎪⎩12NTr(Toep(u))+12Tr(T),s.t.[Toep(u)XHXT]⪰0⎫⎪ ⎪⎬⎪ ⎪⎭, (\theequation)

where is a complex vector whose first entry is real, denotes the Hermitian Toeplitz matrix whose first column is , and is a Hermitian matrix. Based on \eqrefeq:z, and paralleling the arguments outlined in the previous sub-section, the -th iteration achieves estimates, and say, of and by processing and solving the optimization problem:

 (^X(l),^v(l))=minX,v(l)λ∥X∥A+γ∥v(l)∥1 +~N2∑k=−~N112∣∣z(l)(k)−⟨fk,HAv(l)⟩−⟨X,¯dkeHk⟩∣∣2,

where and are the weight factors. In practice, we set and . In the light of \eqrefeq:atomic, the above can be transformed into the following Semi-Definite Programming (SDP):

 (^X(l),^v(l))=argminX,T,u,v(l)λ2NTr(Toep(u)) +~N2∑k=−~N112∣∣z(l)(k)−⟨fk,HAv(l)⟩−⟨X,¯dkeHk⟩∣∣2 +λTr(T)2+γ∥v(l)∥1, (\theequation) s.t.[Toep(u)XHXT]⪰0,

where denotes the Toeplitz matrix whose first column is the input vector. The above problem is convex, and can be solved by using a convex solver. The corresponding computational load is per iteration if the interior point method is applied. We name the algorithm based on solving \eqrefeq:atomicnorm as the CS Atomic-Norm (CS-AN) based algorithm. Similar to the CS-L1 algorithm, CS-AN iterates \eqrefeq:SDP and \eqrefeq:estb until . From now on, denotes the estimate of when the algorithm terminates. We hasten to underline here that the atomic norm can be exploited to enforce sparsity in the continuous domain without any discretization [tang2013compressed, yang2016super], thus leading to better estimation performance compared to on-grid CS techniques. Atomic-norm-based optimization has also been shown to achieve better de-noising performance than subspace methods such as MUSIC [bhaskar2013atomic]. Notice that if has low rank, it is amenable to nuclear-norm-based estimation. However, the nuclear norm-based approach cannot fully take advantage of the signal structure (e.g., in our scenario is a combination of multiple complex exponentials with unknown modulation ), whereby its performance has been shown to be inferior to that of atomic norm-based optimization [fernandez2016demixing, chen2014robust]. Solving \eqrefeq:SDP does not directly provide estimates of the delays of the active radars, . Notice however that each row of is a linear combination of several complex exponentials, in that, denoting the -th row of , we have . Hence, and can be obtained by MUSIC [naha2015determining] or prony’s method [stoica1997introduction] with as input. Denoting the operation of the MUSIC algorithm, it outputs components with different amplitudes and delays111Here is not necessarily the same as because can be zero while for some . For example, the radar has candidate waveforms, and select one for transmission. In such case, is a vector having zero elements and one element with magnitude 1. As a result, the number of output delays should satisfy ., i.e.,

 Tk={¯βj(k),¯τ′j(k)}j=1,2,...,¯Jk=MUSIC(^Xk,1:N) (\theequation)

for . An association procedure is thus needed to combine the different s and come up with the estimates, say, of the radar signals delays. In practice, the calculation of may not be accurate due to computational errors. Hence, if the estimated delays for different are closer than a small threshold , then they are regarded as generated by the same radar, and the corresponding s are combined. For clarity, we summarize the association process in Algorithm 1, where is the set of the estimated radar delays and waveform parameters. In the algorithm, for , if there exists such that and , then both components belongs to the same radar and the radar delay is updated by the weighted summation of and . Otherwise, an additional radar with the estimated delay and amplitude is added to the set . Notice that, once again, the inherent scaling ambiguity between each and the corresponding . Hence, we only estimate via \eqrefeq:hatg with replaced by . {algorithm}[htbp] Radar delay and path gain estimation Input for , . 1, , . For      For         If there exists such that         and             2,