RFI Mitigation for One-bit UWB Radar Systems
Tianyi Zhang, Jiaying Ren, Jian Li, Lam H. Nguyen, Petre Stoica
11 RFI Mitigation for One-bit UWB Radar Systems
Tianyi Zhang, Jiaying Ren, Jian Li,
Fellow, IEEE, and Lam H. Nguyen
Abstract —Radio frequency interference (RFI) mitigation iscritical to the proper operation of ultra-wideband (UWB) radarsystems since RFI can severely degrade the radar imagingcapability and target detection performance. In this paper, weaddress the RFI mitigation problem for one-bit UWB radarsystems. A one-bit UWB system obtains its signed measurementsvia a low-cost and high rate sampling scheme, referred to asthe Continuous Time Binary Value (CTBV) technology. Thissampling strategy compares the signal to a known thresholdvarying with slow-time and therefore can be used to achievea rather high sampling rate and quantization resolution withrather simple and affordable hardware. This paper establishes aproper data model for the RFI sources and proposes a novel RFImitigation method for the one-bit UWB radar system that usesthe CTBV sampling technique. Specifically, we first model theRFI sources as a sum of sinusoids with frequencies fixed duringthe coherent processing interval (CPI) and we exploit the sparsityof the RFI spectrum. We extend a majorization-minimizationbased 1bRELAX algorithm, referred to as 1bMMRELAX, to es-timate the RFI source parameters from the signed measurementsobtained by using the CTBV sampling strategy. We also devise anew fast frequency initialization method based on the AlternatingDirection Method of Multipliers (ADMM) methodology for theextended 1bMMRELAX algorithm to significantly improve itscomputational efficiency. Moreover, an ADMM-based sparsemethod is introduced to recover the desired radar echoes usingthe estimated RFI parameters. Both simulated and experimentalresults are presented to demonstrate that our proposed algorithmoutperforms the existing digital integration method, especially forsevere RFI cases.
Index Terms —Signed measurements, one-bit sampling, time-varying thresholds, one-bit UWB radar, RFI mitigation,majorization-minimization (MM), RELAX, one-bit Bayesian in-formation criterion (1bBIC), sparse recovery, ADMM
I. I
NTRODUCTION U LTRA-WIDEBAND (UWB) radar has been used ina wide range of applications, including, for example,landmine and unexploded ordinance (UXO) detection usingground penetrating radar (GPR) [1], hidden object imagingvia foliage penetrating (FOPEN) radar [2], as well as humandetection [3] and non-contact human vital sign monitoringusing UWB radar [4]. Due to the large bandwidth, whichcan be over 10 GHz for an impulse UWB radar system, ananalog-to-digital converter (ADC) with a high-sampling rateof over 20 GHz is needed at its receiver. However, a radarsystem using such an ADC, especially with high-resolutionquantization, may be too expensive to be commercially viable.
This work was supported in part by the National Science Foundation underGrant 1704240, and in part by the U.S. Army Research Laboratory and theU.S. Army Research Office under Grant W911NF-16-2-0223 (correspondingauthor: Jian Li).T. Zhang, J. Ren, and J. Li are with the department od Electrical andComputer Engineering, University of Florida, Gainesville, FL 32611, USA(e-mail: tianyi.zhang@ufl.edu, jiaying.ren@ufl.edu, [email protected]fl.edu).L.H.Nguyen is with the U.S. Army Research Laboratory, Adelphi, MD20783 USA (e-mail: [email protected])
Indeed, high rate ADC with large quantization depth, evenif available, can significantly increase the cost and powerconsumption of the UWB radar system. In contrast, an ADCwith low-resolution quantization can be attractive due to itslow-cost and low power consumption advantage and its abilityto achieve ultra-high sampling rates [5], [6]. For instance,the NVA6100 impulse radar system [7], a low-cost and lowpower consumption, single-chip UWB radar from Novelda,utilizes the so-called Continuous Time Binary Value (CTBV)technology [7], [8] to achieve a very high sampling rate of39 GHz and a 13-bit quantization resolution with a simplecircuit design. CTBV is an efficient one-bit sampling strategy,which obtains its signed measurements via comparing thereceived signal to a known threshold varying with slow-time,i.e., varying from one pulse repetition interval (PRI) to another.High-precision samples can be obtained from these signedmeasurements via a simple digital integration (DI) method[8]. The affordable NVA6100 system can be used for diverseapplications, including vital sign monitoring [4], through-wallimaging and object tracking [7]. We refer to the CTBV-basedUWB radar system in this paper as the one-bit UWB radarsystem.One of the most significant challenges of ensuring theproper operations of UWB radar systems is to mitigate the se-vere radio frequency interferences (RFIs) they encounter sincethere are many competing users within the ultra-widebandfrequency range they operate in. Typical RFI sources includeFM radio transmitters, TV broadcast transmitters, cellularphones, and other radiation devices. Their operating frequencybands tend to overlap with those of the UWB radar systems[9]. These RFI sources pose a significant hindrance to theproper operations of the UWB radar systems in terms ofreduced signal-to-noise ratio (SNR) and degraded radar imag-ing quality. Therefore, effective RFI mitigation is criticallyimportant for the proper operations of the UWB radar systems.RFI mitigation is a notoriously challenging problem sinceit is difficult to predict and model RFI signals accuratelydue to their dynamic range and diverse modulation schemes.Many RFI mitigation methods, such as RFI suppression viafiltering techniques [9]–[13] and RFI extraction based onRFI estimation methods [14]–[17], have been developed forradar systems using high-precision ADCs. However, it appearsthat RFI mitigation for one-bit UWB radar systems has notbeen considered in the literature before and the existinghigh-resolution quantization based methods are not directlyapplicable.In this paper, we introduce a RFI mitigation method forthe one-bit UWB radar systems using the CTBV samplingtechnique, in particular the NVA6100 impulse radar system.Our main contributions can be summarized as follows:1) We present a novel RFI mitigation framework for one-bit a r X i v : . [ ee ss . SP ] F e b UWB radar systems, which obtain their measurements usinga one-bit sampling strategy that varies the known quantizationthreshold with slow-time.2) We establish a proper data model for the RFI sourcesand extend the recently-developed majorization-minimization(MM) based 1bMMRELAX [18] method for sinusoidal pa-rameter estimation, i.e., for single-PRI based signed measure-ments, to deal with multiple-PRI based signed measurements.Specifically, since the desired UWB radar echo signals arerelatively weak with a flat spectrum and the RFI sourcesare typically strong with sparse, narrow spectral peaks in thefast-time frequency domain (see Figure 1, for example), weconsider modeling the RFI signal as a sum of sinusoids. Thesinusoidal frequencies are assumed fixed within the coherentprocessing interval (CPI). The extended 1bMMRELAX algo-rithm [19] can be used to obtain the maximum likelihood(ML) estimates of the parameters of the RFI sources fromthe multiple-PRI based signed measurements.3) To further reduce the computational cost of the ex-tended 1bMMRELAX, we also devise an Alternating Direc-tion Method of Multipliers (ADMM) [20] based fast frequencyinitialization method which exploits the sparsity of the RFIspectrum.4) Since the number of RFI sources, i.e., the model orderof RFI signals, is unknown, we extend the single-PRI based1bBIC [21] to the multiple-PRI based cases. We use theextended 1bMMRELAX with the extended 1bBIC to simulta-neously estimate the RFI parameters and determine the numberof RFI sources.5) We model the desired UWB radar echoes as sparseimpulses in the fast-time domain due to the sparsity of strongtargets, and introduce an ADMM-based sparse method toefficiently and effectively recover the desired UWB radarechoes based on the estimated RFI parameters.6) Both simulated and measured RFI examples are presentedin this paper to demonstrate the effectiveness of the proposedmethods, especially when the RFI problem is severe.
Slow-time Index F r e qu e n cy ( G H z ) -50-45-40-35-30-25-20-15-10-50(dB) Fig. 1: An example of fast-time RFI spectrum vs. slow-timeindex, for the RFI-only data measured by the experimentalARL radar receiver.The rest of this paper is organized as follows. In SectionII, we introduce the one-bit UWB radar system using the CTBV sampling technique and formulate the RFI mitigationproblem for the system. Next, in Section III, we presentthe extended 1bMMREALX algorithm along with the fastfrequency initialization method and the extended 1bBIC toestimate the RFI parameters and determine the number ofRFI sources from multiple-PRI based signed measurementsobtained by using the CBTV sampling strategy. Then, usingthe estimated RFI parameters, we introduce the ADMM basedsparse echo signal recovery method in Section IV. Finally, inSection V, we provide both simulated and experimental resultsto demonstrate the effectiveness of the proposed algorithm forRFI mitigation for one-bit UWB radar systems.
Notation:
We denote vectors and matrices by boldfacelower-case and upper-case letters, respectively. ( · ) T denotesthe transpose operation. ˆ( · ) refers to the estimated result ofthe related value. X ∈ R N × M denotes a real-valued N × M matrix and x ∈ R N denotes a real-valued vector with N elements. X [ n, m ] means the ( n, m ) th element of matrix X . X [ n, :] and X [: , m ] means the n th row and m th column ofthe matrix X , respectively. x [ n ] denotes the n th element ofthe vector x . For a matrix or a vector, || · || p means the (cid:96) p element-wise norm of this matrix or vector, i.e., || X || p =( (cid:80) Mm =1 (cid:80) Nn =1 | X [ n, m ] | p ) /p or || x || p = ( (cid:80) Nn =1 | x [ n ] | p ) /p . || X || , = (cid:80) Nn =1 || X [ n, :] || denotes the (cid:96) , norm of thematrix X . I N denotes the N × N identity matrix.II. P ROBLEM F ORMULATION
A. One-bit UWB Radar System
Consider an exemplary one-bit UWB radar system, theNVA6100 system, which is a recently-developed UWB radarsystem on a single chip with rather simple circuit designs [7].It can achieve an ultra-high sampling frequency of 39 GHz anda 13-bit quantization resolution with low-cost and low powerconsumption hardware by using the CTBV sampling scheme,a kind of one-bit sampling technique. The radar transmits asuper-narrow pulse repeatedly and uses a one-bit ADC withdifferent thresholds to obtain the signed measurements. Morespecifically, each reflected signal will be compared with aknown threshold, which varies uniformly over different PRIswithin the CPI, and whether the sampling data is larger orsmaller than the threshold will be recorded. Thus, in the ab-sence of RFI and other interferences, the signed measurementmatrix Y ∈ R N × M obtained by NVA6100 can be expressedas follows: Y = sign( S − H ) , (1)where N and M denote the number of fast-time samples perPRI and the number of PRIs or slow-time samples within theCPI, respectively, S denotes the desired radar echo signal and H denotes the known threshold matrix, whose columns varylinearly with slow-time, i.e., H [ n, m ] = − h +2( m − h/ ( M − , h > , n = 1 , , · · · , N, m = 1 , , · · · , M . sign( · ) is theelement-wise sign operator defined as: sign( x ) = (cid:40) , x ≥ , − , x < . (2)With the assumption that the reflected signal will be the same,i.e., S [: ,
1] = S [: ,
2] = · · · = S [: , M ] = s , within a small time window, for example, within a CPI, the high-precisionmeasurements can be obtained by using the simple digitalintegration (DI) method [8] from the signed measurementmatrix Y . The output of the one-bit system using the DImethod, ˆ s DI , can be written in the following form: ˆ s DI [ n ] = (cid:34) ∆ h M (cid:88) m =1
12 ( Y [ n, m ] + 1) (cid:35) − h − ∆ h, ∆ h = 2 h/ ( M − , n = 1 , . . . , N. (3)The structure of the NVA6100 receiver and the procedure ofthe DI method are shown in Figure 2. The DI method assumes (a) Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e Time Index M a gn i t ud e ⇒⇒⇒⇒⇒ ⇒ (b) Fig. 2: a) Structure of the receiver of a one-bit UWB radarsystem, b) an illustration of the DI method.no interference or weak interferences, and thus it cannotprovide satisfactory performance for strong RFI mitigation.Since strong RFI problems exist in practical applications, aneffective RFI mitigation technique is needed for the properoperations of the one-bit UWB radar systems.
B. Data Model
In the presence of RFI and other noise and disturbances,the signed measurement matrix Y can be written as follows: Y = sign( R θ + S + E − H ) , S [: ,
1] = S [: ,
2] = · · · = S [: , M ] = s , (4) where E denotes the noise and other disturbances and R θ denotes the RFI matrix. By using the fact that the RFI sourcestend to have strong, narrow peaks in the fast-time frequencydomain and the frequencies of the RFI sources change onlyvery slightly over the CPI (see Figure 1), the RFI sources canbe modeled as a sum of sinusoids with their frequencies fixedover the slow-time within the CPI [17]. Thus, each elementof R θ can be expressed as follows: R θ [ n, m ] = K (cid:88) k =1 A k,m sin( ω k ( n −
1) + φ k,m )= K (cid:88) k =1 a k,m cos( ω k ( n − b k,m sin( ω k ( n − n = 1 , · · · , N, m = 1 , · · · , M, (5)where K denotes the number of sinusoids or RFI sources, ω k ∈ [0 , π ) denotes the frequency of the k th RFI source, and A k,m ∈ R + and φ k,m ∈ [0 , π ) denote the amplitude andphase of the k th RFI source during the m th PRI, respectively.The unknown parameter vector of the RFI is denoted by θ =[ a , , b , , . . . , a ,M , b ,M , ω , . . . , a K, , b K, , . . . , a K,M ,b K,M , ω K ] T ∈ R (2 M +1) K with a k,m = A k,m sin φ k,m ∈ R and b k,m = A k,m cos φ k,m ∈ R . Our goal is to recover thedesired radar echo vector s from the signed measurementmatrix Y while mitigating the impact of the RFI.III. 1 B MMRELAX F OR RFI P
ARAMETER E STIMATION
A. Maximum likelihood estimation
We first assume that the desired UWB radar echoes togetherwith the noise and other disturbances, i.e., S + E , obey i.i.d.Gaussian distribution with zero-mean and unknown variance σ . The numerical and experimental examples in Section Vshow that the proposed algorithm is robust to this assumption.We consider the maximum likelihood (ML) estimator forthe RFI parameter estimation problem due to its desirableproperties including consistency and asymptotic efficiency.The ML estimate of the parameter vector β = [ θ T , σ ] T can beobtained by minimizing the following negative log-likelihoodfunction [22]: ˆ˜ β = arg min ˜ β l ( ˜ β )= arg min ˜ β M (cid:88) m =1 N (cid:88) n =1 − log (cid:34) Φ (cid:32) Y [ n, m ] (cid:32) K (cid:88) k =1 ˜ a k,m cos( ω k ( n − b k,m sin( ω k ( n − − λ H [ n, m ] (cid:33)(cid:33)(cid:35) , (6)where Φ( x ) denotes the cumulative distribution func-tion of the standard normal distribution, λ = σ , ˜ a k,m = a k,m σ , ˜ b k,m = b k,m σ , and ˜ β = [ ˜ θ T , λ ] T is the modified unknown parameter vector with ˜ θ =[˜ a , , ˜ b , , . . . , ˜ a ,M , ˜ b ,M , ω , . . . , ˜ a K, , ˜ b K, , . . . , ˜ a K,M , ˜ b K,M , ω K ] T . Note that the cost function in (6) is highly nonlinear and non-convex and is difficult to optimize withrespect to β .Let ω = [ ω , ω , . . . , ω K ] T be a vector composed of the fre-quencies of R θ . For given ω , the above optimization problemis convex with respect to { ˜ a k,m } K,Mk =1 ,M =1 , { ˜ b k,m } K,Mk =1 ,M =1 and λ (see Section 3.54 of [23] for a detailed proof). Thus, the MLproblem can be solved by first performing a K -dimensionalsearch over ω on the feasible space of frequencies [0 , π ) K and then using globally optimal algorithms, for example, theNewton’s Method [23], to compute the corresponding optimalvalues { ˆ˜ a k,m } K,Mk =1 ,M =1 , { ˆ˜ b k,m } K,Mk =1 ,M =1 and ˆ λ as functions of ω . The details on the implementation of the K -dimensionalfrequency search can be found in [24].Note that as the number of sinusoids K and the size ofthe signed measurement matrix increase, this method canbecome computationally prohibitive and thus, more efficientalgorithms, such as the relaxation based methods in [18], [24],may be considered. B. 1bREALX and 1bMMRELAX
Inspired by the RELAX algorithm [25], which is a con-ceptually and computationally simple method for sinusoidalparameter estimation for infinite-precision data, two algo-rithms, 1bRELAX [24] and 1bMMRELAX [18], were recentlyproposed as relaxation-based approaches to estimate the sinu-soidal parameters from single-PRI based signed measurementsby maximizing the likelihood function iteratively. Since 1bRE-LAX and 1bMMRELAX are both cyclic algorithms, their localconvergence is guaranteed under wild conditions [26].We extend below 1bRELAX and 1bMMRELAX to the caseof multiple-PRI based signed measurements. We first extendthe 1bRELAX algorithm to the multiple-PRI cases and thedetailed steps are depicted in Table I. 1bRELAX maximizesthe likelihood function iteratively by conducting a sequenceof one-dimensional frequency searches on [0 , π ) rather thanthe K -dimensional search on [0 , π ) K used by the ML method,and therefore greatly reduces the computational cost. However,since the frequency update procedure is still implemented bymeans of an exhaustive search and a large number of iterationsis required, 1bRELAX is still rather time-consuming.To enhance the computational efficiency of 1bRELAX,the majorization-minimization [27] based 1bMMRELAX isintroduced in [18] for the single-PRI based cases. By usingthe MM technique, 1bMMRELAX transforms the likelihoodmaximization problem in each step into a sequence of simpleinfinite-precision sinusoidal parameter estimation problems,which can be efficiently solved via simple FFT operations. Wenext extend the single-PRI based 1bMMRELAX algorithm tothe multiple-PRI based case, and we still refer to the extendedalgorithm as 1bMMRELAX for simplicity [19].We start by applying the MM-based method to minimize thenegative log-likelihood function l ( ˜ β ) in (6). The majorizingfunction can be obtained by using Lemma 1 in [18]. Theoptimization problem at the ( i + 1) th MM iteration thatminimizes the majorizing function at ˜ β i = [( ˜ θ i ) T , λ i ] T , which TABLE I: 1bRELAX for Multiple-PRI Based Case Input:
Signed measurement matrix Y , the desired or estimatedmodel order ˆ K , and the maximum number of update iteration T R .2: Assume ˜ K = 1 . Obtain {{ ˆ˜ a ,m , ˆ˜ b ,m } Mm =1 , ˆ ω } and ˆ λ bysolving (6) via the exhaustive search (over ω ).3: for ˜ K = 2 : ˆ K t = 0
5: Repeat:6: Obtain {{ ˆ˜ a ˜ K,m , ˆ˜ b ˜ K,m } Mm =1 , ˆ ω ˜ K } by solving (6)via the exhaustive search with {{ ˜ a q,m , ˜ b q,m } Mm =1 , ω q } ˜ K − q =1 and λ replaced by their most recent estimates {{ ˆ˜ a q,m , ˆ˜ b q,m } Mm =1 , ˆ ω q } ˜ K − q =1 and ˆ λ ;7: Redetermine {{ ˆ˜ a ,m , ˆ˜ b ,m } Mm =1 , ˆ ω } and ˆ λ by solving (6)via the exhaustive search with {{ ˜ a q,m , ˜ b q,m } Mm =1 , ω q } ˜ Kq =2 replaced by their most recent estimates {{ ˆ˜ a q,m , ˆ˜ b q,m } Mm =1 , ˆ ω q } ˜ Kq =2 ;8: if ˜ K > for k = 2 : ˜ K −
10: Update {{ ˜ a k,m , ˜ k q,m } Mm =1 , ω k } by by solving (6) viathe exhaustive search with {{ ˜ a q,m , ˜ b q,m } Mm =1 , ω q } ˜ Kq =1 ,q (cid:54) = k and λ replace by their most recent estimates {{ ˆ˜ a q,m , ˆ˜ b q,m } Mm =1 , ˆ ω q } ˜ Kq =1 ,q (cid:54) = k and ˆ λ .11: end end t = t + 1 ;14: Until practical convergence or t reaches the maximum number T R .15: end Output: { ˆ a k,m } ˆ K,Mk =1 ,m =1 = { ˆ˜ a k,m } ˆ K,Mk =1 ,m =1 / ˆ λ , { ˆ b k,m } ˆ K,Mk =1 ,m =1 = { ˆ˜ b k,m } ˆ K,Mk =1 ,m =1 / ˆ λ , { ω k } ˆ Kk =1 ,and ˆ λ . is the estimate of ˜ β obtained at the i th MM iteration, can besimplified as: min ˜ θ ,λ G (cid:16) ˜ θ , λ | ˜ β i (cid:17) = M (cid:88) m =1 N (cid:88) n =1 (cid:104) R ˜ θ [ n, m ] − λ H [ n, m ] − ˜ Z ˜ β i [ n, m ] (cid:105) , (7)where ˜ Z ˜ β i [ n, m ] = Y [ n, m ] (cid:16) X ˜ β i [ n, m ] − f (cid:48) ( X ˜ β i [ n, m ]) (cid:17) is the ( n, m ) th element of an auxiliary matrix ˜ Z ˜ β i , with f ( x ) = − log Φ( x ) , and X ˜ β i [ n, m ] = Y [ n, m ] (cid:0) R ˜ θ i [ n, m ] − λ i H [ n, m ] (cid:1) , n = 1 , · · · , N , m = 1 , · · · , M .The minimization problem in (7) can be conveniently solvedby using a cyclic algorithm [28]. In the cyclic algorithm,we conduct the following two steps iteratively: (1) minimize G (cid:16) ˜ θ , λ | ˜ β i (cid:17) with respect to λ for fixed ˜ θ and (2) minimize G (cid:16) ˜ θ , λ | ˜ β i (cid:17) with respect to ˜ θ for given λ . The closed-formsolution of the first step can be readily obtained as follows: λ i +1 j +1 =max , (cid:80) Mm =1 H T [: , m ] (cid:104) R ˜ θ i +1 j [: , m ] − ˜ Z ˜ β i [: , m ] (cid:105) || H || , (8)where the subscript j denotes the iteration number in thecyclic minimization performed at the i th MM iteration. In thesecond step of the cyclic algorithm, by regarding { V [ n, m ] = λ i +1 j +1 H [ n, m ] + ˜ Z ˜ β i [ n, m ] } N,Mn =1 ,m =1 as the input data, theminimization problem with respect to ˜ θ can be interpretedas the infinite-precision direction-of-arrival (DOA) estimation TABLE II: MM Steps for Solving (6) Input:
Signed measurement matrix Y , known threshold matrix H ,initialization ˜ θ and λ , maximum number of MM iterations T M ,maximum number of inner loop iterations T C , model order K ,and i = 0 .2: Repeat:3: Update: ˜ β i = [ ˜ θ T , λ i ] T ,4: ˜ Z ˜ β i [ n, m ] = Y [ n, m ] (cid:16) X ˜ β i [ n, m ] − f (cid:48) ( X ˜ β i [ n, m ]) (cid:17) n = 1 , . . . , N, m = 1 , . . . , M ,5: ˜ θ i +10 = ˜ θ i +1 , j = 0 , k = K .6: Repeat:7: Update λ i +1 j +1 =max , (cid:80) Mm =1 (cid:80) Nn =1 H [ n,m ] (cid:34) R ˜ θ i +1 j [ n,m ] − ˜ Z ˜ β i [ n,m ] (cid:35)(cid:80) Mm =1 (cid:80) Nn =1 H [ n,m ] ; V k [ n, m ] = ˜ Z ˜ β i [ n, m ] + λ i +1 j +1 H [ n, m ] − (cid:80) Kq =1 ,q (cid:54) = k { ˜ a i +1 q,m,j cos( ω i +1 q,j ( n − b i +1 q,m,j sin( ω i +1 q,j ( n − } n = 1 , . . . , N, m = 1 , . . . , M ;9: Update {{ ˜ a i +1 k,m,j +1 , ˜ b i +1 k,m,j +1 } Mm =1 , ω i +1 k,j +1 } by using infinite-precision RELAX [29] for multiple-PRI based signals on { V k [ n, m ] } N,Mn =1 ,m =1 .10: ˜ θ i +1 j +1 = [˜ a i +11 , ,j +1 , ˜ b i +11 , ,j +1 , . . . , ˜ a i +11 ,M,j +1 , ˜ b i +11 ,M,j +1 , ω i +11 ,j +1 ,. . . , ˜ a i +1 K, ,j +1 , ˜ b i +1 K, ,j +1 , . . . , ˜ a i +1 K,M,j +1 , ˜ b i +1 K,M,j +1 , ω i +1 K,j +1 ] T .11: k = mod( k, K ) + 1 , where mod( · ) denotes the modulo operation;12: j = j + 1 ;13: Until practical convergence or j reaches the maximum number T C .14: i = i + 1 .15: Until practical convergence or i reaches the maximum number T M .16: Output: ˆ˜ θ and ˆ˜ λ . problem encountered in array processing [29]: min ω − M (cid:88) m =1 | a ( ω ) V [: , m ] | , a ( ω ) = [1 , e − jω , . . . , e − jω ( N − ] . (9)Then G (cid:16) ˜ θ , λ i +1 j +1 | ˜ β i (cid:17) can be decreased efficiently by using theinfinite-precision RELAX algorithm in [29] for the multiple-PRI based signals. The infinite-precision RELAX in [29] canbe easily implemented by first using an N -point ( N > N )zero-padded FFT for each PRI and followed by the subsequentfine search of minimizing (9) via the Matlab fminbnd functionover the interval [ˆ ω FFT k − πN , ˆ ω FFT k + πN ] , where ˆ ω FFT k is thefrequency estimate of the k th RFI source obtained by usingFFT. Since there exists a simple closed-form solution for λ ,we re-determine λ after updating each sinusoid.The MM approach for solving the optimization problemin (6) is summarized in Table II. Then, the 1bMMRELAXalgorithm is obtained by replacing the update procedure ofthe extended 1bRELAX (see Steps − of Table I) bythe MM steps proposed in Table II. When initializing theMM approach via the exhaustive search, however, a M -dimensional convex optimization problem needs to be solved N times (see Step 6 of Table I). For large M , this step is stillcomputationally expensive and therefore a faster frequencyinitialization algorithm is needed. C. Fast frequency initialization
To further improve the computational efficiency of theextended 1bMMRELAX, we devise a fast frequency initializa-tion algorithm to estimate the frequencies of the RFI sourcesby exploiting the sparse property of the RFI spectrum.Let { w q = ( q − πQ } Qq =1 denote a grid that covers [0 , π ) .Assuming that the grid is fine enough such that the frequencies (normalized by the sampling frequency) corresponding to theRFI sources are on this grid (or practically, close to the grid)[17], the ( n, m ) th element of the RFI signal R can then berewritten as follows: R [ n, m ] = Q (cid:88) q =1 ˘ a q,m cos( w q ( n − b q,m sin( w q ( n − . (10)Then, denote F = . . . . . . w ) . . . cos( w Q ) sin( w ) . . . sin( w Q ) ... ... ... ... ... ... cos( w ( N − . . . cos( w Q ( N − w ( N − . . . sin( w Q ( N − , (11)and A = ˘ a , · · · ˘ a ,M ... . . . ... ˘ a Q, · · · ˘ a Q,M ˘ b , · · · ˘ b ,M ... . . . ... ˘ b Q, · · · ˘ b Q,M . (12)With the assumption that S + E obeys i.i.d. Gaussian dis-tribution with zero-mean and unknown variance σ , we canestablish an optimization problem as follows, by making useof the group sparsity of the RFI spectrum [17] and the negativelog-likelihood function: min ˜ A ,λ ζ || ˜ A || , + || f ( Y (cid:12) ( F ˜ A − λ H )) || , (13)where λ = 1 /σ , ˜ A = A /σ , (cid:12) denotes the element-wisematrix product, ζ is a user-parameter controlling the balancebetween the two penalty terms, and f ( x ) = − log Φ( x ) .Since the optimization problem in (13) is difficult to solvedirectly, similar to Section III-B, the MM technique is used toobtain a simplified solution. By using the MM technique, theupdate formula at the ( i + 1) th MM iteration can be writtenas: min ˜ A ,λ ζ || ˜ A || , + || F ˜ A − λ H − ˜ Z i FI || , ˜ Z i FI = Y (cid:12) ( ˜ X i FI − f (cid:48) ( ˜ X i FI )) , ˜ X i FI = Y (cid:12) ( F ˜ A i − λ i H ) . (14)Note that in (14), we minimize the majorization functionof (13) at { ˜ A i , λ i } , where { ˜ A i , λ i } is the estimate of { ˜ A , λ } obtained at the i th MM iteration. To solve (14) by ADMM[20] efficiently and effectively, we rewrite (14) as follows: min ˜ A , ˜ B FI ,λ ζ || ˜ A || , + || F ˜ B FI − λ H − ˜ Z i FI || , s . t . ˜ A = ˜ B FI . (15) Then the augmented Lagrangian for (15), with the Lagrangemultiplier Υ FI , is given as: L FI ( ˜ A , ˜ B FI , λ, Υ FI ) = ζ || ˜ A || , + || F ˜ B FI − λ H − ˜ Z i FI || + M (cid:88) m =1 Υ FI [: , m ] T ( ˜ A [: , m ] − ˜ B FI [: , m ]) + ρ FI || ˜ A − ˜ B FI || , (16)where ρ FI > is the penalty parameter. At each iteration ofADMM, we minimize L FI with respect to { ˜ A , ˜ B FI , λ, Υ FI } sequentially. The detailed update steps of ADMM at the ( i +1) th MM iteration are described as follows.
1) Update of ˜ A : The subproblem with respect to ˜ A is: min ˜ A ζ || ˜ A || , + M (cid:88) m =1 Υ i FI [: , m ] T ( ˜ A [: , m ] − ˜ B i FI [: , m ])+ ρ FI || ˜ A − ˜ B i FI || . (17)For given { Υ i FI , ˜ B i FI } , the solution to this subproblem is: ˜ A i +1 = 1 ρ diag( c )( ρ FI ˜ B i FI − Υ i FI ) , c[ n ] = max (cid:32) , − ζ || ρ FI ˜ B i FI [ n, :] − Υ i FI [ n, :] || (cid:33) , (18)where diag( · ) denotes the diagonal matrix whose diagonalelements are composed of the corresponding vector.
2) Update of λ : The subproblem with respect to λ is: min λ || F ˜ B i FI − λ H − ˜ Z i FI || . (19)The solution to this subproblem is: λ i +1 = max (cid:32) , (cid:80) Mm =1 H [: , m ] T ( F ˜ B i FI [: , m ] − ˜ Z i FI [: , m ]) || H || (cid:33) . (20)
3) Update of ˜ B FI : With the estimated { ˜ A i +1 , λ i +1 } and Υ i FI , the subproblem with respect to ˜ B FI is: min ˜ B FI || F ˜ B FI − λ i +1 H − ˜ Z i FI || + M (cid:88) m =1 Υ i FI [: , m ] T ( ˜ A i +1 [: , m ] − ˜ B FI [: , m ])+ ρ FI || ˜ A i +1 − ˜ B FI || . (21)The solution to this subproblem is: ˜ B i +1FI =( F T F + ρ FI I N ) − ( λ i +1 F T H + F T ˜ Z i FI + Υ i FI + ρ FI ˜ A i +1 ) . (22)
4) Update of Υ FI : The Lagrange multiplier Υ FI can beupdated in a manner of gradient descent: Υ i +1FI ← Υ i FI + ρ FI ( ˜ A i +1 − ˜ B i +1FI ) . (23) TABLE III: Fast Frequency Initialization (FI) Input:
Signed measurement matrix Y , known threshold matrix H ,initial values ˜ A , λ , and ε absFI , ε relFI , ε λ FI .2: MM Repeat:3: Initialize ρ FI for ADMM iteration.4: Construct ˜ Z i FI based on ˜ A i , λ i by (14).5: ADMM Repeat:6: Update ˜ A by (18);7: Update λ by (20);8: Update ˜ B FI by (22);9: Update Υ FI by (23);10: Calculate r priFI , r dualFI , r λ FI , ε priFI , ε dualFI by (24) ∼ (26) and (28);11: Update ρ FI by (29).12: Until stopping criteria of ADMM iterations (27) are satisfied.13 Until stopping criteria of MM iterations are satisfied.14: Output: ˆ˜ A FI and ˆ λ FI , ˆ A FI = ˆ λ FI ˆ˜ A FI .
5) Stopping Criteria:
We stop the ADMM iterations whenthe primal, dual residuals and the λ residual are all sufficientlysmall. At the ( j +1) th iteration of ADMM, the primal residualis: r j +1priFI = || ˜ A j +1 − ˜ B j +1FI || , (24)the dual residual is: r j +1dualFI = ρ FI || ˜ B j +1FI − ˜ B j FI || , (25)and the λ residual is: r j +1 λ FI = || λ j +1 − λ j || /λ j +1 . (26)The stopping criteria are: r j priFI ≤ ε j priFI , r j dualFI ≤ ε j dualFI and r jλ FI ≤ ε jλ FI , (27)where ε j priFI , ε j dualFI and ε jλ FI are tolerances for primal, dualand λ residuals for the j th ADMM iteration, respectively.In our problem, ε jλ FI can be set as a small constant and ε j priFI , ε j dualFI can be set as follows: ε j priFI = √ N M ε absFI + ε relFI max( || ˜ A j || , || ˜ B j FI || ) ,ε j dualFI = √ N M ε absFI + ε relFI || Υ j FI || , (28)where ε absFI and ε relFI are absolute relative tolerances, respec-tively, which are set depending on the scale of the problem.
6) Update of ρ FI : To improve the convergence rate andreduce the influence of the initial value of the penalty param-eter ρ FI , we update the penalty parameter for each ADMMiteration using the following updating scheme [20]: ρ j +1FI = ρ j FI , r j priFI > r j dualFI , . ρ j FI , r j dualFI > r j priFI ,ρ j FI , otherwise , (29)where ρ j FI denotes the penalty parameter in j th ADMMiteration.The fast frequency initialization can be sum-marized in Table III. From ˆ A FI , we obtain ˚ A = (cid:113) ˆ A FI [1 : Q, :] + ˆ A FI [ Q + 1 : 2 Q, :] , where thesquare and square root are both element-wise operations.Suppose that the maximum possible model order is K max .Since the frequencies of the RFI sources are known to benot close to zero, the K max fast frequency initializations ofthe RFI sources, referred to as { ˆ ω FI k } K max k =1 , correspond to the K max row indices of ˚ A [2 : Q, :] (note that the row index corresponds to the zero frequency) with the K max largest (cid:96) norms. The K max initial frequency estimates are then sortedbased on arranging the corresponding (cid:96) norms in descendingorder.For the ˜ K th step of 1bMMRELAX, the coarse initialfrequency estimate of the ˜ K th sinusoid is taken as ˆ ω FI˜ K , in lieuof the one obtained via the exhaustive search. Also, ˆ λ FI can beused as an initial value of λ in the first step of 1bMMRELAX. D. 1bBIC
For the case that the number of RFI sources, i.e., the truemodel order K , is unknown, we consider extending the single-PRI based one-bit Bayesian information criterion (1bBIC) in[21] to the multiple-PRI based cases so that the extended1bBIC can be used with the extended 1bMMRELAX todetermine the RFI model order. Suppose that ˆ β ˜ K is parametervector estimated by the extended 1bMMRELAX algorithm foran assumed model order ˜ K , where ≤ ˜ K ≤ K max . It isproven in Appendix A that the 1bBIC cost function for themultiple-PRI based signed measurements has the followingform: K ) = − M (cid:88) m =1 N (cid:88) n =1 log (cid:32) Φ (cid:32) Y [ n, m ] ˆ R ˆ θ ˜ K [ n, m ] − H [ n, m ]ˆ σ (cid:33)(cid:33) + ˜ K (3 + 2 M ) log N. (30)The estimate ˆ K of K is selected as the integer that minimizesthe above 1bBIC cost function with respect to the assumednumber of sinusoids ˜ K .IV. R ADAR E CHO S IGNAL R ECOVERY
Due to the sparsity of strong targets in a scene of interest,the desired UWB radar echo vector s is in general quite sparse.Thus, using the estimated RFI sources, s can be recoveredby exploiting its sparsity and minimizing the correspondingnegative log-likelihood function. Since the desired signal isinvariant over all PRIs within the CPI, (4) can be rewrittenapproximately as follows: Y [: , m ] ≈ sign( D γ − U [: , m ]) , m = 1 , . . . , M, U = H − ˆ R ˆ θ ∈ R N × M , (31)where ˆ R ˆ θ is the RFI estimate, D ∈ R N × N denotes thedictionary whose columns are time-shifted digitized versionsof the transmitted impulse, and γ ∈ R N is a sparse vectorcontaining the information of the magnitudes and positions ofthe radar echoes. The estimate of γ can be obtained by solvingthe following convex optimization problem [30]: min ˜ γ ,λ ζ M || ˜ γ || + M (cid:88) m =1 || f ( Y [: , m ] (cid:12) ( D ˜ γ − λ U [: , m ])) || , (32)where ζ is a user-parameter and ˜ γ = γ /σ, λ = 1 /σ , and f ( x ) = − log Φ( x ) . Similar to (13), the convex objective function (32) can beminimized efficiently by using the MM and ADMM tech-niques. The optimization problem at the ( i +1) th MM iterationof (32) can be written as: min ˜ γ ,λ ζ M || ˜ γ || + M (cid:88) m =1 || D ˜ γ − λ U [: , m ] − ˜ Z i ER [: , m ] || (33)where ˜ Z i ER = Y (cid:12) ( ˜ X i ER − f (cid:48) ( ˜ X i ER ))˜ X i ER [: , m ] = Y [: , m ] (cid:12) ( D ˜ γ i − λ i U [: , m ]) , m = 1 , , . . . , M. (34)Here, { ˜ γ i , λ i } denotes the estimate of { ˜ γ , λ } at the i th MMiteration. Then, (33) can be solved effectively and efficientlyby ADMM. To use ADMM technique, we rewrite (14) asfollows: min ˜ γ , ˜ B ER ,λ ζ M || ˜ γ || + M (cid:88) m =1 || D ˜ B ER [: , m ] − λ U [: , m ] − ˜ Z i ER [: , m ] || s . t . ˜ γ = ˜ B ER [: , m ] , m = 1 , , . . . , M. (35)Then the augmented Lagrangian for (35) with the Lagrangemultiplier Υ ER is given as: L ER (˜ γ , ˜ B ER , λ, Υ ER ) = ζ M || ˜ γ || + M (cid:88) m =1 || D ˜ B ER [: , m ] − λ U [: , m ] − ˜ Z i ER [: , m ] || + M (cid:88) m =1 Υ ER [: , m ] T (˜ γ − ˜ B ER [: , m ])+ ρ ER M (cid:88) m =1 || ˜ γ − ˜ B ER [: , m ] || , (36)where ρ ER > is the penalty parameter and in each iterationof ADMM, we minimize L ER by updating { ˜ γ , ˜ B ER , λ, Υ ER } sequentially. The detailed update steps are described as fol-lows.
1) Update of ˜ γ : Given { Υ i ER , ˜ B i ER } , the subproblem withrespect to ˜ γ is: min ˜ γ ζ M || ˜ γ || + M (cid:88) m =1 Υ i ER [: , m ] T (˜ γ − ˜ B i ER [: , m ])+ ρ ER M (cid:88) m =1 || ˜ γ − ˜ B i ER [: , m ] || . (37)The solution to this subproblem is: ˜ γ i +1 =softthreshold (cid:32) M M (cid:88) m =1 (cid:18) ˜ B i ER [: , m ] − ρ ER Υ i ER [: , m ] (cid:19) , ζ ρ ER (cid:33) , (38)where softthreshold( x, a ) = x + a, x ≤ − a,x − a, x ≥ a, , otherwise . (39)
2) Update of λ : The subproblem with respect to λ is: min λ M (cid:88) m =1 || D ˜ B i ER [: , m ] − λ U [: , m ] − ˜ Z i ER [: , m ] || . (40)The solution to this subproblem is: λ i +1 = max (cid:32) , (cid:80) Mm =1 U [: , m ] T ( D ˜ B i ER [: , m ] − ˜ Z i ER [: , m ]) || U || (cid:33) . (41)
3) Update of ˜ B ER : The subproblem with respect to ˜ B ER is: min ˜ B ER M (cid:88) m =1 || D ˜ B ER [: , m ] − λ i +1 U [: , m ] − ˜ Z i ER [: , m ] || + M (cid:88) m =1 Υ i ER [: , m ] T (˜ γ i +1 − ˜ B ER [: , m ])+ ρ ER M (cid:88) m =1 || ˜ γ i +1 − ˜ B ER [: , m ] || . (42)For given { ˜ γ i +1 , λ i +1 , Υ i ER } , the solution to this subproblemis: ˜ B ER =( D T D + ρ ER I N ) − ( λ i +1 D T U + D T ˜ Z i ER + Υ i ER + ρ ER ˜ Γ i +1 ) , ˜ Γ = [˜ γ , . . . , ˜ γ ] ∈ R N × M . (43)
4) Update of Υ ER : The Lagrange multiplier Υ ER can beupdated via gradient descent: Υ i +1ER [: , m ] ← Υ i ER [: , m ] + ρ ER (˜ γ i +1 − ˜ B i +1ER [: , m ]) ,m = 1 , , . . . , M. (44)
5) Stopping Criteria:
Similar to Section III-C5, we stop theADMM iterations when the primal, dual residuals and the λ residual are all sufficiently small. At the ( j + 1) th iteration ofADMM, the three values can be expressed as: primal residual : r j +1priER = || ˜ Γ j +1 − ˜ B j +1ER || , dual residual : r j +1dualER = ρ ER || ˜ B j +1ER − ˜ B j ER || ,λ residual : r j +1 λ ER = || λ j +1 − λ j || /λ j +1 . (45)Then the stopping criteria are: r j priER ≤ ε j priER , r j dualER ≤ ε j dualER and r jλ ER ≤ ε jλ ER , (46)where ε jλ ER can be set as a small constant and ε j priER , ε j dualER can be set as follows: ε j priER = √ N M ε absER + ε relER max( || ˜ Γ j || , || ˜ B j ER || ) ,ε j dualER = √ N M ε absER + ε relER || Υ j ER || , (47)with ε absER and ε relER being absolute and relative tolerances,respectively, which are chosen depending on the scale of theproblem. TABLE IV: Echo Recovery Input:
Signed measurement matrix Y , known threshold matrix H ,estimated RFI signal ˆ R ˆ θ , initial values ˜ γ and ε absER , ε relER , ε λ ER .The initial value of λ can be set as the value obtained by 1bMMREALX.2: U = H − ˆ R ˆ θ .3: MM Repeat:4: Initialize ρ ER for ADMM iteration.5: Construct ˜ Z i ER based on ˜ B i ER , λ i by (33).6: ADMM Repeat:7: Update ˜ γ by (38);8: Update λ by (41);9: Update ˜ B ER by (43);10: Update Υ ER by (44);11: Calculate r priER , r dualER , r λ ER , ε priER , ε dualER by (45) and (47);12: Update ρ ER by (48).13: Until stopping criteria of ADMM iterations (46) are satisfied.14 Until stopping criteria of MM iterations are satisfied.15: Output: ˆ˜ γ and ˆ λ . Recovered Radar Echo ˆ s = D ˆ˜ γ / ˆ λ .
6) Update of ρ ER : Similar to Section III-C6, we update ρ ER by the following scheme: ρ j +1ER = ρ j ER , r j priER > r j dualER , . ρ j ER , r j dualER > r j priER ,ρ j ER , otherwise , (48)where ρ j ER denotes the penalty parameter in the j th ADMMiteration. The echo recovery algorithm is summarized in TableIV. V. S IMULATED AND E XPERIMENTAL E XAMPLES
In this section, we evaluate the RFI mitigation performanceof the proposed algorithm for one-bit UWB radar systemsusing both simulated and measured RFI data sets. Hereafter,the proposed algorithm includes using the fast frequency ini-tialization with the extended 1bMMRELAX for RFI parameterestimation and the extended 1bBIC for RFI order determi-nation (see Section III) and the sparse method to recoverthe radar echoes (see Section IV). We conduct experimentsusing simulated RFI-free UWB radar data and two differentRFI data sets: simulated RFI data set and measured RFI dataset. The measured RFI data set was collected by the ARLexperimental radar receiver from a real-world environmentwith the antenna pointing toward Washington DC (see [31]–[33] for more details about the data collection using the ARLradar). Because the sampling rate of the ARL radar receiveris GHz, we assume that all data sets used in this section areobtained with an GHz sampling rate.All data sets contain slow-time samples within a CPIand fast-time samples per PRI, i.e., M = 8192 , N = 512 .The transmitted radar impulse is shaped as the first orderderivative of a Gaussian pulse with samples coveringthe frequency range of ∼ MHz (see Figures 3(a)and 3(b)). The simulated RFI-free and noise-free UWB radarechoes are generated by targets at different ranges withdifferent amplitudes (see Figure 3(c)).All signed measurements are obtained by sampling the RFI-contaminated data via the CTBV sampling technology. The DImethod described in Section II-A is used as a benchmark. Theextended 1bMMRELAX algorithm without the fast frequencyinitialization is not considered herein because of its computa-tionally prohibitive complexity. Fast-time Index -1.5-1-0.500.511.5 M a gn i t ud e (a) Frequency (GHz) -3.5-3-2.5-2-1.5-1-0.50 M a gn i t ud e ( d B ) (b)
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e (c) Fig. 3: (a) Simulated transmitted radar pulse; (b) Spectrum ofthe transmitted radar pulse; (c) Simulated RFI-free and noise-free radar echoes for one PRI.
A. Evaluation Metric
Note that the measured RFI data set inevitably containsnoise and other disturbances whereas the simulated RFI dataset is noise-free. We add white Gaussian noise E to thesimulated RFI data set, but we do not add extra noise tothe measured (and hence already noisy) RFI data set. Dif-ferent interference-to-noise ratios (INRs) of the simulated RFIpowers to the additive Gaussian noise powers, measured by
20 log (cid:107) R (cid:107) (cid:107) E (cid:107) (dB), will be considered in Section V-C to showthe performance of the proposed algorithm.Denote the signal-to-interference-plus-noise ratio (SINR)for the simulated RFI data sets and radar echoes as follows: SINR = 20 log || S || || R + E || (dB) . (49)For the case of simulated radar echoes and measured RFIdata sets, with the measured RFI already containing noise,the SINR is computed as follows: SINR = 20 log || S || || R || (dB) . (50)We fix the desired radar echo signal and add the scaledsimulated RFI plus noise or the scaled measured RFI toobtain contaminated data sets with various SINR values. Themaximum threshold h is set as for all cases according tothe magnitude of the desired radar echo signal in Figure 3(c).The radar echo signal recovery performance is measuredusing the normalized recovery error (NRE) of the recoveredsignal vector: NRE = 20 log || s − ˆ s || || s || (dB) , (51)where ˆ s is the recovered UWB radar echo signal. B. Implementation Details
In our implementation of 1bMMRELAX, the MM iterationsare terminated if the relative change of the negative log-likelihood function l ( ˜ β ) between two consecutive iterationsis less than − or a maximum number of the MM it-erations T M = 20 is reached. Within each MM iteration,we terminate the inner loop if the relative change of theobjective function in (7) is less than − or a maximumnumber of iterations T C = 20 is reached. When using the N -point FFT in the infinite-precision RELAX, N is set to N . When using the fast frequency initialization, we set Q = N, ε absFI = 10 − , ε relFI = 10 − and ε λ FI = 10 − ,and set the user-parameter ζ = 1 for all cases. We terminatethe inner loop, i.e., the ADMM iteration, when the stoppingcriteria in (27) is satisfied or the iteration number is larger than10. The outer loop, i.e., the MM iteration, will be terminatedif the relatively change of the values of the objective function(13) is smaller than − or the iteration number is largerthan 5. In the final radar echo signal recovery step, we set ε absER = 10 − , ε relER = 10 − and ε λ ER = 10 − , and ζ = 0 . for all cases. We terminate the inner loop, i.e., theADMM iteration, when the stopping criteria in (46) is satisfiedor the iteration number is larger than 100. The outer loop, i.e.,the MM iteration, will be terminated if the relative change of the values of the objective function (32) is smaller than − or the iteration number is larger than 50. C. Simulated RFI cases
In this section, we present the results of the proposedalgorithm for the simulated radar echoes contaminated by thesimulated RFI. The magnitudes of the RFI sources usuallydo not change greatly among different PRIs within a CPI.Thus we simulate the RFI sources as a sum of sinusoids withamplitudes and frequencies fixed from one PRI to anotherwithin a CPI and the phases varying independently (withuniform distribution between and π ) with slow-time. Thesimulated RFI can be written as follows: R simulated [ n, m ] = K (cid:88) k =1 A k sin( ω k ( n −
1) + φ k,m ) (52)where A k A , k = 2 , . . . , K , is the amplitude ratio of the k -th RFIsource relative to the first one. The detailed parameter settingsof the simulated RFI data set are shown in Table V and Figure4. When generating the contaminated data sets with differentSINR values, the desired RFI can be obtained by varying A while fixing { A k A } Kk =2 .TABLE V: Simulated RFI Parameter Settings RFI Frequencies (MHz)
500 350 700 900 1050
RFI Amplitude Ratios .
95 0 . .
87 0 . Slow-time Index F r e qu ec n y ( G H z ) -50-45-40-35-30-25-20-15-10-50(dB) Fig. 4: Spectrum of the simulated RFI.The RFI mitigation results for the simulated RFI data setobtained by the proposed algorithm and the DI method areshown in Figures 5 - 7. It is clear that the proposed algorithmoutperforms the DI method as SINR varies from − dBto − dB. Additionally, for the challenging cases of SINR = − dB, for example, most of the targets are retained bythe proposed algorithm, while few targets are distinguishableby the DI method. The results also show that using theextended 1bMMRELAX with the extended 1bBIC providesaccurate model order estimates for the multiple-PRI basedsigned measurements for a wide range of SINR values. Notethat the choices of the user parameters ζ and ζ work welleven as SINR varies from − dB to − dB. -40 -35 -30 SINR (dB) -14-12-10-8-6-4-20 O u t pu t NR E ( d B ) Proposed AlgorithmDI Method (a) -40 -35 -30
SINR (dB) -14-12-10-8-6-4-20 O u t pu t NR E ( d B ) Proposed AlgorithmDI Method (b)
Fig. 5: Results of the proposed method and the DI method forsimulated RFI data sets with fixed ζ = 0 . . The INR is a) dB and b) dB. D. Measured RFI cases
In this section, we present RFI mitigation results using theaforementioned simulated RFI-free UWB radar echoes and themeasured RFI data set collected by the ARL experimentalradar receiver [31], [32]. We use the first × subset ofthe original RFI data set measured by the ARL radar receiver.The spectrum of the measured RFI data set is shown in Figure1. The RFI mitigation results for the measured RFI data setby the proposed algorithm and the DI method are shown inFigures 8 - 9. Over a wide range of SINR values, the proposedalgorithm outperforms the DI method. More specifically, forthe challenging case of SINR = − dB, most of the targetsare retained by the proposed algorithm, while few targets arediscernible by the DI method.VI. C ONCLUSIONS
In this paper, we have established a RFI mitigation frame-work for a one-bit UWB radar system that obtains its signedmeasurements by using the CTBV sampling technique. Wehave established a proper data model for the RFI sourcesand have extended a majorization-minimization based 1bMM-RELAX algorithm to efficiently estimate the RFI parameters
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e Original Radar Echo SignalRecovered Radar Echo Signal by Proposed Algorithm (a)
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e Original Radar Echo SignalRecovered Radar Echo Signal by DI Method (b)
Fig. 6: Radar echo recovery results for the simulated RFIobtained by using a) the proposed algorithm and b) the DImethod, when SINR = − dB and INR = 10 dB.from multiple-PRI based signed measurements. We have alsointroduced a fast frequency initialization algorithm, based onthe MM and ADMM techniques, to further reduce the com-putational cost of the extended 1bMMRELAX. Moreover, wehave extended 1bBIC to the multiple-PRI based cases so thatthe extended 1bMMRELAX can be used with the extended1bBIC to simultaneously estimate the RFI parameters anddetermine the number of RFI sources. Next, by exploitingthe sparse property of the UWB radar echoes, a sparsesignal recovery method, also based on the MM and ADMMtechniques, has been devised to estimate the desired UWBradar echoes using the estimated RFI parameters. Finally, wehave provided examples of using both simulated and measuredRFI data sets to demonstrate that the proposed algorithm canoutperform the benchmark DI method significantly for radarecho recovery, especially for the severe RFI cases.A PPENDIX AD ERIVATION OF B BIC
FOR M ULTIPLE -PRI B
ASED S IGNED M EASUREMENTS
Rewrite the signed measurement matrix Y as follows: Y = sign( R ˘ θ + S + E − H ) ∈ R N × M , (53)
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e Original Radar Echo SignalRecovered Radar Echo Signal by Proposed Algorithm (a)
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e Original Radar Echo SignalRecovered Radar Echo Signal by DI Method (b)
Fig. 7: Radar echo recovery results for the simulated RFIobtained by using a) the proposed algorithm and b) the DImethod, when SINR = − dB and INR = 0 dB. -40 -35 -30 SINR (dB) -10-9-8-7-6-5-4-3-2-10 O u t pu t NR E ( d B ) Proposed AlgorithmDI Method
Fig. 8: Results of the proposed algorithm and the DI methodfor the measured RFI data sets.
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e Original Radar Echo SignalRecovered Radar Echo Signal by Proposed Algorithm (a)
50 100 150 200 250 300 350 400 450 500
Fast-time Index -400-300-200-1000100200300400 M a gn i t ud e Original Radar Echo SignalRecovered Radar Echo Signal by DI Method (b)
Fig. 9: Radar echo recovery results for the measured RFIobtained by using a) the proposed algorithm and b) the DImethod, when SINR = − dB.where R ˘ θ [ n, m ] = K (cid:88) k =1 A k,m sin( ω k ( n −
1) + φ k,m ) , (54)with ˘ θ = [ ω , . . . , ω K , A , , . . . , A K,M , φ , , . . . , φ K,M ] T ∈ R (2 M +1) K denoting the unknown parameter vector of theRFI sources. Assuming that S + E obeys i.i.d. Gaussiandistribution with zero-mean and unknown variance σ , theunknown parameter vector of Y is denoted by ˘ β = [ ˘ θ , σ ] T .By assuming that the prior probability density function(PDF) of ˘ β is flat around the ML estimate and independentof N and M , the estimate of the signal model order K canbe obtained by minimizing the following criterion [34]: l Y (ˆ˘ β ) + log det(ˆ J ( Y , ˆ˘ β )) , (55)where l Y (ˆ˘ β ) denotes the negative log-likelihood function of Y and ˆ˘ β is the ML estimate of ˘ β . With the aforementionedGaussian distribution assumption, l Y (ˆ˘ β ) can be expressed asfollows: l Y (ˆ˘ β ) = − M (cid:88) m =1 N (cid:88) n =1 log (cid:20) Φ (cid:18) Y [ n, m ] R ˘ θ [ n, m ] − H [ n, m ] σ (cid:19)(cid:21) . (56) The function det( · ) means taking the determinant of a matrix.The matrix ˆ J ( Y , ˆ˘ β ) in (55) is defined as: ˆ J ( Y , ˆ˘ β ) = ∂ l Y (ˆ˘ β ) ∂ ˘ β ∂ ˘ β T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˘ β =ˆ˘ β . (57)Under mild conditions (see, e.g., [21], [34] and the referencestherein), the matrix ˆ J ( Y , ˆ˘ β ) has the following asymptoticrelationship with the Fisher information matrix (FIM) J forthe parameter vector ˘ β [21]: [ P N JP N − P N ˆ J ( Y , ˆ˘ β ) P N ] → , as N → ∞ , (58)where P N is a normalization matrix, and J = E (cid:40) ∂ l Y (ˆ˘ β ) ∂ ˘ β ∂ ˘ β T (cid:41) = E (cid:40) ∂l Y (ˆ˘ β ) ∂ ˘ β ∂l Y (ˆ˘ β ) ∂ ˘ β T (cid:41) . (59)Hence, after proper normalizations, ˆ J ( Y , ˆ˘ β ) in (55) can besubstituted by the FIM in the large sample (i.e., N (cid:29) ) case,a fact that can be used to obtain a simple asymptotically validexpression for the penalty term in (55).Let u n,m = [ ∂ R ˘ θ [ n,m ] ∂ ˘ θ T , − R ˘ θ [ n,m ] − H [ n,m ] σ ] T . After simplecalculations, the FIM can be expressed as a sum of positivesemidefinite matrices [22]: J = 12 πσ N (cid:88) n =1 M (cid:88) m =1 ξ (cid:18) R ˘ θ [ n, m ] − H [ n, m ] σ (cid:19) u n,m u Tn,m , (60)where ξ ( x ) is defined as follows: ξ ( x ) = (cid:20) x ) + 1Φ( − x ) (cid:21) e − x . (61)Let J = σ (cid:80) Nn =1 (cid:80) Mm =1 u n,m u Tn,m . Note that log det( J ) and log det( J ) are asymptotically equivalent (see [21] for adetailed proof), to within a constant that does not depend on N . Therefore, the FIM J can be substituted by the simplermatrix J . Note that J is not equal to the FIM for theconventional BIC. Next, we analyze J and the normalizationmatrix P N needed for P N J P N to be O (1) for our RFI signalmodel.Under the weak assumption that the threshold matrix H isuncorrelated with the RFI R ˘ θ , we have the following results(see [21] and Appendix A of [35]): lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂ω k i ∂ R ˘ θ [ n, m ] ∂ω k j = lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ( n − A k i ,m A k j ,m cos( ω k i ( n −
1) + φ k i ,m ) cos( ω k j ( n −
1) + φ k j ,m )= δ i,j σ M (cid:88) m =1 A k i ,m ,δ i,j = (cid:40) , i (cid:54) = j , i = j ≤ i, j ≤ K, (62) lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂ω k i ∂ R ˘ θ [ n, m ] ∂A k j ,m (cid:48) = lim N →∞ σ N N (cid:88) n =1 ( n − A k i ,m (cid:48) cos( ω k i ( n −
1) + φ k i ,m (cid:48) ) sin( ω k j ( n −
1) + φ k j ,m (cid:48) )= 0 , ≤ m (cid:48) ≤ M, (63) lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂ω k i ∂ R ˘ θ [ n, m ] ∂φ k j ,m (cid:48) = lim N →∞ σ N N (cid:88) n =1 ( n − A k i ,m (cid:48) A k j ,m (cid:48) cos( ω k i ( n −
1) + φ k i ,m (cid:48) ) cos( ω k j ( n −
1) + φ k j ,m (cid:48) )= δ i,j σ A k i ,m (cid:48) , (64) lim N →∞ − σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂ω k i R ˘ θ [ n, m ] − H [ n, m ] σ = lim N →∞ − σ N N (cid:88) n =1 M (cid:88) m =1 ( n − A k i ,m cos( ω k i ( n −
1) + φ k i ,m )( R ˘ θ [ n, m ] − H [ n, m ])= 0 , (65) lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂A k i ,m p ∂ R ˘ θ [ n, m ] ∂A k j ,m q , ≤ p, q ≤ M = lim N →∞ σ N N (cid:88) n =1 sin( ω k i ( n −
1) + φ k i ,m p )sin( ω k j ( n −
1) + φ k j ,m p ) δ p,q = 12 σ δ i,j δ p,q , (66) lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂A k i ,m p ∂ R ˘ θ [ n, m ] ∂φ k j ,m q = lim N →∞ σ N N (cid:88) n =1 A k j ,m p sin( ω k i ( n −
1) + φ k i ,m p )cos( ω k j ( n −
1) + φ k j ,m p ) δ p,q = 0 , (67) lim N →∞ − σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂A k i ,m (cid:48) R ˘ θ [ n, m ] − H [ n, m ] σ = lim N →∞ − σ N N (cid:88) n =1 sin( ω k i ( n −
1) + φ k i ,m (cid:48) )( R ˘ θ [ n, m (cid:48) ] − H [ n, m (cid:48) ])= lim N →∞ − σ N N (cid:88) n =1 (cid:34) sin( ω k i ( n −
1) + φ k i ,m (cid:48) ) K (cid:88) k =1 A k,m (cid:48) sin( ω k ( n −
1) + φ k,m (cid:48) ) (cid:35) = − A k i ,m (cid:48) σ , (68) lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂φ k i ,m p ∂ R ˘ θ [ n, m ] ∂φ k j ,m q = lim N →∞ σ N N (cid:88) n =1 A k i ,m p A k j ,m p cos( ω k i ( n −
1) + φ k i ,m p )cos( ω k j ( n −
1) + φ k j ,m p ) δ p,q = A k i ,m p σ δ i,j δ p,q , (69) lim N →∞ − σ N N (cid:88) n =1 M (cid:88) m =1 ∂ R ˘ θ [ n, m ] ∂φ k i ,m (cid:48) R ˘ θ [ n, m ] − H [ n, m ] σ = lim N →∞ − σ N N (cid:88) n =1 A k i ,m (cid:48) cos( ω k i ( n −
1) + φ k i ,m (cid:48) )( R ˘ θ [ n, m (cid:48) ] − H [ n, m (cid:48) ])= 0 , (70)and lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 (cid:18) R ˘ θ [ n, m ] − H [ n, m ] σ (cid:19) = lim N →∞ σ N N (cid:88) n =1 M (cid:88) m =1 (cid:0) R ˘ θ [ n, m ] − H [ n, m ] (cid:1) = 1 σ (cid:34) K (cid:88) k =1 M (cid:88) m =1 A k,m + σ H (cid:35) ,σ H = lim N →∞ N N (cid:88) n =1 M (cid:88) m =1 H [ n, m ] . (71)Define the normalization matrix P N as P N = (cid:20) N / I K N / I KM +1 (cid:21) . (72) Using (62)-(71), we have that lim N →∞ P N J P N = 12 σ ˜ J ωω J ωφ I KM J Aσ ˜ J Tωφ J φφ
00 ˜ J TAσ (cid:80) Kk =1 (cid:80) Mm =1 A k,m +2 σ H σ , ˜ J ωω = (cid:80) Mm =1 A ,m . . . (cid:80) Mm =1 A K,m ∈ R K × K , ˜ J ωφ = A , · · · A ,M · · · · · · . . . . . . . . . . . . . . . · · · · · · A K, · · · A K,M ∈ R K × KM , ˜ J Aσ = [ − A , σ , . . . , − A K,M σ ] T ∈ R KM , (73) ˜ J φφ = A , . . . A K,M ∈ R KM × KM , which implies, after straightforward calculations, that det( lim N →∞ P N J P N ) = 112 K (cid:18) σ (cid:19) K (1+2 M )+1 σ H σ K (cid:89) k =1 (cid:32) M (cid:88) m =1 A k,m (cid:33) K (cid:89) k =1 M (cid:89) m =1 A k,m . (74)Note that the asymptotic matrix in (73) is nonsingular if andonly if σ H (cid:54) = 0 (i.e., the threshold is not zero). Combining (58)and (74), we get log det(ˆ J ( Y , ˆ˘ β )) = log det( P − N ) + log det( P N J P N )= ( K (3 + 2 M ) + 1) log N + O (1) . (75)Therefore, the extended 1bBIC metric is given by − M (cid:88) m =1 N (cid:88) n =1 log (cid:20) Φ (cid:18) Y [ n, m ] R ˘ θ [ n, m ] − H [ n, m ] σ (cid:19)(cid:21) + K (3 + 2 M ) log N. (76)where we only keep the terms that depend on K .R EFERENCES[1] C. Chen, M. B. Higgins, K. O’Neill, and R. Detsch, “Ultrawide-bandwidth fully-polarimetric ground penetrating radar classification ofsubsurface unexploded ordinance,”
IEEE Transactions on Geoscienceand Remote Sensing , vol. 39, no. 6, pp. 1221–1230, June 2001.[2] X. Xu and R. M. Narayanan, “FOPEN SAR imaging using UWBstep-frequency and random noise waveforms,”
IEEE Transactions onAerospace and Electronic Systems , vol. 37, no. 4, pp. 1287–1300, Oct2001.[3] A. G. Yarovoy, L. P. Ligthart, J. Matuzas, and B. Levitas, “UWB radarfor human being detection,”
IEEE Aerospace and Electronic SystemsMagazine , vol. 21, no. 3, pp. 10–14, March 2006.[4] X. Shang, J. Liu, and J. Li, “Multiple object localization and vitalsign monitoring using IR-UWB MIMO radar,”
IEEE Transactions onAerospace and Electronic Systems , vol. 56, no. 6, pp. 4437–4450, 2020. [5] B. Zhao, L. Huang, J. Li, M. Liu, and J. Wang, “Deceptive SAR jammingbased on 1-bit sampling and time-varying thresholds,”
IEEE Journal ofSelected Topics in Applied Earth Observations and Remote Sensing ,vol. 11, no. 3, pp. 939–950, March 2018.[6] B. Zhao, L. Huang, and W. Bao, “One-bit SAR imaging based on single-frequency thresholds,”
IEEE Transactions on Geoscience and RemoteSensing , May 2007, pp.1210–1213.[9] T. Koutsoudis and L. A. Lovas, “RF interference suppression in ultra-wideband radar receivers,” in
Algorithms for Synthetic Aperture RadarImagery II , vol. 2487. International Society for Optics and Photonics,1995, pp. 107–119.[10] D. O. Carhoun, “Adaptive nulling and spatial spectral estimation usingan iterated principal components decomposition,” in , April 1991,pp. 3309–3312 vol.5.[11] H. Subbaram and K. Abend, “Interference suppression via orthogonalprojections: a performance analysis,”
IEEE Transactions on Antennasand Propagation , vol. 41, no. 9, pp. 1187–1194, Sep. 1993.[12] V. T. Vu, T. K. Sj¨ogren, M. I. Pettersson, and L. H˚akasson, “An approachto suppress RFI in ultrawideband low frequency SAR,” in , May 2010, pp. 1381–1385.[13] X. Luo, L. M. H. Ulander, J. Askne, G. Smith, and P. O. Frolind,“RFI suppression in ultra-wideband SAR systems using LMS filters infrequency domain,”
Electronics Letters , vol. 37, no. 4, pp. 241–243, Feb2001.[14] T. Miller, L. Potter, and J. McCorkle, “RFI suppression for ultra wide-band radar,”
IEEE Transactions on Aerospace and Electronic Systems ,vol. 33, no. 4, pp. 1142–1156, Oct 1997.[15] F. Zhou, R. Wu, M. Xing, and Z. Bao, “Eigensubspace-based filteringwith application in narrow-band interference suppression for SAR,”
IEEE Geoscience and Remote Sensing Letters , vol. 4, no. 1, pp. 75–79, Jan 2007.[16] F. Zhou, M. Tao, X. Bai, and J. Liu, “Narrow-band interferencesuppression for SAR based on independent component analysis,”
IEEETransactions on Geoscience and Remote Sensing , vol. 51, no. 10, pp.4952–4960, Oct 2013.[17] J. Ren, T. Zhang, J. Li, L. H. Nguyen, and P. Stoica, “RFI mitigationfor UWB radar via hyperparameter-free sparse SPICE methods,”
IEEETransactions on Geoscience and Remote Sensing , vol. 57, no. 6, pp.3105–3118, June 2019.[18] J. Ren, T. Zhang, J. Li, and P. Stoica, “Sinusoidal parameter estima-tion from signed measurements via majorization-minimization basedRELAX,”
IEEE Transactions on Signal Processing , vol. 67, no. 8, pp.2173–2186, April 2019.[19] T. Zhang, J. Ren, C. Gianelli, and J. Li, “RFI mitigation for one-bitUWB radar systems,” in , 2019, pp. 1545–1549.[20] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributedoptimization and statistical learning via the alternating direction methodof multipliers,”
Foundations and Trends in Machine Learning , vol. 3,no. 1, pp. 1–122, 2011.[21] C. Li, R. Zhang, J. Li, and P. Stoica, “Bayesian information criterionfor signed measurements with application to sinusoidal signals,”
IEEESignal Processing Letters , vol. 25, no. 8, pp. 1251–1255, Aug 2018.[22] C. Gianelli, L. Xu, J. Li, and P. Stoica, “One-bit compressive samplingwith time-varying thresholds: Maximum likelihood and the Cram´er-Raobound,” in , Nov 2016, pp. 399–403.[23] S. Boyd and L. Vandenberghe,
Convex optimization . Cambrige,UK: Cambridge university press, 2009. [Online]. Available: http://stanford.edu/ ∼ {} boyd/cvxbook/bv cvxbook.pdf[24] C. Gianelli, L. Xu, J. Li, and P. Stoica, “One-bit compressive samplingwith time-varying thresholds for multiple sinusoids,” in , Dec 2017, pp. 1–5.[25] J. Li and P. Stoica, “Efficient mixed-spectrum estimation with ap-plications to target feature extraction,” IEEE Transactions on SignalProcessing , vol. 44, no. 2, pp. 281–295, Feb 1996.[26] W. I. Zangwill,
Nonlinear programming: A unified approach . Engle-wood Cliffs, N.J: Prentice-Hall, 1969. [27] D. Hunter and K. Lange, “A tutorial on MM algorithms,” The AmericanStatistician , vol. 58, no. 1, pp. 30–37, 2004.[28] P. Stoica and Y. Selen, “Cyclic minimizers, majorization techniques,and the expectation-maximization algorithm: a refresher,”
IEEE SignalProcessing Magazine , vol. 21, no. 1, pp. 112–114, Jan 2004.[29] J. Li, D. Zheng, and P. Stoica, “Angle and waveform estimation viaRELAX,”
IEEE Transactions on Aerospace and Electronic Systems ,vol. 33, no. 3, pp. 1077–1087, July 1997.[30] P. T. Boufounos and R. G. Baraniuk, “1-bit compressive sensing,” in ,March 2008, pp. 16–21.[31] L. H. Nguyen, T. Tran, and T. Do, “Sparse models and sparse recoveryfor ultra-wideband SAR applications,”
IEEE Transactions on Aerospaceand Electronic Systems , vol. 50, no. 2, pp. 940–958, April 2014.[32] L. H. Nguyen and T. D. Tran, “Efficient and robust RFI extractionvia sparse recovery,”
IEEE Journal of Selected Topics in Applied EarthObservations and Remote Sensing , vol. 9, no. 6, pp. 2104–2117, June2016.[33] M. Ressler, L. Nguyen, F. Koenig, D. Wong, and G. Smith, “TheArmy Research Laboratory (ARL) synchronous impulse reconstruction(SIRE) forward-looking radar,” in
Unmanned Systems Technology IX ,G. R. Gerhart, D. W. Gage, and C. M. Shoemaker, Eds., vol. 6561,International Society for Optics and Photonics. SPIE, 2007, pp. 35 –46.[34] P. Stoica and Y. Selen, “Model-order selection: a review of informationcriterion rules,”
IEEE Signal Processing Magazine , vol. 21, no. 4, pp.36–47, July 2004.[35] P. Stoica, R. L. Moses, B. Friedlander, and T. Soderstrom, “Maximumlikelihood estimation of the parameters of multiple sinusoids from noisymeasurements,”