Analysis of Frequency Agile Radar via Compressed Sensing

# Analysis of Frequency Agile Radar via Compressed Sensing

Tianyao Huang, Yimin Liu,  Xingyu Xu, Yonina C. Eldar,  Xiqin Wang Partial results [Huang2015] of this work were presented at the IEEE China Summit and International Conference on Signal and Information Processing, Chengdu, China, Sep. 2015. The work of T. Huang, Y. Liu, X. Xu and X. Wang was supported by the National Natural Science Foundation of China under Grant 61801258 and 61571260. (Corresponding author: Yimin Liu.) T. Huang, Y. Liu, X. Wang and X. Xu are with the Department of Electronic Engineering, Tsinghua University, Beijing 100084, China (e-mail:{huangtianyao, yiminliu, wangxq_ee}@tsinghua.edu.cn, xy-xu15@mails.tsinghua.edu.cn). Yonina C. Eldar is with the Department of Electrical Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel (e-mail: yonina@ee.technion.ac.il).
###### Abstract

Frequency agile radar (FAR) is known to have excellent electronic counter-countermeasures (ECCM) performance and the potential to realize spectrum sharing in dense electromagnetic environments. Many compressed sensing (CS) based algorithms have been developed for joint range and Doppler estimation in FAR. This paper considers theoretical analysis of FAR via CS algorithms. In particular, we analyze the properties of the sensing matrix, which is a highly structured random matrix. We then derive bounds on the number of recoverable targets. Numerical simulations and field experiments validate the theoretical findings and demonstrate the effectiveness of CS approaches to FAR.

## I Introduction

Frequency agile radars (FARs) are pulse-based radars, in which the carrier frequencies are varied in a random/pseudo-random manner from pulse to pulse as illustrated in Fig. 1. Each transmission occupies a narrow band (). Pulse returns of different frequencies are processed coherently to synthesize a wider band (), which generates high range resolution (HRR) profiles.

Since the works [liu2000sp, Axelsson2007], frequency agility has received increasing attention [Liu2008, Huang2012, Huang2014, Zhang2011, Yang2013, Liu2014, Cohen2017a] in the radar community due to its multi-fold merits. First, frequency agility introduces good electronic counter-countermeasures (ECCM) performance, because the randomly varied frequencies of the pulses are difficult to track and predict. In addition, the flexibility of the narrow band transmission makes it easier to avoid and reject barrage jamming. Second, like stepped frequency radar, FAR can be used for two-dimensional (2D) imaging [Yang2013, Liu2014], while requiring only a narrow band receiver, which significantly lowers the hardware system cost. Third, in contrast to linearly stepped frequency radar, FAR decouples the range-Doppler parameters and produces a thumbtack ambiguity function [Liu2008]. It also mitigates aliasing artifacts in synthetic aperture radar (SAR) [Luminati2004] and extends the unambiguous Doppler window in inverse SAR (ISAR) imaging [Huang2012]. Finally, frequency agility can be utilized to exploit vacant spectral bands [Cohen2017a], and shows potential to increase spectrum efficiency and cope with spectrum sharing issues in a contested, congested and competitive electromagnetic environment.

We consider the problem of joint HRR profile and Doppler estimation of the target. When a traditional matched filter is used for joint range-Doppler estimation, sidelobe pedestal problems occur in FAR [Liu2008]. Therefore, weak targets could be masked by the sidelobe of dominant ones, which restricts the application of FAR in target detection and feature extraction [Liu2008]. By exploiting target sparsity, compressed sensing (CS) techniques [Eldar2012, Eldar2015] have been applied in order to alleviate the sidelobe pedestal problem. Liu et al. [Liu2008] propose the RV-IAP algorithm for joint range-Doppler estimation, which is based on the Orthogonal Matching Pursuit (OMP) method [Tropp2007]. Since then, many practical CS algorithms for FAR have been developed [Huang2012a, Huang2014].

This paper focuses on the theoretical analysis of CS methods for FAR in terms of reconstruction performance. Theoretical conditions that guarantee perfect recovery in general CS have been extensively studied. Randomness plays a key role in many theoretical results, and often leads to good empirical results [Eldar2012]. Near optimal conditions for Gaussian, sub-Gaussian, Bernoulli and random Fourier matrices have been derived [Baraniuk2008a, Krahmer2014] (and references therein). However, the measurement matrix in FAR differs from these random matrices, so that previous theoretical results are not directly applicable.

We start by studying the measurement matrix properties of a FAR system. This matrix is random, due to the randomness of the frequencies. We begin by deriving probability bounds on the spark and coherence of the measurement matrix, extending the results of [Huang2015]. Based on these bounds, we develop recovery guarantees for joint HRR profile and Doppler estimation using FAR systems. Theoretical results show that owing to the randomness of the carrier frequencies, with high probability, one can jointly obtain a HRR profile and Doppler of targets while transmitting narrow-band pulses. The number of recoverable targets is proved to be using minimization, or on the order of using minimization, where is the number of pulses.

We next perform simulations and field experiments to demonstrate the reconstruction performance of FAR using CS methods. We build an X-band FAR prototype with a synthetic bandwidth of 1 GHz, and test the recovery performance in a real environment. The results show that both the HRR profiles and Doppler of the observed target (a moving car) are reconstructed with and .

The rest of paper is organized as follows. Section II introduces the signal model and problem formulation. In Section III, a brief review of CS algorithms and their performance guarantees is provided. We derive conditions for joint range-Doppler recovery using FAR in Section IV. Numerical simulation and field experiment results are shown in Section V. Section VI concludes the paper.

Throughout the paper we use the following notation. The sets , , , refer to complex, real, integer, and natural numbers. Notation is used for the modulus, absolute value or cardinality for a complex, real valued number, or a set, respectively, and . For , (or ) is the largest (smallest) integer less (greater) than or equal to . Uppercase boldface letters denote matrices (e.g., ), and lowercase boldface letters denote vectors (e.g., ). The ,-th element of matrix is written as , and denotes the -th entry of a vector. Given a matrix , a number (or a set of integers, ), () denotes the -th column of (the sub-matrix consisting of the columns of indexed by ). As for a vector , denotes the sub-vector consisting of the elements of indexed by . The complex conjugate operator, transpose operator, and the complex conjugate-transpose operator are , , , respectively. We use , as the norm of an argument, and denotes the probability of an event. Operations and represent the expectation and variance of a random argument, respectively. The real and imaginary part of a complex valued argument are denoted by and , respectively.

## Ii Signal Model

In this section, we introduce FAR, following the presentation in [Huang2014]. A FAR system transmits monotone pulses, where the -th transmitted pulse is written as

 Tx(n,t):=rect(t−nTrTp)ej2πfn(t−nTr), (1)

, where and are the pulse repetition interval and pulse duration, respectively, , and rect represents the rectangular envelope of the pulse

 rect(x):={1, 0≤x≤1,0, otherwise. (2)

The frequency of the -th pulse is randomly varied as , where is the initial frequency, is the -th random frequency-modulation code, , and is the synthetic bandwidth. For a single pulse, the bandwidth () is narrow, , and the coarse range resolution (CRR) is , where is the speed of light. Synthesizing echoes of different frequencies refines the range resolution to . We denote the number of HRR bins inside a CRR bin as

 M:=⌈Tpc2⋅2Bc⌉=⌈TpB⌉∈N. (3)

Received echoes are assumed delays of the transmissions. We begin by assuming that there is a single ideal scatterer with scattering coefficient . The echo of the -th pulse can then be written as

 Rx(n,t):=βTx(n,t−2r(t)c), (4)

where denotes the range of the scatterer with respect to the radar at time instant . We assume that the scatterer is moving along the line of sight at a constant speed , so that . After down conversion, the echo becomes

 Rd(n,t):=Rx(n,t)⋅e−j2πfn(t−nTr)=βrect(t−2r(t)/c−nTrTp)ej2πfn(t−2r(t)/c−nTr)⋅e−j2πfn(t−nTr)=βrect(t−2r(t)/c−nTrTp)e−j2πfn2r(t)c. (5)

Echoes are sampled at the Nyquist rate of a single pulse, , so that each echo pulse is sampled once. Every sample corresponds to a CRR bin, and data from all CRR bins are processed in the same way. Returns of pulses from the same CRR bin are combined to a vector

 [Rd(0,t),Rd(1,Tr+t),…,Rd(N−1,(N−1)Tr+t)], (6)

and processed to generate HRR profiles and Doppler estimates. During the coherent processing interval (CPI), i.e. , we assume that the scatterer does not cross a CRR bin, which means that

 vNTr

Without loss of generality, suppose that the -th CRR bin contains the scatterer, . The corresponding sampling instant for the -th pulse is . Substituting into (5), the sampled echoes are given by

 Rd(n,nTr+l/fs)=βe−j4π(fc+dnB)(r(0)+v(nTr+l/fs))/c≈βe−j4πfcr(0)+vl/fsce−j4πdnBr(0)/ce−j4πfcvTrnζn/c, (8)

where the approximation holds if the term , which requires . Here . Generally, different carrier frequencies imply different Doppler shifts, unless the relative bandwidth is negligible, i.e. . However, in a (synthetic) wideband radar, this approximation does not usually hold, and could give rise to estimation performance deterioration in practice if applied. In the simulations and field experiments in Section V, the signal processing algorithms do not adopt this assumption. However, in the mathematical analysis in Section IV, we assume for theoretical convenience. The impact of the relative bandwidth will be discussed in the simulations.

For brevity, we omit the notation , and write . We further introduce notations

 ⎧⎪ ⎪⎨⎪ ⎪⎩~γ:=βe−j4πfcr(0)+vl/fsc,~p:=−4πBr(0)/(Mc),~q:=−4πfcvTr/c. (9)

With these definitions (8) becomes

 Rd(n)≈~γej~pMdn+j~qnζn. (10)

After the unknowns , and are estimated, the absolute intensity , HRR range and velocity are inferred as , and , respectively.

When there are scatterers occurring inside the CRR cell, radar returns are modeled as a combination of returns from all scatterers,

 Rd(n)=K−1∑k=0~γkej~pkMdn+j~qknζn, (11)

where , and in (10) are replaced with , and for the -th scatterer, respectively.

To avoid grating lobes in the HRR profiles (which are also called ghost images in the literature) [Liu2009, Liu2014a], the frequency codes are required to satisfy , . When codes are discrete, we denote by the set of available frequency codes and by the number of codes. The codes are often uniformly spaced, e.g. . It is required that and a typical choice is . When codes belong to a continuous set (for some of our theoretical results), the requirement is usually easy to satisfy with . We assume that in both discrete and continuous cases, the codes are identically, independently, and uniformly distributed.

### Ii-B Signal Model in Matrix Form

We can rewrite (11) in matrix form as

 y=Φx, (12)

where the measurement vector has entries . The vector corresponds to the scattering intensities . The pair defines the key target parameters, range and Doppler, and belongs to a continuous 2D domain. The resolutions for and are and , respectively. Consider the unambiguous continuous region , and discretize and at the Nyquist rates, and , respectively. Thus, one obtains and , , . Denote the sets containing HRR grids and Doppler grids as and , respectively, and assume that the targets are located precisely on the grid. Define the matrix with entries

 [X]m,n={~γk, if ∃k, (~pk,~qk)=(pm,qn),0, otherwise, (13)

representing the 2D scattering coefficients in the range-Doppler domain, and . We vectorize to obtain with entries .

To introduce the measurement matrix , we define the matrices and , corresponding to HRR range and Doppler parameters, respectively, with entries

 [R]n,m:=ejpmMdn, (14)
 [D]n,l:=ejqlnζn, (15)

, and . If , then is a Fourier matrix. Define , where denotes the Khatri-Rao product. Then the elements of are given by

 [Φ]n,l+mN:=[R]n,m[D]n,l=ejpmMdn+jqlnζn, (16)

and . When echoes are corrupted by additive noise , (12) becomes

 y=Φx+w. (17)

The sensing matrix in (17) has more columns than rows, , which shows that joint range and Doppler estimation in FAR is naturally an under-determined problem. When is -sparse, which means there are non-zeroes in , and , CS algorithms can be used to solve (17). The targets’ parameters can then be recovered from the support set of .

### Ii-C Discussion on the Signal Model

Note that when there is only one scatterer observed, the matched filter that maximizes the signal to noise ratio (SNR) works well in FAR. However, when there are multiple scatterers, sidelobe pedestal problem occurs and weak targets can be masked by the dominant targets’ sidelobe. The matched filter estimates the scattering intensities by

 ^x:=ΦHy=ΦHΦx+ΦHw. (18)

In such an under-determined model, , spurious responses emerge in even if there is no noise, i.e. . These spurious responses are the sidelobe pedestal.

To better interpret the sidelobe pedestal problem, we compare the signal model of FAR with that of an instantaneous wideband radar (IWR). In such a hypothetical radar, we assume that the radar transmits/receives all of its sub-bands (with as the set of frequency codes, ), and processes the echoes individually for each band. In FAR, the same set is also applied with . In analogy to (11), the return of the -th frequency in the -th pulse can be written as

 RIWR(m,n):=K−1∑k=0~γkej~pkm+j~qknηm, (19)

where , , and . For notational brevity and simplicity, we assume and for IWR and FAR, respectively. In this case, (19) can be rewritten in matrix form as

 Z=FXDT, (20)

where has entries , and is a Fourier matrix with entries , , . Equivalently,

 z=(F⊗D)x, (21)

where and denotes the Kronecker product. The sensing matrix in the IWR is orthogonal, i.e. , and the sidelobe pedestal problem vanishes.

The measurements in FAR can be regarded as sampling111Since an instantaneous narrowband waveform is used in FAR, it naturally enjoys low data rate in comparison with IWR. However, this paper does not aim at minimizing the data rate. It may have the potential to further reduce the data rate by combining frequency agility with approaches like sub-Nyquist sampling in the fast-time domain [Bar-Ilan2014, Cohen2017a, Xi2014], omitting some pulses or frequency bands [Hu2011, Cohen2017]. of the IWR measurements, i.e.

 [y]n=[Z]Mdn,n, n=0,1,…,N−1. (22)

Only one sub-band data is acquired for each pulse. Therefore the sensing matrix of FAR consists of partial rows of that of the IWR, i.e.,

 (ΦT)n=(ΨT)n+MdnN, n=0,1,…,N−1, (23)

and becomes an under-determined matrix. This interpretation suggests that the sidelobe pedestal of FAR results from the information loss in the frequency domain. The spectral incompleteness leads to an under-determined problem (12). In Section IV, we prove that, owing to the randomness of the frequencies, the scatterers can still be correctly reconstructed via CS methods with high probability.

## Iii Review of Compressed Sensing

In Section IV, we prove that using CS methods, FAR can provably recover the HRR profiles and Doppler. Before deriving the results, we review some basic notions of CS [Eldar2015].

Consider an under-determined linear regression problem, e.g. (12), where is sparse. The sparsest solution can be obtained via

 minx∥x∥0, s.t. y=Ax, (P0)

where denotes “norm” of a vector, i.e. the number of non-zeroes. This solution is the true vector, when the sensing matrix has the spark property.

###### Definition 1 (Spark, [Eldar2015]).

Given a matrix , Spark() is the smallest possible number such that there exists a subgroup of columns from that are linearly dependent.

Unique recovery of can be ensured if the following condition is satisfied.

###### Theorem 2.

The equation is uniquely solved by () if and only if .

The above theorem provides a fundamental limit on the maximum sparsity that leads to successful recovery. In general, optimization is NP-hard. A widely used alternative is basis pursuit, which solves the problem

 minx∥x∥1,s.t. y=Ax. (P1)

In noisy cases, variants like basis pursuit denoising, LASSO and Dantzig selector can be applied. Many greedy methods have also been suggested to approximate ().

Sufficient conditions that guarantee uniqueness using these methods are extensively studied. Bounds on the mutual incoherence property (MIP) and restricted isometry property (RIP) are widely applied conditions to ensure sparse recovery. In this paper, we rely on the MIP. A matrix has MIP if its coherence is small, where coherence is defined as the maximum correlation between two columns, i.e.

 μ(A):=maxl≠k∣∣AHlAk∣∣∥Al∥2∥Ak∥2. (24)
###### Theorem 3 ([Fuchs2004]).

If a matrix has coherence , then for any of sparsity , is the unique solution to ().

The condition in Theorem 3 ensures recovery in the presence of noise and also recovery using a variety of computationally efficient methods[Eldar2012].

## Iv Sensing Matrix Properties of FAR

In this section, we analyze the spark and MIP properties of the FAR’s sensing matrix. These results are then used together with Theorems 2 and 3 to establish performance guarantees for FAR. In the following derivations, we assume that .

### Iv-a Spark Property

The following theorem proves that the sensing matrix of FAR almost surely has the spark property.

###### Theorem 4.

Consider defined in (16) with drawn independently from a uniform continuous distribution over , . Then, with probability , Spark.

###### Proof.

See Appendix A. ∎

Since has rows, there must be a linearly dependent submatrix with columns. Owing to the randomness of the carrier frequencies, Theorem 4 shows that a sub-matrix built from any columns of is of full rank almost surely. The result is based on the assumption of a continuous distribution on . An immediate consequence of Theorems 2 and 4 is the following corollary.

###### Corollary 5.

Consider a FAR whose frequency modulation codes are drawn independently from a uniform continuous distribution over , . Then, with probability 1, scatterers can be exactly recovered by (), where is the number of pulses.

### Iv-B Mutual Incoherence Property

To obtain performance guarantees using minimization or greedy CS methods under noiseless/noisy environments, we derive the MIP for FAR. We start by analyzing the asymptotic statistics of the FAR’s sensing matrix. Then invoking Theorem 3, we obtain the maximum number of scatterers that FAR guarantees to exactly reconstruct with high probability.

Assume in this subsection that , , and recall that the parameters and are on a grid, i.e. and . First, consider the Gram matrix , which links to the coherence. Define as the modulus matrix of the Gram matrix, i.e.

 [G]k,l=∣∣∣[ΦHΦ]k,l∣∣∣,k,l=0,1,…,NM−1. (25)

We then have the following results, some of which are partially inspired by [Rossi2014].

###### Lemma 6.

The rows of the modulus matrix are permutations of elements in its first row.

###### Proof.

Denote by and , , two columns in , corresponding to and , respectively, , . Then

 ΦHl1Φl2=N−1∑n=0e−jpm1Mdn−jqk1n⋅ejpm2Mdn+jqk2n=N−1∑n=0e−j(pm1−pm2)Mdn−j(qk1−qk2)n=N−1∑n=0e−j2π(m1−m2)dn−j2π(k1−k2)Nn. (26)

Clearly (26) depends only on the difference between grid points, i.e. and , and is independent of the particular indices and . In addition, . Therefore, for any element in , one can find an element in the first row of with the same value. ∎

Consider now the -th element in the -th row of , , which corresponds to the -th column of . Note that each column of relies on a specific parameter pair . For notational simplicity, we drop the subscripts of related to , and define

 χl:=1NΦH0Φl=1NN−1∑n=0ejpMdn+jqn, l=1,…,NM−1. (27)

We now analyze the mutual coherence, . Since is random, is also random unless , in which case reduces to a constant for . This constant does do not affect the value of and is thus ignored. Define a set excluding these constants as

 (28)

Then is a random variable, , and has the following statistical characteristics.

###### Lemma 7.

As , the real and imaginary parts of , and , , have a joint Gaussian distribution,

 [Re(χl)Im(χl)]∼N([00],[12N0012N]), (29)

except in the special case that the corresponding parameters . In this setting, the joint Gaussian distribution becomes

 [Re(χl)Im(χl)]∼N([00],[1N000]). (30)
###### Proof.

See Appendix B. ∎

The special case corresponds to a specific grid point of the plane. Considering the generic case and this special case separately leads to the following conclusions.

###### Corollary 8.

When and , with ,

 P(|χl|>ϵ)≤e−Nϵ2/2. (31)
###### Proof.

Lemma 7 proves that when and do not equal simultaneously, the real and imaginary parts of asymptotically obey independently. Therefore, the magnitude obeys a Rayleigh distribution with probability density function , , and cumulative distribution function , . Thus, , which yields

 P(|χl|>ϵ)=e−Nϵ2≤e−Nϵ2/2. (32)

In the special case that , the real part of asymptotically obeys independently and the imaginary part vanishes. Then elementary estimates of the Gaussian error function yield

 P(|χl|>ϵ)≤√2πNϵ2e−Nϵ2/2, (33)

which is less than if , i.e. . ∎

###### Lemma 9.

The maximum , , satisfies the following as ,

 P(μ>ϵ)≤(MN−N)e−Nϵ2/2. (34)
###### Proof.

For fixed , we have as . According to the union bound

 P(μ>ϵ)≤∑l∈ΞP(|χl|>ϵ)≤(MN−N)e−Nϵ2/2, (35)

since there are indices in . ∎

We next derive a condition for FAR to meet the requirement of Theorem 3 with high probability.

###### Theorem 10.

The coherence of , defined in (16), obeys with a probability higher than , when

 K≤12√2√Nlog(MN−N)−logδ+12. (36)
###### Proof.

Let . From (36), we have that

 Nϵ2≥2(log(MN−N)−logδ). (37)

Assume that and . Then

 Nϵ2≥2log2>2π. (38)

Using (34) and (37), we finally obtain

 P(μ≤ϵ)>1−(MN−N)e−Nϵ2/2≥1−δ, (39)

completing the proof. ∎

Theorem 10 shows that the sensing matrix of FAR has the MIP; thus, according to Theorem 3, HRR range-Doppler reconstruction is guaranteed if the number of targets satisfies , where in the numerator represents the number of measurements, and in the denominator links to the number of grid points. The work [Rossi2014] for direction of arrival estimation using random array MIMO radar also proposes a bound based on MIP, which guarantees recovery of a number of targets on the order of , where and are the number of measurements and grid points, respectively. In [Rossi2014], a bound for non-uniform recovery based on RIPless theory is also provided, and is relaxed to . However, RIPless is not directly applicable to FAR, because for each row of , i.e. , is rank deficient and the isotropy property does not hold. Dorsch and Rauhut [Dorsch2017] analyze the joint angle-delay-Doppler recovery performance using MIMO based on RIP. They assume a periodic random probing signal with independent samples. In this case the recoverable number of scatterers is on the order of . Though the RIP leads to a tighter bound than MIP, the RIP of the FAR system matrix is still an open question.

Sufficient conditions that guarantee uniform recovery are usually pessimistic. It is well known that CS algorithms often outperform the theoretical uniform recovery guarantees. In the next section, we evaluate the practical performance of FAR using CS methods.

## V Simulation and Experimental Results

In this section, simulations and field experiments are executed to demonstrate the properties of the sensing matrix of FAR and the effectiveness of CS algorithms to reconstruct the targets’ HRR range and Doppler.

### V-a Spark Property

First, the spark property of the sensing matrix is discussed. We construct a sub-matrix of , where the set and , and calculate the minimum singular value of . We check whether it is equal to zero (rank deficient). Concretely, we set and , which are small to make it possible to enumerate all the sub-matrices of . We record the minimum singular value (normalized by ) of each sub-matrix; thus, we obtain results, among which the minimum is denoted as . The frequency codes are distributed uniformly on a continuous set. We further assume the relative bandwidth satisfies . We perform 2000 Monte-Carlo trials. The histograms of and are depicted in Fig. 2 and Fig. 3, respectively. The minimum of and is . The results indicate that a continuous distribution of codes results in good properties of the sensing matrix. For comparison, we also perform simulations with codes distributed on the discrete set , and count the number of minimum singular values , i.e. , less than , which leads to Pr. Therefore, a continuous distribution leads to better spark performance than a discrete distribution.

### V-B Mip

Next, we consider the MIP of the sensing matrix. The parameters are set to and . The frequency codes are uniformly distributed over the discrete set . The relative bandwidths are set to , where means that the assumption holds. We also simulate the continuous case with and . Curves are obtained with Monte-Carlo trials. For each trial, we calculate the mutual coherence of the sensing matrix and depict the corresponding cumulative distribution function. The theoretical bound in (34) is also displayed. The results are shown in Fig. 4. It can be seen that the theoretical upper bound (34) is tight under the assumption that the relative bandwidth is negligible (thus holds). When the relative bandwidth is large, the actual mutual coherence could exceed the predicted one, e.g., in the case that . However, the curve of is under that of , which indicates that a larger relative bandwidth does not necessarily result in worse mutual incoherence.

### V-C Recovery Performance in Noiseless Cases

In Fig. 5, we consider the recovery performance in noiseless cases. In particular, we plot the probability of exact recovery using the basis pursuit algorithm () and matched filter (18), where exact recovery means the support set of the unknown vector is exactly estimated. In the simulations, the pulse number is , the number of frequencies is , and the amplitudes of scattering coefficients are all set to 1. The number of scatterers, , is varied. The initial carrier frequency is GHz and the bandwidth is MHz. For each point on the curve, we perform 200 Monte-Carlo trials, where the frequency codes are randomly drawn obeying , the support set of is random, and the phases of non-zeros in are i.i.d . We solve () using CVX [cvx, gb08]. In both methods, we assume the number of scatterers, , is known, and the support set is obtained as the indices of the largest magnitudes in . The magnitudes are also compared with a threshold . Those not exceeding the threshold are removed from the support set. From Fig. 5, it is seen that CS dramatically outperforms the traditional matched filter. When , the support set recovery probabilities using matched filter drops significantly, because some of the scatterers are masked by the sidelobes. When basis pursuit is applied, the region leading to exact support set recovery in FAR is fairly broad. However, the theoretical bound (36) is 1.5 with , and is quite pessimistic.

### V-D Recovery Performance in Noisy Cases

We next consider noisy cases, and choose the probability of successful recovery to evaluate the performance of different CS algorithms. A successful recovery is defined as exact recovery of the support set. In the simulations, GHz, MHz, , and the number of scatterers . The scatterers have identical amplitudes of 1 with random phases. The noise in (17) is assumed Gaussian white noise with a covariance matrix , and varies from -15 dB to 15 dB. Subspace pursuit [Dai2009] and Lasso are compared. In subspace pursuit, the number of scatterers is assumed known a priori. The Lasso algorithm solves

 minx12∥y−Φx∥22+λ∥x∥1 (40)

with , and is implemented with CVX. When the magnitude of the estimate is larger than , the corresponding index is put into the estimated support set. The results are shown in Fig. 6 with 200 Monte-Carlo trials. Both algorithms have high successful recovery probabilities for an FAR in high SNR, dB. Subspace pursuit outperforms Lasso with the genie-aided information on the cardinal number of the support set. We also note that the selection of the parameters and has a significant impact on the performance of Lasso.

### V-E Field Experiments

Next, we show field experiments from a true FAR prototype. We use separated antennas for transmitting and receiving, respectively, so that returns from short range objects are not eclipsed by the transmission. The radar works at an initial frequency of GHz. Frequencies are varied pulse by pulse. There are frequencies with a minimum gap of MHz, which results in a synthetic bandwidth of MHz. In each pulse, the carrier frequency is randomly chosen. Specifically, and . The HRR is m. We set ms and the equivalent pulse duration ns. The CRR is m, and the number of HRR bins in a CRR bin is , which satisfies the condition to eliminate ghost images [Liu2009, Liu2014a]. The number of pulses is . The moving target does not cross a CRR bin during the CPI, which requires m/s. For a slowly moving car, the velocity is lower than 10 m/s.

Field experiments are executed to evaluate performance. The target is a household car (see Fig. 7) with two corner reflectors and four small metal spheres upon its roof to enhance the SNR. When the radar is operated, the car moves in front of the radar at a nearly constant speed along the road, surrounded by static objects including a big stone, roadside trees and iron barriers. Returns from all scatterers located in , i.e. km, are collected. There are CRR bins. We perform static clutter canceling [Axelsson2006] over all CRR bins. Then in these CRR bins, we find the CRR bin which has the maximum amplitude of radar echoes, and infer that the car is located in that bin. With data in that CRR bin, we perform range-Doppler processing of the target. We apply the matched filter (18) and the OMP algorithm to jointly estimate the HRR profile and Doppler of the target. In OMP, the algorithm iterates 50 times, i.e., assuming . The results are shown in Fig. 8 and Fig. 9, respectively. To better demonstrate the sidelobe, we project the three dimensional images in Fig. 8 and Fig. 9 onto amplitude-range dimensions; see Fig. 10. In Fig. 10, only the maximum amplitudes in using matched filter are found and shown. In both methods, the velocity estimation is 6.9 m/s, and the span of the car is around 1.9 m. Comparing the recovery performance, we see that the matched filter suffers from sidelobe pedestal while OMP demonstrates a clearer reconstruction.

## Vi Conclusion

In this paper, sparse recovery for a frequency agile radar with random frequency codes is studied. We analyzed the spark and mutual incoherence properties of the radar sensing matrix, which guarantee reliable reconstruction of the targets. Using minimization, FAR exactly recovers scatterers in noiseless cases almost surely, where is the number of pulses. When we apply minimization or greedy CS methods and there is noise, the number of scatterers that is guaranteed to be reliably reconstructed by FAR is on the order of , where is the number of HRR bins in a CRR bin. Numerical simulations and field experiments were executed to validate the theoretical results and also demonstrate the practical recovery performance of FAR.

## Appendix A Proof of Theorem 4

To avoid confusion, in this appendix we will index the row of a matrix by and the column by . We need to prove that any columns of are almost surely linearly independent. Following the form in (16), we first fix some columns , with running through . These columns constitute a new matrix whose elements are . We need to show that is almost surely invertible. For brevity, we set

 zξ:=exp(j2πdξ),cξ,η:=exp(j2πlηξ/N). (41)

With the notation above, we may write . The proof follows from the following three lemmas. The first one, being more abstract, is a strengthened version of the well-known fact that the set of zeros of a nonzero polynomial has measure zero.

###### Lemma 11.

Let be a nonzero complex polynomial in -variables. Denote by the set of zeros of . Consider the -torus with its obvious embedding . Let be the Haar measure on . We have

 σn(N∩Tn)=0. (42)
###### Proof.

The left hand side of the above equality is well-defined, since is closed and, consequently, is a closed set in . Note that coincides with the product measure . This enables us to prove the result by induction via Fubini’s theorem.

For the set is discrete, thus the proposition is trivial. Suppose for the proposition is true. Let be the projection onto the first component. We now make the natural identification . With the result for in hand, we can find a -negligible set (a set of measure zero) such that for any , the polynomial (viewed as a polynomial in ) has at least one nonzero coefficient, i.e. is a nonzero polynomial in . Then by the result for , for , the measure of the slice

 σ1(π(N∩Tk+1∩{(z2,⋯,zk+1)=z(k)}))=0. (43)

We write this fact as . Bearing in mind that , we apply Fubini’s theorem to obtain

 σk+1(N∩Tk+1)=∫Tk+11Ndσk+1=∫Tk∖Odσk∫S11N⋅μ(dz1)+∫Odσk∫S11N⋅σ1(dz1)=∫Tk∖Oσ1(π(N|z(k)))σk(dz(k))+0=0, (44)

which completes the proof. ∎

###### Lemma 12.

Let be independent random variables with continuous distributions for . Fix some complex numbers , positive real constants and nonnegative integers , . Let . Then the following two statements are equivalent:

1. The random matrix with elements

 [A]ξ,η=cξ,ηzmηξ (45)

is almost surely invertible.

2. There exists a vector such that the deterministic matrix with elements

 [A(w)]ξ,η=cξ,ηwmηξ (46)

is invertible.

###### Proof.

(i)(ii) is obvious. We prove that (ii)(i). Note that is a polynomial in variables . By (ii), this polynomial is nonzero, since