Interference Removal for Radar/Communication Co-existence: the Random Scattering Case

1 Interference Remo v al for Radar/Communication Co-e xistence: the Random Scattering Case Y inchuan Li, Le Zheng, Member , IEEE , Marco Lops, F ellow , IEEE and Xiaodong W ang, F ellow , IEEE Abstract In this paper we consider an un-cooperati ve spectrum sharing scenario, wherein a radar system is to be ov erlaid to a pre-existing wireless communication system. Given the order of magnitude of the trans- mitted po wers in play , we focus on the issue of interference mitigation at the communication recei v er . W e explicitly account for the reverberation produced by the (typically high-po wer) radar transmitter whose signal hits scattering centers (whether targets or clutter) producing interference onto the communication receiv er , which is assumed to operate in an un-synchronized and un-coordinated scenario. W e first show that receiv er design amounts to solving a non-con ve x problem of joint interference remov al and data demodulation: next, we introduce two algorithms, both e xploiting sparsity of a proper representation of the interference and of the vector containing the errors of the data block. The first algorithm is basically a relaxed constrained Atomic Norm minimization, while the latter relies on a tw o-stage processing structure and is based on alternating minimization. The merits of these algorithms are demonstrated through e xtensiv e simulations: interestingly , the two-stage alternating minimization algorithm turns out to achie ve satisf actory performance with moderate computational complexity . Index T erms Y . Li is with the School of Information and Electronics, Beijing Institute of T echnology , Beijing 100081, China, the Beijing Ke y Laboratory of Embedded Real-time Information Processing T echnology , Beijing 100081, China and Electrical Engineering Department, Columbia Univ ersity , New Y ork, USA, 10027 (e-mail: yinchuan.li.cn@gmail.com). L. Zheng and X. W ang are with Electrical Engineering Department, Columbia Univ ersity , New Y ork, USA, 10027 (e-mail: le.zheng.cn@gmail.com; wangx@ee.columbia.edu). M. Lops is with the Department of Electrical Engineering and Information T echnologies, Univ ersit ` a di Napoli ”Federico II”, V ia Claudio, 21 - I-80125 Naples (Italy) (e-mail: lops@unina.it). 2 Radar/communication co-e xistence, multi-path, atomic norm, compressed sensing, non-con ve x, blind decon volution, off-grid, sparsity . I . I N T RO D U C T I O N The e ver increasing demand for high data rates in wireless communications has forced co- existence of communication and radar systems in the same frequency bands [1]: this can be achie ved by either allo wing only one system to be equipped by an acti ve properly designed transmitter - see, e.g. the information embedding strategies [2] to transmit information through a radar wav eform, and the approach in [3], [4], which can somehow be classified as a passi ve radar [1], to accomplish sensing functions through communication signals - or considering architectures with multiple transmitters operating in spectral overlap [5]–[7]. The latter scenario, which is the one considered in this paper , requires proper transcei ver design: the strategies proposed so far range from a geometrical approach, aimed at mitigating the interference produced by one system on the other through suitable projection operations [8], [9], to a cognition-based radar wav eform design [10]–[12]. A more comprehensiv e approach is co-design [13]–[15], wherein the radar wav eform(s) and the communication code-book are jointly designed by minimizing a measure of the mutual interference under certain constraints. A common point of these strategies is some form of coordination between the two activ e systems, and a remarkable degree of prior cognition, to be possibly acquired or updated through the periodic transmission of pilot signals to handle dynamic scenarios. In some situations, howe ver , such a cooperation is either un-feasible - due, e.g., to security reasons - or too costly , whereby the radar and the communication systems should operate with little or no coordination. Such scenarios hav e been considered, e.g., in [16], wherein a blind null space estimation method is proposed as an extension of the results of [17]. A dif ferent approach to handle un-coordinated co-existence is the one proposed in [18], considering full bandwidth ov erlap between a pre-e xisting communication system and multiple o verlaid radars: assuming that the interfering radar wa veforms li ve in the subspace of a kno wn dictionary , the communication performance is guaranteed by joint interference remov al/data demodulation iterati ve procedures, le veraging ideas from compressed sensing and atomic norm (AN) minimization techniques. A major limitation of [18] is that the clutter induced by random scatterers disseminated in the con- trolled scene and reflecting the radar signal towards the communication recei ver is not accounted for: this is a signal-dependent interference which, if not properly handled, typically produces 3 dramatic ef fects on the radar performance and could totally pre vent reliable communication. Additionally , synchronism between the radar and the communication system is assumed, as well as prior knowledge of the afore-mentioned dictionary . The present contribution is aimed at extending the results of [18] by explicitly accounting for the re verberation produced by a single radar transmitter onto the communication receiv er . In particular , we consider an Orthogonal Frequenc y Division Multiplexing (OFDM) communication system co-e xisting with a short-range radar using a sophisticated wa v eform: the Pulse Repetition Interv al (PRI) of the radar coincides with the duration of the communication data symbol block, and a totally un-synchronizated and un-coordinated scenario is considered, nor any assumption is made on the radar code structure. It is note worthy that the PRI of the radar coincides with the duration of the communication data symbol block is possible in practice, as detailed in the next section. Similar to [18], we focus on the communication receiv er performance, which is justified in the light of sev eral considerations: first, the order of magnitudes of the po wers transmitted by the communication and the radar transmitters is typically very different; additionally , while the communication transmitter points at the communication receiv er whose location is typically known, whereby its effect on the radar receiver can be mitigated through beam-forming techniques [19], a search radar employs rather wide and rotating beams, which produce random and time-v arying re verberation onto the communication receiv er . T o this end, we propose two dif ferent algorithms, both exploiting two types of sparsity: on one hand, indeed, as scatterers are sparsely distributed in space, the interfering signals hitting the communication RX are sparse ; on the other , an iterati ve demodulation algorithm should require that the vector containing the demodulation errors of a data block be itself sparser and sparser as the iterations go. Since the delays with which the interferers arriv e at the communication receiv er are con- tinuous parameters, mere application of compressed sensing theory [20], [21] would produce unsatisfactory performance [22] in a situation where these signals cannot be sparsely represented by a finite discrete dictionary [23]–[25]. W e consider instead the recently de veloped mathematical theory of continuous sparse recovery for super-resolution [26]–[28], and especially of the AN minimization techniques which are successfully used for continuous frequency recovery , line spectral estimation and direction-of-arriv al estimation [28]–[30]. As an alternati ve, based on the fact that the radar code is unknown and the radar interferences impinge on the communication RX with unknown multiple delays and coupling coefficients, estimating the interfering code and the multiple delays is inherently linked to solving a blind decon volution problem [31]–[33], 4 C o m m u n i c a t i o n s i g n a l C y c l i c p r e f i x R a d a r s i g n a l D a t a C o m p l e t e P R I s i g n a l o v e r l a p p i n g c y c l i c s h i f t U n - s y n c h r o n i z e d e r r o r P R I = P R I = Fig. 1. Transmitted radar and communication signals which is non-con ve x and ill-posed without further constraints: this moti v ated us to also explore the recently dev eloped mathematical theory of blind decon volution [33]–[35] to improv e the estimation accuracy of the interfering wa veform. The remainder of this paper is or ganized as follows. In Section II, we present the signal models of the co-existed radar and communication system and set up the problem. In Section III, we dev elop the proposed con ve x relaxation method using both the AN and the ` 1 -norm. In Section IV , the proposed two-stage alternating minimization algorithm is dev eloped. Simulation results are presented in Section V . Finally , in Section VI, we draw conclusions from the results obtained in this paper . I I . S Y S T E M D E S C R I P T I O N S & P RO B L E M F O R M U L A T I O N A. T r ansmitted Signals W e consider an OFDM communication system coexisting with a radar system. Assume that the OFDM system consists of N = N d + N p sub-carriers, with N d data sub-carriers and N p cyclic prefix (CP) sub-carriers. The duration of an OFDM block is N T , T being the “sub-pulse duration. ” Denote the n c -th normalized data symbol block as b n c ( k ) , k = 0 , ..., N d − 1 , such that 5 E [ b n c ( k ) b n c ( k ) ∗ ] = 1 with ( · ) ∗ denoting the complex conjugate operator . Then, the transmitted baseband OFDM signal is given by s c ( t ) = ∞ X n c = −∞ N d − 1 X k =0 b n c ( k ) e i 2 π k t N d T u c ( t − n c N T ) , (1) where u c ( t ) =    1 , t ∈ [ − N p T , N d T ] , 0 , otherwise . (2) As for the radar signal, we assume that the communication and the radar systems are in full bandwidth overlap , and the PRI equals the duration of data symbols N d T 1 . This assumption implies that, as sho wn in Fig. 1, the sub-pulse duration of the radar system and that of the communication system are the same; and at each PRI, a block of N d data symbols are transmitted. On the other hand, when the communication receiver processes one OFDM block of N d symbols, there is a complete PRI radar signal ov erlapping therewith after the cyclic shift, regardless of whether or not the two systems are synchronized. W e assume that the radar transmits a single sophisticated (i.e., with large duration-bandwidth product) pulse in any gi ven PRI, which consists of L amplitude-modulated sub-pulses. Denoting by g = [ g (0) , g (1) , ..., g ( L − 1)] T ∈ C L × 1 the wa veform code and by ξ ( t ) the basic sub-pulse wa veform, the transmitted baseband radar signal is giv en by s r ( t ) = ∞ X n r = −∞ L − 1 X ` =0 g ( ` ) ξ ( t − `T − n r N d T ) , (3) where T << PRI is approximately the in verse of the bandwidth and is related to the radar range resolution. W e remind here that the duty cycle δ = LT PRI  1 is typically lo w in order to guarantee a proper hearing period [39]. 1 This assumption is possible in practice. For example, according to the December 2017 3GPP first release of the 5G Ne w Radio standard, the data symbols of the 5G signal hav e a duration on the order of 10 µs , while short-range civilian radars (e.g., automotiv e radar) typically hav e a PRI in the order of 10 µs [36]. In addition, some WLAN systems use OFDM wav eform with a data symbol duration on the order of 1 µs [37], while short-range impulse radars for high speed moving targets detection and through-the-wall radars may have a PRI in the order of 1 µs [38], [39]. 6 B. Received Signal W e assume that the communication system operates on a block-fading channel whose coher- ence time is much larger than the OFDM blocklength, whereby the useful component at the communication receiver is giv en by y c ( t ) = s c ( t ) ∗ h ( t ) = s c ( t ) ∗ M c X m =1 α m δ ( t − τ c m ) , (4) In the previous equation ∗ denotes the con volution operator , h ( t ) is the channel impulse response, M c is the total number of propagation paths, α m and τ c m are the m -th path’ s complex gain and delay , respectiv ely . The presence of a co-e xisting radar system produces additional interference on the communi- cation receiv er . In particular , if we assume that there are M r scatterers, whether from clutter or targets, located in as many different range cells, the signal scattered to wards the communication recei ver can be modeled as y r ( t ) = s r ( t ) ∗ M r X m =1 c m e j 2 π f m t δ ( t − τ R − τ r m ) , (5) where, since the radar and communication systems are un-synchronized, 0 ≤ τ R ≤ N d T is the corresponding delay at a reference interv al (i.e., for n c = n r = 0 ), while c m , τ r m and f m denote the scattering coefficient, the delay and the Doppler shift of the m -th reflector , respectiv ely . On the recei ver side, we assume that the communication receiv er processes one OFDM block of N d symbols at a time. Since the duration of an OFDM block is usually small, we hav e f m N d T  1 , then the phase rotation due to the Doppler shift ov er a block duration can be approximated as constant [40], and is thus not measurable and uninfluential, hence it is ignored from now on. The CP is removed assuming that its length is no less than the maximum communication multi-path delay . Let n c = 0 with no loss of generality and thus the subscript n c is also omitted. Focusing the attention on the interval [0 , N d T ] , we thus have, for the receiv ed signal, the model: r ( t ) = ∞ X n r = −∞ M r X m =1 c m L − 1 X ` =0 g ( ` ) ξ ( t − `T − n r N d T − τ R − τ r m ) + M c X m =1 α m N d − 1 X k =0 b ( k ) e i 2 π k t − τ c m N d T + ˜ w ( t ) , t ∈ [0 , N d T ] , (6) where ˜ w ( t ) is a white, complex circularly symmetric Gaussian noise process. 7 C. Pr oblem F ormulation The communication receiv er is assumed to undertake the standard OFDM operations on each OFDM packet of duration N d T . In particular , we focus on the first packet occupying the interv al [0 , N d T ] [41]. Let ∆ m ,  − τ R − τ r m N d T  , (7) τ m , − τ R − τ r m N d T − ∆ m ∈ [0 , 1) , (8) where b·c is the floor function. W e hav e ¯ r ( k ) for k = 0 , ..., N d − 1 , ¯ r ( k ) = 1 N d T Z N d T 0 r ( t ) e − i 2 πk t N d T dt = M r X m =1 c m L − 1 X ` =0 g ( ` ) 1 N d T Z N d T 0 ∞ X n r = −∞ ξ ( t − `T − n r N d T + τ m N d T + ∆ m N d T ) e − i 2 πk t N d T dt + M c X m =1 α m N d − 1 X n =0 b ( n ) 1 N d T Z N d T 0 e i 2 π n t − τ c m N d T e − i 2 πk t N d T dt + w ( k ) = M r X m =1 c m e i 2 π kτ m L − 1 X ` =0 g ( ` ) e − i 2 πk ` N d 1 N d T Z N d T − τ m N d T − `T − τ m N d T − `T ∞ X n r = −∞ ξ ( t − ( n r − ∆ m ) N d T ) e − i 2 πk t N d T dt | {z } ¯ ξ ( ω ) | ω = 2 πk N d T + M c X m =1 α m e − i 2 π k τ c m N d T N d − 1 X n =0 b ( n ) 1 N d T Z N d T 0 e i 2 π ( n − k ) t N d T dt | {z } N d T · δ ( n − k ) + w ( k ) (9) = M r X m =1 c m e i 2 π kτ m L − 1 X ` =0 g ( ` ) e − i 2 πk ` N d | {z } ¯ g ( k ) ¯ ξ  2 π k N d T  + M c X m =1 α m e − i 2 π k τ c m N d T | {z } H ( k ) b ( k ) + w ( k ) , (10) where we note that in (9) P ∞ n r = −∞ ξ ( t − ( n r − ∆ m ) N d T ) is a periodic signal with period N d T , and each period is composed of ξ ( t ) , t ∈ [0 , T ] ; therefore the first integral is its Fourier transform, i.e., ¯ ξ ( ω ) = 1 N d T Z N d T 0 ξ ( t ) e − iω t dt (11) e valuated at ω = 2 π k N d T . In (10), w ( k ) , 1 N d T Z N d T 0 ˜ w ( t ) e − i 2 πk t N d T dt ∼ C N (0 , σ 2 w ); (12) ¯ g = [ ¯ g (0) , ¯ g (1) , ..., ¯ g ( N d − 1)] T = F L g ∈ C N d × 1 (13) 8 is the discrete Fourier transform (DFT) of g , with F L denoting the first L columns of the N d -points DFT matrix F ; and H ( k ) = Z N d T 0 h ( t ) e − i 2 πk t N d T dt = M c X m =1 α m e − i 2 π k τ c m N d T (14) is the channel transfer function at frequency k N d T , which can be estimated using pilot signals [42]. Let us now define ¯ ξ = [ ¯ ξ (0) , ¯ ξ ( 2 π N d T ) , ..., ¯ ξ ( 2 π ( N d − 1) N d T )] T ∈ C N d × 1 (15) and H = diag([ H (0) , H (1) , ..., H ( N d − 1)] T ) ∈ C N d × N d , (16) i.e., an N d × N d diagonal matrix with elements of [ H (0) , H (1) , ..., H ( N d − 1)] T on the diagonal. W e also introduce the vectors ¯ r = [ ¯ r (0) , ¯ r (1) , ..., ¯ r ( N d − 1)] T ∈ C N d × 1 , (17) b = [ b (0) , b (1) , ..., b ( N d − 1)] T ∈ C N d × 1 , (18) w = [ w (0) , w (1) , ..., w ( N d − 1)] T ∈ C N d × 1 (19) and ν τ = M r X m =1 c m a ( τ m ) ∈ C N d × 1 , (20) with a ( τ ) = [1 , e i 2 π τ , ..., e i 2 π ( N d − 1) τ ] T . (21) Then, (10) can be given the following compact vector form ¯ r = H b + ¯ ξ  ( F L g )  ν τ + w , (22) where  denotes the pointwise product. Assume that an estimate of the data symbols, ˆ b , is av ailable by directly performing demodulation using ¯ r . W e subtract the demodulated data from ¯ r , to obtain z = ¯ r − H ˆ b = H v + ¯ ξ  ( F L g )  ν τ + w , (23) where v = b − ˆ b ∈ C N d × 1 . (24) 9 Our main problem is to estimate g , ν τ and v from the noisy measurements z . T o this end, we first notice that, in a realistic scenario, the number of scatterers is much lo wer than the number of OFDM symbols in a packet, i.e. M r  N d in (20) ; secondly , we want the demodulation error rate to be lo w , i.e. we want to force the vector v to have a small number of non-zero entries: both are sparsity conditions that we can exploit. Notice howe ver that the delays τ m in (20) take on continuous values, whereby using traditional compressed sensing techniques would entail heavy losses due to the off-grid problem: as a consequence, we resort here to Atomic Norm (AN) minimization instead [28], [43]. Conv ersely , the second type of sparsity simply results in a suitable constraint in the optimization problem. T o be more precise, define the set of atoms A = { a ( τ ) : τ ∈ [0 , 1) } . Then the ` 0 -atomic norm [44] associated to ν τ is given by k ν τ k A , 0 = inf c m ∈ C ,τ m ∈ [0 , 1) ( M : ν τ = M X m =1 c m a ( τ m ) ) . (25) Our problem can be formulated as ( ˆ g , ˆ ν τ , ˆ v ) = arg min g ∈ C L × 1 , ν τ ∈ C N d × 1 v ∈ C N d × 1 k ν τ k A , 0 + λ k v k 0 , (26) s.t.   z − H v − ¯ ξ  ( F L g )  ν τ   2 2 ≤ , k g k 2 = 1 , where λ > 0 is a weight factor ,  > 0 is the error tolerance and k v k 0  N d denotes the ` 0 -norm of v . For the case that the radar signal is strong, we can perform iterati ve demodulation and radar interference estimation: in each iteration, after solving (26), we make use of ˆ v and the current ˆ b to obtain a refined demodulation ˜ b = arg min b ∈B N d k b − ˆ b − ˆ v k 2 , (27) where B is the modulation symbol constellation set. Then we update z in (23) by setting ˆ b ← ˜ b and solve (26) again. Note that in (26) the objectiv e function is non-con ve x since it in volv es the ` 0 -atomic norm and the ` 0 -norm. The first constraint is also non-con vex, because ( F L g )  ν τ is the DFT of the con v olution g ~ ( F − 1 ν τ ) with ~ the circular con volution operator, and it is kno wn that the blind decon v olution problem is non-con v ex [31]–[33]. 10 I I I . T H E C O N V E X R E L A X A T I O N M E T H O D Define D = diag( ¯ ξ ) F L ∈ C N d × L , and let d H k ∈ C 1 × L be the k -th row of D . Then, (23) can be rewritten as z ( k ) = e T k ( H v ) + d H k g e T k ν τ + w ( k ) = e H k ( H v ) + d H k ( g ν T τ ) e k + w ( k ) = h H v , e k i +  g ν T τ , d k e H k  + w ( k ) , k = 0 , ..., N d − 1 , (28) where z ( k ) denotes the k -th element of z , e k is the k -th column of the N d × N d identity matrix and h X , Y i = T r ( Y H X ) . Notice that the original problem in (26) entails estimating g and ν τ separately . In the new formulation, we are interested in estimating g ν T τ = g P M r m =1 c m a ( τ m ) T instead. In particular , we relax (28) by introducing X = M r X m =1 c m g m a ( τ m ) T ∈ C L × N d , (29) i.e., a mixture of atoms from the atom set ˜ A =  g a ( τ ) T : τ ∈ [0 , 1) , k g k 2 = 1 , g ∈ C L × 1  . (30) as the quantity of interest. Further , we replace the ` 0 -atomic norm and the ` 0 -norm in the objectiv e function of (26) by the ` 1 -atomic norm [43] and the ` 1 -norm, respectiv ely . The ` 1 -atomic norm seeks the tightest con v ex relaxation of enforcing sparsity in the atom set ˜ A , and is defined as k X k ˜ A , 1 = inf n η > 0 : X ∈ η con v ( ˜ A ) o = inf c m ∈ C ,τ m ∈ [0 , 1) k g m k 2 =1 ( X m | c m | : X = X m c m g m a ( τ m ) T ) , (31) where con v ( · ) denotes the con v ex hull of the input atom set. It is kno wn that (31) has the follo wing equi v alent form [43]: k X k ˜ A , 1 = inf u ∈ C N d × 1 , T ∈ C L × L          1 2 N d T r(T o ep( u )) + 1 2 T r( T ) , s . t .   T o ep( u ) X H X T    0          , (32) where u ∈ C N d × 1 is a complex vector whose first entry is real, T o ep( u ) denotes the N d × N d Hermitian T oeplitz matrix whose first column is u , and T is a Hermitian L × L matrix. Equations (31) and (32) are related through the relationship T o ep( u ) = X m | c m | a ( τ m ) a ( τ m ) H , (33) T = X m | c m | g m g H m . (34) 11 Finally , we can relax the original problem in (26) to the following semi-definite programming (SDP) [18], [33], [45] ( ˆ X , ˆ v ) = arg min X ∈ C L × N d , T ∈ C L × L u ∈ C N d × 1 , v ∈ C N d × 1 T r (T oep( u )) 2 N d + T r ( T ) 2 + ¯ λ k v k 1 , (35) s.t. N d − 1 X k =0   z ( k ) − h H v , e k i −  X , d k e H k    2 ≤ ,   T o ep( u ) X H X T    0 , where ¯ λ > 0 is a weight factor . Since problem (35) is con ve x, it can be solved with standard con v ex solvers, e.g., CVX [46]. Once an estimate ˆ X of X is obtained, estimates of the delays { τ m } and of the radar code g can be obtained by either solving the dual problem of (35) as in [18], [33], or using the MUltiple SIgnal Classifier (MUSIC) method as in [47]. Note that by relaxing the original problem of estimating g ν T τ = g P M r m =1 c m a ( τ m ) T to estimating X = P M r m =1 c m g m a ( τ m ) T , we may obtain spurious scatterers in solving the relaxed problem. As an example, suppose that the true code is g = [ 1 √ 2 , 1 √ 2 , 0 , ..., 0] T ∈ C L × 1 , and there are two scatters with τ 1 = τ ∈ [ 1 N d , 1 − 1 N d ) , c 1 = √ 2 , τ 2 = τ − 1 N d , c 2 = √ 2 . (36) Then the following is a spurious solution to the relaxed problem: g 1 = [ 1 √ 3 , 1 √ 3 , 1 √ 3 , 0 , ..., 0] T ∈ C L × 1 , g 2 = [0 , 0 , 1 , 0 , ..., 0] T ∈ C L × 1 , (37) τ 0 1 = τ , c 0 1 = √ 3 , τ 0 2 = τ + 1 N d , c 0 2 = 1 , in the sense that we hav e  g ν T τ , d k e H k  = ¯ ξ ( k )  e i 2 π kτ + 2 e i 2 π k ( τ − 1 N d ) + e i 2 π k ( τ − 2 N d )  =  X , d k e H k  , k = 0 , ..., N d − 1 , (38) where ν τ = c 1 a ( τ 1 ) + c 2 a ( τ 2 ) , X = c 0 1 g 0 1 a ( τ 0 1 ) T + c 0 2 g 0 2 a ( τ 0 2 ) T . (39) 12 I V . A T W O - S T A G E A LT E R N A T I N G M I N I M I Z A T I O N A L G O R I T H M Here we propose a new method to solve the non-con ve x problem in (26). The basic idea is to alternativ ely solve with respect to (w .r .t) g and ( ν τ , v ) ; and in solving w .r .t. g , we use the conjugate gradient search on Riemannian manifold; in solving w .r .t. ( ν τ , v ) , we take the matching-pursuit and greedy demixing approach. Moreov er , we solve the problem in (26) twice in a two-stage fashion: the first stage obtains a local optimum and the second stage makes use of the first-stage solution in forming the initial condition and solves a higher dimensional problem that leads to an approximate global optimum. A. Stage 1 - Obtaining Local Optimum W e first obtain a locally optimal solution to the non-con ve x problem (26) by iterati vely solving w .r .t. g and ( ν τ , v ) as follows: S-1: Let ˆ g be the estimate - a vailable from the previous iteration - of g , and define Φ , diag( ¯ ξ  ( F L ˆ g )) . (40) Then the new estimates ( ˆ ν τ , ˆ v ) can be obtained by solving the problem: ( ˆ ν τ , ˆ v ) = arg min ν τ ∈ C N d × 1 v ∈ C N d × 1 k ν τ k A , 0 + λ k v k 0 , s.t. k z − H v − Φ ν τ k 2 2 ≤ . (41) S-2: W ith the estimates ˆ ν τ and ˆ v of the previous step, defining ¯ z , z − H ˆ v , (42) W , diag( ¯ ξ  ˆ ν τ ) F L . (43) g can be easily updated by solving: ˆ g = arg min g ∈ C L × 1 k ¯ z − W g k 2 2 , s.t. k g k 2 = 1 . (44) The abov e alternating minimization procedure can be initialized by a random radar code ˆ g . W e next present the details of the two steps. 1) Greedy-demixing for solving S-1: Since ν τ = P M r m =1 c m a ( τ m ) , estimating ν τ in (41) im- plies estimating M r as well as tw o v ectors c = [ c 1 , c 2 , ..., c M r ] T ∈ C M r and τ = [ τ 1 , τ 2 , ..., τ M r ] T ∈ [0 , 1) M r in ν τ = Θ ( τ ) c , where Θ ( τ ) = [ a ( τ 1 ) , a ( τ 2 ) , ..., a ( τ M r )] ∈ C N d × M r . (45) 13 If the delays are on-grid, it is easy to estimate v and τ in (41) using an orthogonal matching pursuit (OMP) algorithm [48]. Howe ver , since the delays here are of f-grid, step S-1 in volv es not only demixing v and ν τ , but also super-resolving the delays in ν τ . In the spirit of the matching- pursuit method [49] and the greedy-demixing approach of [50], we adopt the following procedure for solving (41) in step S-1. S-1(a) Initialization: Let S be the set of support of v , and T be the set of delays, and initialize them as S ← ∅ , T ← ∅ . Define a residual r res ∈ C N d × 1 and initialize it as r res ← z . S-1(b) Selection: Find the highest correlation with the current residual r res and update either S or T . In particular , compute k  = arg max k ∈{ 0 , 1 ,...,N d − 1 } |h H (: , k ) , r res i| , (46) τ  = arg max τ ∈ [0 , 1) |h Φ a ( τ ) , r res i| . (47) If ˜ λ |h H (: , k  ) , r res i| > |h Φ a ( τ  ) , r res i| , then S ← S ∪ { k  } otherwise T ← T ∪ { τ  } , where ˜ λ is a weight factor . T o compute (47), we first search ov er a fine grid on [0 , 1) with N f > N d points. Then we perform a local search around the best grid point τ  grid . In particular, it is shown in Appendix A that (47) has the follo wing equiv alent form τ  = arg min τ ∈ [0 , 1) T r { A ⊥ ( τ ) R res } , (48) where R res = ( Φ − 1 r res )( Φ − 1 r res ) H ∈ C N d × N d , (49) A ⊥ ( τ ) = I N d − 1 N d a ( τ ) a ( τ ) H ∈ C N d × N d . (50) Problem (48) can be solved using Newton’ s method as τ i +1 = τ i − µ i K ( τ i ) − 1 p ( τ i ) , i = 0 , 1 , ... (51) where τ 0 = τ  grid ; µ i is the step size which is chosen according to the backtracking line search [51], gi ven in Appendix B; p ( τ ) and K ( τ ) are the gradient and Hessian, giv en respecti vely by [52], [53] p ( τ ) = ∇ τ  T r { A ⊥ ( τ ) R res }  = − 2Re  1 N d a H ( τ ) R res A ⊥ ( τ ) ∂ a ( τ ) ∂ τ  ∈ R , (52) K ( τ ) = ∇ 2 τ  T r { A ⊥ ( τ ) R res }  ≈ 2Re  ( ∂ a ( τ ) ∂ τ ) H A ⊥ ( τ ) ∂ a ( τ ) ∂ τ  a ( τ ) H R res a ( τ ) N 2 d  ∈ R , (53) 14 where ∂ a ( τ ) ∂ τ =  1 , i 2 π e i 2 π τ , ..., i 2 π ( N d − 1) e i 2 π ( N d − 1) τ  T ∈ C N d × 1 . (54) The iteration in (51) stops when | K ( τ i ) − 1 p ( τ i ) | < δ , where δ is the error tolerance, or the maximum iteration number I is reached. S-1(c) Updating τ using Newton’ s method: If T is updated in step S-1(b), for the current ˆ v , using ¯ z = z − H ˆ v , we refine the estimates of the delays in T by solving the follo wing problem min τ ∈ [0 , 1) |T | , c ∈ C |T | k ¯ z − ΦΘ ( τ ) c k 2 2 . (55) Substituting the solution ˆ c = Θ ( τ ) † Φ − 1 ¯ z , where ( · ) † denotes the pseudo-in v erse, i.e., Y † = ( Y H Y ) − 1 Y H , back to (55) results in the follo wing optimization problem [52], [53]: ˆ τ = arg min τ ∈ [0 , 1) |T | T r { P ⊥ ( τ ) R } , (56) where R = ( Φ − 1 ¯ z )( Φ − 1 ¯ z ) H ∈ C N d × N d , (57) P ⊥ ( τ ) = I N d − Θ ( τ ) Θ ( τ ) † ∈ C N d × N d . (58) Problem (56) can be solved using Newton’ s method as τ i +1 = τ i − ¯ µ i K ( τ i ) − 1 p ( τ i ) , i = 0 , 1 , ... (59) where τ 0 is taken as the current elements in T ; ¯ µ i is the step size which is chosen according to the backtracking line search [51], giv en in Appendix B; p ( τ ) and K ( τ ) are the gradient and Hessian matrix, gi ven respectiv ely by [52], [53] p ( τ ) = ∇ τ  T r { P ⊥ ( τ ) R }  = − 2Re  vec-diag  Θ † ( τ ) R P ⊥ ( τ ) T ( τ )  ∈ R |T |× 1 , (60) K ( τ ) = ∇ 2 τ  T r { P ⊥ ( τ ) R }  ∈ R |T |×|T | ≈ 2Re  ( T ( τ ) H P ⊥ ( τ ) T ( τ ))  ( Θ ( τ ) † R Θ ( τ ) † H ) T  , (61) 15 where vec-diag [ Y ] , with Y being a square matrix, denotes a column vector formed by the diagonal elements of Y , and T ( τ ) is gi ven by T ( τ ) = " ∂ a ( τ ) ∂ τ     τ = τ 1 , ∂ a ( τ ) ∂ τ     τ = τ 2 , ..., ∂ a ( τ ) ∂ τ     τ = τ |T | # ∈ C N d ×|T | =        1 1 · · · 1 i 2 π e i 2 π τ 1 i 2 π e i 2 π τ 2 · · · i 2 π e i 2 π τ |T | . . . . . . . . . . . . i 2 π ( N d − 1) e i 2 π ( N d − 1) τ 1 i 2 π ( N d − 1) e i 2 π ( N d − 1) τ 2 · · · i 2 π ( N d − 1) e i 2 π ( N d − 1) τ |T |        . (62) The iteration in (59) stops when k K ( τ i ) − 1 p ( τ i ) k 2 < δ or the maximum iteration number I is reached. S-1(d) Updating ( v , c ) using least-squar es: W ith the current S and T , estimate v and c by solving the following least-squares problem: ( ˆ v , ˆ c ) = arg min v ( S ) ∈ C |S | c ∈ C |T | k z − H (: , S ) v ( S ) − ΦΘ ( ˆ τ ) c k 2 2 , (63) where H (: , S ) and v ( S ) denote the columns and elements of H and v respectiv ely index ed by S . Then, we remov e any atoms in T whose corresponding coef ficients are smaller than a small threshold ˜ δ . S-1(e) Residual update: r res = z − H ˆ v − ΦΘ ( ˆ τ ) ˆ c (64) and repeat steps S-1(b) to S-1(e) until k r res k 2 2 ≤  , or the maximum iteration number I 0 is reached. 2) Conjugate gradient descent for solving S-2: The constraint k g k 2 = 1 in (44) can be regarded as forcing g on a unit sphere, which belongs to the Riemannian manifolds. W e thus resort to the conjugate gradient method on Riemannian manifold [54] to update g . The Euclidean gradient of the cost function in (44) is q ( g ) = ∇ g k ¯ z − W g k 2 2 = − 2 W H ( ¯ z − W g ) . (65) Projecting the Euclidean gradient to the tangent space of Riemannian manifold yields the Rie- mannian gradient [54], [55] q R ( g ) = q ( g ) − Re( q ( g )  g ∗ )  g . (66) 16 Then the search direction in the i -th iteration is gi ven by [54] q C ( g i ) = γ i [ q C ( g i − 1 ) − Re( q C ( g i − 1 )  ( g i ) ∗ )  g i ] − q R ( g i ) , i = 1 , 2 ..., (67) with q C ( g 0 ) = − q R ( g 0 ) , where γ i = q R ( g i ) T ( q R ( g i ) − q R ( g i − 1 )) k q R ( g i − 1 ) k 2 2 . (68) Then, problem (44) can be solved by the follo wing iterations g i +1 = g i + ˜ µ i q C ( g i ) k g i + ˜ µ i q C ( g i ) k 2 2 , i = 0 , 1 , ..., (69) where ˜ µ i is a step size, which is also chosen according to the backtracking line search [51], gi ven in Appendix B. The iteration in (69) stops when k g i +1 − g i k 2 < ¯ δ , where ¯ δ is the error tolerance, or the maximum iteration number ¯ I is reached. B. Stage 2 - Inferring the Global Optimum After Stage 1, we obtain a locally optimum solution ( ˆ g , ˆ τ , ˆ c ) to problem (26). T o further search for the global optimum, we make use of a theoretical result in [35]. Recall that ( F L g )  ν τ = F ( g ~ ( F − 1 ν τ )) . When g ∈ C L × 1 with L  N d , and F − 1 ν τ ∈ C N d is a sparse vector , then g ~ ( F − 1 ν τ ) is the so-called short-and-sparse (SaS) con volution. It is shown in [35] that, for the SaS blind decon v olution problem, if F − 1 ν τ follo ws the Bernoulli-Gaussian (BG) model, then any local optimum ˆ g is close to certain signed shift truncation of the global optimum g ? with high probability . The signed shift truncation is the result of truncation, circular shift and sign changes on a sequence (see the two examples in Fig. 2). Hence, we speculate that the estimated code ˆ g obtained by Stage 1 is close to a signed shift truncation of the global optimum g ? . In fact, this conjecture is corroborated by e xtensi ve simulations. For e xample, the landscape of the objecti ve function in (26) when v = 0 and L = 3 is shown in Fig. 3. In particular , for a gi ven point on the sphere k g k 2 = 1 , we calculate the corresponding min k ν τ k A , 0 via steps S-1(a)-(e). Dark blue represents small values while dark red represents large values. The landscape clearly sho ws that (26) is non-con ve x. Furthermore, we calculate all the signed shift truncations of the ground truth g ? and mark them on the sphere. W e can see that the local optima are very close to certain signed shift truncations of the ground truth. As the local optimum still captures a considerable portion of the global optimum, then in a higher dimensional space, the zero-padded local optimum should be close to one cyclic shift 17 0 1 2 3 4 G l o b a l o p t i m u m 0 1 2 3 4 0 1 2 3 4 S i g n e d s h i f t t r u n c a t i o n Fig. 2. T wo examples of signed shift truncation. Fig. 3. Geometry of the objective function of (26) on the ` 2 ball when error v = 0 . Dark blue represents small values indicating a local optimum. All local optima are close to signed shift truncations of the ground truth g ? . of the zero-padded global optimum (see Fig. 4). Hence, we first estimate a cyclic shifted zero- padded global optimum instead of estimating the global optimum directly , and the zero-padded local optimum serv es as a significantly better initialization than a random initial value in a higher dimensional space [34]. The estimated ˆ g of Stage 1 after zero-padding is ˜ g 0 =  ¯ 0 T , ˆ g T , ¯ 0 T  T ∈ C (3 L − 2) × 1 , where ¯ 0 ∈ C ( L − 1) × 1 is the all-zero vector 2 . Problem (26) in a higher dimensional space is gi ven by ( ˜ g h , ˜ ν , ˜ v ) = arg min ν ∈ C N d × 1 , v ∈ C N d × 1 g h ∈ C (3 L − 2) × 1 k ν k A , 0 + λ k v k 0 , (70) s.t.   z − H v − ¯ ξ  ( F 3 L − 2 g h )  ν   2 2 ≤ , k g h k 2 = 1 , where F 3 L − 2 denotes the first 3 L − 2 columns of F . Solving (70) by using the same alternating minimization algorithm outlined in Section IV .A with initial v alue ˜ g 0 yields the estimates ˜ g h , ˜ c , 2 Note that if ˆ g is shifted by more than its own length L , there will be no truncation of global optimum g ? retained. Therefore, the length of zero padding vector ¯ 0 is set as L − 1 . 18 0 1 2 3 4 G l o b a l o p t i m u m L o c a l o p t i m u m Z e r o - p a d d e d l o c a l o p t i m u m 0 1 2 3 4 0 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 Z e r o - p a d d e d g l o b a l o p t i m u m C y c l i c s h i f t S o m e p a r t c l o s e H i g h e r d i m e n s i o n a l 0 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 Fig. 4. An example of local optimum and global optimum, and their relationship in the higher dimensional space. Some part of the zero-padded local optimum and the cyclic shift of the zero-padded global optimum are close. 0 1 2 3 4 5 6 7 D e l a y 0 . 5 0 1 0 . 5 1 = 0 1 2 3 4 D e l a y 0 . 5 0 1 = G l o b a l o p t i m u m Z e r o - p a d d e d g l o b a l o p t i m u m Fig. 5. Shift ambiguity in short-and-sparse deconv olution. ˜ τ and ˜ v . Since ˜ g h is an estimate of the cyclic shifted zero-padded global optimum g ? , we need to extract the estimate of g ? by detecting the first element that is larger than a small threshold δ h , i.e., ˜ ` = arg min ` ∈{ 0 , 1 ,..., 3 L − 2 } `, s.t. ˜ g h ( ` ) > δ h , (71) where ˜ g h ( ` ) is the ` -th element of ˜ g h . Then the estimated global optimum code is ˜ g ? = [ ˜ g h ( ˜ ` ) , ˜ g h ( ˜ ` + 1) , ..., ˜ g h ( ˜ ` + L − 1)] T ∈ C L × 1 . Note that con volution has the shift ambiguity property . Extracting g ? from ˜ g h is equi valent to a cyclic shift of length − ˜ ` N d on τ in the con v olution g ~ ( F − 1 ν τ ) (see Fig. 5). Hence, the estimated global optimum delay is ˜ τ ? = ˜ τ + ˜ ` N d . Finally , we summarize the proposed two-stage alternating minimization (2-AltMin) method in Algorithm 1. The main computational load of Algorithm 1 is the calculation of gradient and Hessian in Newton’ s method, with complexities O ( N 3 d ) and O ( M r N 2 d ) per iteration, respectiv ely . 19 Algorithm 1 T wo-stage alternating minimization procedure for solving (26). Input z , H , ¯ ξ , ˜ λ ,  , δ , ¯ δ , ˜ δ , I , I 0 , ¯ I and δ h . 1, Initialize ˆ g as a random code. Repeat (Stage 1) 2, Obtain ˆ v , ˆ τ and ˆ c via S-1. 3, Obtain ˆ g via S-2. 4, ˆ b ← ˜ b = arg min b ∈B N d k b − ˆ b − ˆ v k 2 . 5, z = ¯ r − H ˆ b . Until ˜ b = ˆ b . 6, ˜ g 0 =  ¯ 0 T , ˆ g T , ¯ 0 T  T . Repeat (Stage 2) 7, Obtain ˜ g h , ˜ τ , ˜ c and ˜ v by performing steps 2-5 with the higher dimensional g h initialized as ˜ g 0 . Until ˜ b = ˆ b . 8, Obtain ˜ ` by using (71). 9, ˜ g ? = [ ˜ g h ( ˜ ` ) , ˜ g h ( ˜ ` + 1) , ..., ˜ g h ( ˜ ` + L − 1)] T . 10, ˜ τ ? = ˜ τ + ˜ ` N d . Return ˜ g ? , ˜ τ ? , ˜ c ? and ˜ v ? . Hence the computational comple xity of Algorithm 1 is O ( N 3 d ) per iteration. On the other hand, the complexity of the con ve x relation (CR) method discussed in Section III is O (( N d + L ) 6 ) per iteration if the interior point method is used [18]. Hence the proposed 2-AltMin method is both computationally more efficient and more accurate as sho wn by simulation results in the next section. V . S I M U L A T I O N R E S U L T S A. Baseline for Comparison: On-grid Method As a baseline of comparison, we consider the on-grid method for estimating the continuous delays { τ } , by using an ov ercomplete dictionary matrix ˜ A = [ a 0 , a 1 , ..., a ˜ M − 1 ] ∈ C N d × ˜ M , (72) where ˜ M ≥ N d and a m = a ( m ˜ M ) , m = 0 , 1 , ..., ˜ M − 1 . For sufficiently large ˜ M , the delay is densely sampled. Follo wing the con ve x relaxation used in Section III, define ς = [ c 1 g T 1 , c 2 g T 2 , ..., c ˜ M g T ˜ M ] T ∈ C ˜ M L × 1 (73) 20 as the sparse vector whose non-zero elements correspond to c m g in (28). The original problem (26) can be relaxed to the following on-grid optimization problem ( ˆ ς , ˆ v ) = arg min ς ∈ C ˜ M L × 1 , v ∈ C N d × 1 k ς k 1 + ¯ η k v k 1 , (74) s.t. k z − H v − Υ ς k 2 2 ≤ , where ˜ η is a weight factor and Υ is gi ven by Υ =        a H 0 e 0 d H 0 a H 1 e 0 d H 0 · · · a H ˜ M − 1 e 0 d H 0 a H 0 e 1 d H 1 a H 1 e 1 d H 1 · · · a H ˜ M − 1 e 1 d H 1 . . . . . . . . . . . . a H 0 e N d − 1 d H N d − 1 a H 1 e N d − 1 d H N d − 1 · · · a H ˜ M − 1 e N d − 1 d H N d − 1        ∈ C N d × ˜ M L . (75) Since problem (74) is con ve x, it can be solved with standard con ve x solvers, e.g., CVX [46]. And the complexity in each iteration is O (( ˜ M L + N d ) 3 ) if the interior point method is used [18]. Then, the radar delays and code can be identified by locating the non-zero entries of ˆ ς , i.e., if [ ˆ ς mL , ˆ ς mL +1 , ..., ˆ ς ( m +1) L − 1 ] has elements larger than a pre-set small threshold, then a radar delay exists at m ˜ M N d T and normalizing [ ˆ ς mL , ˆ ς mL +1 , ..., ˆ ς ( m +1) L − 1 ] yields the corresponding estimated radar code. Note that this on-grid method is also a relaxed method, and similar to the example gi ven in Section III, it can be sho wn that some columns of Υ can be identical. Hence, Υ is coher ent [20] and many delay false alarms could be generated in ς , which will be illustrated in the simulations. B. Simulation Setup In order to demonstrate the performance of the proposed algorithms, we consider a scenario where a radar transmitter produces multiple reflections to wards a communication receiv er . The communication system uses an OFDM signal with N d = 256 , N p = 64 and a total bandwidth of 2 . 56 MHz, i.e, the frequency spacing between adjacent subcarriers is 10 kHz. Hence the duration of data symbols N d T = 100 µs and a quadrature phase-shift keying (QPSK) modulation is used. The transmitted OFDM signal is generated according to (1) with normalized data symbols. Since the communication takes place o ver a multi-path Rayleigh-fading channel (see eq. (4)), the path gains { α m } are i.i.d. complex Gaussian distributed, α m ∼ C N (0 , σ 2 h ) . Based on (10), we define the SNR at the communication RX as SNR = E {| P M c m =1 α m e − i 2 π k τ c m N d T | 2 } σ 2 w = P M c m =1 E {| α m | 2 } σ 2 w = M c σ 2 h σ 2 w , (76) 21 0 20 40 60 80 100 0 0.5 magnitude radar interference 0 20 40 60 80 100 0 0.5 magnitude communication signal 0 20 40 60 80 100 0 0.5 magnitude radar interference + communication signal (a) 0 50 100 150 200 250 -5 0 5 real part H -1 radar interference 0 50 100 150 200 250 -1 0 1 real part communication data 0 50 100 150 200 250 sample in frequency -5 0 5 real part H -1 radar interference + communication data (b) Fig. 6. Plots of signal wav eforms in (a) time domain and (b) frequency domain. In (a), the magnitude of the radar interference, communication signal and the receiv ed signal of communication RX are plotted against time. In (b), the real part of the interference on communication data, communication data and their combination are plotted versus frequency sample. where σ 2 w is the variance of the Gaussian noise sample w ( k ) in (10). In the following simulations, we set M c = 10 and σ 2 h = 0 . 1 . In (3), the radar uses a pulse coded wa veform and pulse uses Gaussian random code with length L , and then we normalize the code to let k g k 2 = 1 for the simplicity of ev aluation. The sub-pulse of radar signal ξ ( t ) , t ∈ [0 , T ] is set as the normalized rectangular pulse of duration T . The reference delay τ R and the delays τ r m of the scatters are randomly generated between 0 and 100 µs . The radar PRI is set as N d T = 100 µs . The scatterers are modeled as point sources in our simulations, and the complex scattering coef ficient c m of the m -th scatter is generated with fixed magnitude c 0 and random phase for con venience of ev aluation. Specifically , based on (10), we define the interference-to-signal ratio (ISR), which is the av erage power ratio of the radar interference and the communication signal, at the communication RX as ISR = 1 N d N d − 1 P k =0 E (     ¯ g ( k ) ¯ ξ ( 2 π k N d T ) M r P m =1 c m e i 2 π kτ m     2 ) M c σ 2 h = M r | c 0 | 2 N d − 1 P k =0    ¯ g ( k ) ¯ ξ ( 2 π k N d T )    2 N d M c σ 2 h . (77) W e ev aluate the mean absolute error (MAE) of the radar delay estimate and the mean-squared- error (MSE) 3 of the estimated radar code for the on-grid method, the CR method and the 2-AltMin 3 W e use the relative MSE rather than the RMSE to ev aluate the accuracy because it reflects the loss in energy . 22 algorithm. Note that in each case the algorithm returns a b unch of delays, which can be either true detections or false alarms. In the simulations, for each estimated delay ˆ τ r ` , ` = 1 , ..., |T | , we calculate the minimum absolute error AE ` with the ground truth delays τ r m , m = 1 , ..., M r , i.e., AE ` = min( { ˆ τ r ` − τ r m } M r m =1 ) . Then, the delay MAE and the relativ e wav eform MSE are respecti vely calculated as MAE τ = 1 MC MC X n MC =1 1 |T | |T | X ` =1 AE ( n MC ) ` , (78) MSE g = 1 MC MC X n MC =1    | g ( n MC ) | − | ˜ g ( n MC ) ? |    2 2 k| g ( n MC ) |k 2 2 , (79) where MC is the number of Monte Carlo runs; AE ( n MC ) ` is the minimum absolute error of the ` -th estimate in the n MC -th run; g ( n MC ) and ˜ g ( n MC ) ? are the radar code and the estimated radar code at the n MC -th run, respectiv ely . The error tolerance is usually set smaller than  w 0 . 05 k z k 2 2 , which implies that the iteration stops when the relativ e error is smaller than 5% [56]. For the proposed algorithms, we set the error tolerance in (26) as  w 0 . 01 k z k 2 2 for better performance. The weight factors for the on-grid method in (74) and the CR method in (35) are respectiv ely set as ¯ η = 1 and ¯ λ = 1 √ N d [18]. And the weight factor for the 2-AltMin algorithm in step S-1(b) is set as ˜ λ = 6 √ N d M c σ 2 h . The grid number ˜ M in (72) is set as ˜ M = 512 . The error tolerances for Ne wton’ s method and conjugate gradient method are both set as δ = ¯ δ = 10 − 6 , and the threshold in step S-1(d) and (71) are respecti vely set as ˜ δ = 0 . 05 and δ h = 0 . 05 . The maximum iteration numbers for Ne wton’ s method, step S-1(e) and the conjugate gradient method are respectiv ely set as I = 10 , I 0 = 50 and ¯ I = 10 . In addition, ρ and ¯ ρ for the backtracking line search in Algorithm 2 are respectiv ely set as 0 . 5 and 0 . 01 . In order to show the performance of the proposed methods, we compare the symbol error rate (SER) of the proposed methods with the SER of directly performing demodulation using ¯ r , which is named “Iteration 0” because its result is the initial point of the iterati ve algorithms. In addition, we compare the performance of the Stage 1 of the 2-AltMin method, which is named “Stage 1”. 23 0 5 10 15 SNR [dB] 10 -2 10 -1 10 0 SER Iteration 0 On-grid CR Stage 1 2-AltMin (a) 0 5 10 15 SNR [dB] 10 -3 10 -2 10 -1 SER Iteration 0 On-grid CR Stage 1 2-AltMin (b) Fig. 7. SER performance comparison when the ISR of communication is (a) 5 dB and (b) -5 dB. (a) 0 5 10 15 20 25 n MC 0 1 2 3 4 5 Stage 1 2-AltMin (b) Fig. 8. Realizations of the delay minimum absolute errors. (a) On-grid and CR methods, (b) Stage 1 and 2-AltMin methods. C. P erformance In the first simulation, the number of scatterers is set as M r = 2 and the length of radar pulse is set at L = 10 . Fig. 6 gi ves the signal at the communication RX and the interference and data for demodulation when the ISR is set at − 5 dB. W e can find that the effect of interference is significant e ven if there are only two multi-paths radar echo and the ISR is moderate. Then, we compare the SER performance of v arious algorithms. In Fig. 7, the ef fect of the SNR is studied: 24 the on-grid, CR and 2-AltMin methods all provide better SER performance than Iteration 0. The 2-AltMin method also outperforms the on-grid and CR methods in all situations. It is worth mentioning that there is a significant improv ement when the Stage 2 of 2-AltMin is used in all cases. W e then ev aluate the relativ e code MSE and delay MAE of the proposed methods. W e first plot the delay minimum absolute errors AE ( n MC ) ` of different methods when the SNR is 15 dB in Fig. 8. W e can clearly see that the on-grid, CR and Stage 1 methods all reach many local optima, while the 2-AltMin method reaches the global optimum with high probability . And the on-grid method produces a large number of delay false alarms, because some columns of Υ in (75) are coher ent . Then, the relativ e wa v eform MSE and delay MAE are plotted against the SNR in Fig. 9. Note that when the SNR is low , there are some very large delay minimum absolute errors, which af fect the analysis of the av erage. Hence we remov e the minimum absolute errors that are larger than 5 µs . Then the delay MAE in Fig. 9(b) is calculated according to (78). As expected, the 2-AltMin method provides much better accuracy than other methods in all situations. In addition, we can see that the interference estimation accuracy may not necessarily improv e with the ISR. When the SNR is lo w , the delay estimation performance is better when ISR = 5 dB, because strong radar interference can pre v ail the noise. While when the SNR is large, the delay estimation performance is better when ISR = − 5 dB because large SNRs guarantee good demodulation performance, with a beneficial effect on the radar interference estimation due to the coupling. The effects of M r and L are shown in Fig. 10. The simulations are run with an SNR of 15 dB and an ISR of 5 dB. In Fig. 10(a), we set L = 6 and plot the SER against the number of scatterers: As M r increases, the sparsity of the problem is reduced, and the sources of interference - with the respectiv e unknown parameters to be estimated - increase, which obviously results in a visible performance degradation for all algorithms. In Fig. 10(b), the number of scatterers is set as M r = 2 and we examine the SER beha vior for v arying radar pluse length L . A performance degradation is also evident for all algorithms. Finally , we giv e an example of the con ver gence behavior of the three algorithms, which is sho wn in Fig. 11. The number of scatterers, the length of radar pulse, the ISR and the SNR are respecti vely set as M r = 2 , L = 10 , ISR = 5 dB and SNR = 15 dB. The on-grid method takes 1616.3 seconds with 4 iterations by using CVX [46]. The CR method takes 1203.5 seconds with 5 iterations by using CVX, while the 2-AltMin method only takes 2.3 seconds with 12 iterations 25 0 2 4 6 8 10 12 14 16 18 SNR[dB] 10 -1 10 0 Relative Code MSE On-grid,ISR=5dB CR,ISR=5dB Stage 1,ISR=5dB 2-AltMin,ISR=5dB On-grid,ISR=-5dB CR,ISR=-5dB Stage 1,ISR=-5dB 2-AltMin,ISR=-5dB (a) 0 2 4 6 8 10 12 14 16 18 SNR[dB] 0 0.5 1 1.5 2 2.5 3 3.5 On-grid,ISR=5dB CR,ISR=5dB Stage 1,ISR=5dB 2-AltMin,ISR=5dB On-grid,ISR=-5dB CR,ISR=-5dB Stage 1,ISR=-5dB 2-AltMin,ISR=-5dB (b) Fig. 9. (a) Relative wav eform MSE performance comparisons. (b) Delay MAE performance comparisons. 2 4 6 8 10 12 14 16 18 20 M r 10 -3 10 -2 10 -1 10 0 SER Iteration 0 On-grid CR Stage 1 2-AltMin (a) 2 4 6 8 10 12 14 16 18 20 L 10 -4 10 -3 10 -2 10 -1 10 0 SER Iteration 0 On-grid CR Stage 1 2-AltMin (b) Fig. 10. SER performance against (a) M r , and (b) L . total (8 iterations in Stage 1 and 4 iterations in Stage 2). The experiments were carried out on a MacBook Pro computer with a 2.3 GHz Intel Core i5 CPU and 8 GB of RAM. The proposed 2-AltMin method is substantially faster than the CR method and the on-grid method and appears much well suited for real-time implementations. 26 0 2 4 6 8 10 12 Iteration number 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 SER On-grid CR 2-AltMin Fig. 11. Conv ergence behavior of dif ferent algorithms. V I . C O N C L U S I O N S In this paper , we have proposed two algorithms for removing the radar interference to facilitate more reliable data demodulation in a communication system ov erlaid with a radar system. The first one is based on forcing an atomic norm constraint, and estimating the combination of the radar parameters by solving a con vex problem under some relaxations. The second algorithm estimates the radar parameters and the communication demodulation errors by two- stage processing. The first stage obtains a local optimum by alternating minimization, and the second stage infers the global optimum in a higher dimensional space by using the estimates of the first stage. The atomic norm and the ` 0 -norm are used to exploit the sparsity of the radar signal components and the sparsity of the demodulation error , respectiv ely . Simulation results sho w that both algorithms provide much better SER performance compared to the con ventional on-grid method. Moreover , the proposed 2-AltMin algorithm of fers superior performance and is computationally efficient. 27 A P P E N D I X A. Pr oof of (48) By noting that ( Φ − 1 ) H ( Φ − 1 )  0 we ha ve arg max τ ∈ [0 , 1) |h Φ a ( τ ) , r res i| = arg max τ ∈ [0 , 1) | r H res Φ a ( τ ) | = arg max τ ∈ [0 , 1) | r H res ( Φ − 1 ) H ( Φ − 1 ) Φ a ( τ ) | = arg max τ ∈ [0 , 1) | ( Φ − 1 r res ) H a ( τ ) | = arg max τ ∈ [0 , 1) (( Φ − 1 r res ) H a ( τ )) H ( Φ − 1 r res ) H a ( τ ) = arg max τ ∈ [0 , 1) T r { a ( τ ) a ( τ ) H ( Φ − 1 r res )( Φ − 1 r res ) H } = arg min τ ∈ [0 , 1) T r { ( I N d − a ( τ ) a ( τ ) H N d ) | {z } A ⊥ ( τ ) ( Φ − 1 r res )( Φ − 1 r res ) H | {z } R res } . (80) Since a ( τ ) = [1 , e i 2 π τ , ..., e i 2 π ( N d − 1) τ ] T , then 1 N d a ( τ ) H = ( a ( τ ) H a ( τ )) − 1 a ( τ ) H = a ( τ ) † . Thus we have A ⊥ ( τ ) = I N d − a ( τ ) a ( τ ) † in the last line. B. Backtrac king Line Sear ch The backtracking line search approach ensures the selected step size is small enough to guar - antee a sufficient decrease of the cost function b ut not too small. For simplify , define the objectiv e functions for (51), (59) and (69) respecti vely as L ( τ ) = Tr { A ⊥ ( τ ) R res } , L ( τ ) = T r { P ⊥ ( τ ) R } and L ( g ) = k ¯ z − W g k 2 2 . And define their search directions respectiv ely as D ( τ ) = K ( τ ) − 1 p ( τ ) , D ( τ ) = K ( τ ) − 1 p ( τ ) and D ( g ) = − q C ( g ) . As an example, in Algorithm 2 we summarize the backtracking line search for calculating µ i in (51). Then ¯ µ i in (59) and ˜ µ i in (69) can Algorithm 2 Backtracking line search Input L ( τ ) , τ i , D ( τ ) , ρ ∈ (0 , 1) and ¯ ρ ∈ (0 , 1 / 2) . 1, Initialize µ i = 1 . 2, Repeat 3, µ i = ρµ i 4, Until L ( τ i − µ i D ( τ i )) ≤ L ( τ i ) − ¯ ρµ i kD ( τ i ) k 2 2 . Return µ i . 28 be obtained with Algorithm 2 by replacing ( L ( τ ) , D ( τ ) , τ i , µ i ) with ( L ( τ ) , D ( τ ) , τ i , ¯ µ i ) and ( L ( g ) , D ( g ) , g i , ˜ µ i ) , respectively . R E F E R E N C E S [1] H. Griffiths, L. Cohen, S. W atts, E. Mokole, C. Baker , M. W icks, and S. Blunt, “Radar spectrum engineering and management: technical and regulatory issues, ” Proc. IEEE , vol. 103, no. 1, pp. 85–102, Jan. 2015. [2] A. Hassanien, M. G. Amin, Y . D. Zhang, and F . Ahmad, “Dual-Function Radar-Communications: Information Embedding Using Sidelobe Control and W aveform Div ersity .” IEEE T rans. Sig. Pr oc , vol. 64, no. 8, pp. 2168–2181, 2016. [3] P . Kumari, N. Gonzalez-Prelcic, and R. W . Heath, “Inv estigating the ieee 802.11ad standard for millimeter wav e automotive radar , ” in V ehicular T echnology Conference (VTC F all), 2015 IEEE 82nd . IEEE, 2015, pp. 1–5. [4] E. Grossi, M. Lops, L. V enturino, and A. Zappone, “Opportunistic radar in 802.11ad networks, ” IEEE T ransactions on Signal Pr ocessing , vol. 66, no. 9, pp. 2441–2454, 2018. [5] A. R. Chiriyath, B. Paul, G. M. Jacyna, and D. W . Bliss, “Inner bounds on performance of radar and communications co-existence, ” IEEE T rans. Signal Pr ocess. , vol. 64, no. 2, pp. 464–474, Jan. 2016. [6] F . Hessar and S. Roy , “Spectrum sharing between a surveillance radar and secondary W i-Fi networks, ” IEEE T rans. Aer osp. Electr on. Syst. , v ol. 52, no. 3, pp. 1434–1448, Jun. 2016. [7] Z. Ding, B. Shu, W . Y in, T . Zeng, and T . Long, “ A modified frequency domain algorithm based on optimal azimuth quadratic factor compensation for geosynchronous SAR imaging, ” IEEE J. Sel. T opics Appl. Earth Observ . , vol. 9, no. 3, pp. 1119–1131, 2016. [8] A. Babaei, W . H. T ranter , and T . Bose, “ A practical precoding approach for radar/communications spectrum sharing, ” in Pr oc. 8th Int. Conf. Cognitive Radio Oriented Wir eless Netw . , 2013, pp. 13–18. [9] S. Sodagari, A. Khawar , T . C. Clancy , and R. McGwier, “ A projection based approach for radar and telecommunication systems coexistence, ” in Pr oc. Global Commun. Conf. , 2012, pp. 5010–5014. [10] H. Deng and B. Himed, “Interference mitigation processing for spectrum-sharing between radar and wireless communica- tions systems, ” IEEE T rans. Aerosp. Electr on. Syst. , vol. 49, no. 3, pp. 1911–1919, 2013. [11] A. Aubry , A. De Maio, M. Piezzo, and A. Farina, “Radar waveform design in a spectrally crowded en vironment via noncon ve x quadratic optimization, ” IEEE T rans. Aerosp. Electr on. Syst. , vol. 50, no. 2, pp. 1138–1152, 2014. [12] K.-W . Huang, M. Bic ˘ a, U. Mitra, and V . K oivunen, “Radar wav eform design in spectrum sharing environment: Coexistence and cognition, ” in Proc. Radar Conf. , 2015, pp. 1698–1703. [13] B. Li, A. Petropulu, and W . Trappe, “Optimum co-design for spectrum sharing between matrix completion based MIMO radars and a MIMO communication system, ” IEEE T rans. Signal Process. , vol. 64, no. 17, pp. 4562–4575, 2016. [14] L. Zheng, M. Lops, X. W ang, and E. Grossi, “Joint design of overlaid communication systems and pulsed radars, ” IEEE T r ans. Signal Pr ocess. , vol. 66, no. 1, pp. 139–154, 2018. [15] A. T urlapaty and Y . Jin, “ A joint design of transmit wav eforms for radar and communications systems in coexistence, ” in Pr oc. Radar Conf . , 2014, pp. 0315–0319. [16] A. Khawar , A. Abdel-Hadi, and T . C. Clanc y , “Spectrum sharing between S-band radar and L TE cellular system: A spatial approach, ” in Proc. Int. Symp. Dyn. Spectrum Access Netw . , 2014, pp. 7–14. [17] A. Manolakos, Y . Noam, K. Dimou, and A. J. Goldsmith, “Blind null-space tracking for MIMO underlay cognitiv e radio networks, ” in Pr oc. Global Commun. Conf. , 2012, pp. 1223–1229. [18] L. Zheng, M. Lops, and X. W ang, “ Adaptive interference removal for uncoordinated radar/communication coexistence, ” IEEE J. Sel. T op. Signal Pr oces. , vol. 12, no. 1, pp. 45–60, 2018. 29 [19] J. Liu, H. Li, and B. Himed, “Joint optimization of transmit and receive beamforming in acti ve arrays, ” IEEE Signal Pr ocess. Lett. , v ol. 21, no. 1, pp. 39–42, Jan. 2014. [20] E. J. Candes, Y . C. Eldar , D. Needell, and P . Randall, “Compressed sensing with coherent and redundant dictionaries, ” Appl. Comput. Harmon. Anal. , vol. 31, no. 1, pp. 59–73, 2011. [21] X. Zhang, W . Cui, and Y . Liu, “Recovery of structured signals with prior information via maximizing correlation, ” IEEE T r ans. Signal Pr ocess. , vol. 66, no. 12, pp. 3296–3310, 2018. [22] Y . Chi, L. L. Scharf, A. Pezeshki, and A. R. Calderbank, “Sensitivity to basis mismatch in compressed sensing, ” IEEE T r ans. Signal Pr ocess. , vol. 59, no. 5, pp. 2182–2195, 2011. [23] L. Stankovi ´ c, I. Orovi ´ c, S. Stankovi ´ c, and M. Amin, “Compressive sensing based separation of nonstationary and stationary signals overlapping in time-frequency , ” IEEE T rans. Signal Process. , vol. 61, no. 18, pp. 4562–4572, 2013. [24] B. Jokanovic and M. Amin, “Reduced interference sparse time-frequency distributions for compressed observations, ” IEEE T r ans. Signal Pr ocess. , vol. 63, no. 24, pp. 6698–6709, 2015. [25] C. Studer , P . Kuppinger , G. Pope, and H. Bolcskei, “Recov ery of sparsely corrupted signals, ” IEEE T rans. Inf. Theory , vol. 58, no. 5, pp. 3115–3130, 2012. [26] E. J. Cand ` es and C. Fernandez-Granda, “Super -resolution from noisy data, ” J. F ourier Anal. Appl. , vol. 19, no. 6, pp. 1229–1254, 2013. [27] E. J. Cand ` es and C. Fernandez-Granda, “T o wards a mathematical theory of super-resolution, ” Commun. Pur e Appl. Math. , vol. 67, no. 6, pp. 906–956, 2014. [28] G. T ang, B. N. Bhaskar, P . Shah, and B. Recht, “Compressed sensing off the grid, ” IEEE T rans. Inf. Theory , vol. 59, no. 11, pp. 7465–7490, Nov . 2013. [29] B. N. Bhaskar , G. T ang, and B. Recht, “ Atomic norm denoising with applications to line spectral estimation, ” IEEE T rans. Signal Pr ocess. , vol. 61, no. 23, pp. 5987–5999, Dec. 2013. [30] Z. T an, Y . C. Eldar , and A. Nehorai, “Direction of arriv al estimation using co-prime arrays: A super resolution vie wpoint, ” IEEE T rans. Signal Pr ocess. , vol. 62, no. 21, pp. 5565–5576, Nov . 2014. [31] S. Ling and T . Strohmer , “Blind decon volution meets blind demixing: Algorithms and performance bounds, ” IEEE T rans. Inf. Theory , vol. 63, no. 7, pp. 4497–4520, 2017. [32] A. Ahmed, B. Recht, and J. Romberg, “Blind decon v olution using con ve x programming, ” IEEE T rans. Inf. Theory , vol. 60, no. 3, pp. 1711–1732, 2014. [33] Y . Chi, “Guaranteed blind sparse spikes decon volution via lifting and con ve x optimization, ” IEEE J . Sel. T op. Signal Pr oces. , vol. 10, no. 4, pp. 782–794, 2016. [34] Y . Zhang, Y . Lau, H.-w . Kuo, S. Cheung, A. P asupathy , and J. Wright, “On the global geometry of sphere-constrained sparse blind decon volution, ” in Pr oc. IEEE Conf. Comput. V ision P attern Reco gni. , 2017, pp. 4894–4902. [35] Y . Zhang, H.-w . Kuo, and J. Wright, “Structured local minima in sparse blind deconv olution, ” in Adv . Neural Inf. Process. Syst. , 2018, pp. 2324–2333. [36] F . Engels, P . Heidenreich, A. M. Zoubir , F . K. Jondral, and M. Wintermantel, “ Adv ances in automotive radar: A framew ork on computationally efficient high-resolution frequency estimation, ” IEEE Signal Proc. Mag. , vol. 34, no. 2, pp. 36–46, 2017. [37] A. Doufe xi, S. Armour , M. Butler , A. Nix, D. Bull, J. McGeehan, and P . Karlsson, “ A comparison of the hiperlan/2 and ieee 802.11 a wireless lan standards, ” IEEE Commun. mag. , vol. 40, no. 5, pp. 172–180, 2002. [38] W . Zhiguo, L. Xi, and F . Y uanchun, “Moving target position with through-wall radar , ” in Int. Conf. Radar . IEEE, 2006, pp. 1–4. 30 [39] G. W ei, Y . Zhou, and S. W u, “Detection and localization of high speed moving targets using a short-range uwb impulse radar , ” in IEEE Radar Conf. IEEE, 2008, pp. 1–4. [40] C. R. Berger , B. Demissie, J. Heckenbach, P . W illett, and S. Zhou, “Signal processing for passi ve radar using OFDM wa veforms, ” IEEE J. Sel. T op. Signal Proces. , vol. 4, no. 1, pp. 226–238, 2010. [41] A. F . Molisch, W ir eless Communications . John W iley & Sons, 2012. [42] D. Hu, L. He, and X. W ang, “ An efficient pilot design method for OFDM-based cognitiv e radio systems, ” IEEE T rans. W ir eless Commun. , vol. 10, no. 4, pp. 1252–1259, 2011. [43] D. Y ang, G. T ang, and M. B. W akin, “Super-resolution of complex exponentials from modulations with unkno wn wa veforms, ” IEEE T rans. Inf. Theory , vol. 62, no. 10, pp. 5809–5830, 2016. [44] Z. Y ang and L. Xie, “Enhancing sparsity and resolution via reweighted atomic norm minimization, ” IEEE T rans. Signal Pr ocess. , vol. 64, no. 4, pp. 995–1006, 2016. [45] S. Li, D. Y ang, G. T ang, and M. W akin, “ Atomic norm minimization for modal analysis from random and compressed samples, ” IEEE T rans. Signal Process. , 2018. [46] S. Boyd and L. V andenberghe, Con vex Optimization . Cambridge Univ ersity Press, 2004. [47] A. Naha, A. K. Samanta, A. Routray , and A. K. Deb, “Determining autocorrelation matrix size and sampling frequency for MUSIC algorithm, ” IEEE Signal Process. Lett. , vol. 22, no. 8, pp. 1016–1020, 2015. [48] J. A. T ropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit, ” IEEE T rans. Inf. Theory , vol. 53, no. 12, pp. 4655–4666, 2007. [49] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P . E. Wright, “Con ver gence properties of the Nelder-Mead simplex method in low dimensions, ” SIAM J ournal on Optimization , vol. 9, no. 1, pp. 112–147, 1998. [50] C. Fernandez-Granda, G. T ang, X. W ang, and L. Zheng, “Demixing sines and spikes: Robust spectral super-resolution in the presence of outliers, ” Information and Infer ence: A Journal of the IMA , vol. 7, no. 1, pp. 105–168, 2017. [51] D. P . Bertsekas, Nonlinear Pr ogramming . Athena Scientific Belmont, 1999. [52] M. V iber g, B. Ottersten, and T . Kailath, “Detection and estimation in sensor arrays using weighted subspace fitting, ” IEEE T r ans. Signal Pr ocess. , vol. 39, no. 11, pp. 2436–2449, 1991. [53] M. V iberg and B. Ottersten, “Sensor array processing based on subspace fitting, ” IEEE T rans. Signal Process. , v ol. 39, no. 5, pp. 1110–1121, 1991. [54] J.-C. Chen, “Low-P APR precoding design for massi ve multiuser MIMO systems via Riemannian manifold optimization, ” IEEE Commun. Lett. , vol. 21, no. 4, pp. 945–948, 2017. [55] P .-A. Absil, R. Mahony , and R. Sepulchre, Optimization Algorithms on Matrix Manifolds . Princeton Univ ersity Press, 2009. [56] G. Li, H. Zhang, X. W ang, and X.-G. Xia, “ISAR 2-D imaging of uniformly rotating targets via matching pursuit, ” IEEE T r ans. Aer osp. Electr on. Syst. , vol. 48, no. 2, pp. 1838–1846, 2012.