Random filtering structure-based compressive sensing radar

Recently with an emerging theory of ‘compressive sensing’ (CS), a radically new concept of compressive sensing radar (CSR) has been proposed in which the time-frequency plane is discretized into a grid. Random filtering is an interesting technique for efficiently acquiring signals in CS theory and can be seen as a linear time-invariant filter followed by decimation. In this paper, random filtering structure-based CSR system is investigated. Note that the sparse representation and sensing matrices are required to be as incoherent as possible; the methods for optimizing the transmit waveform and the FIR filter in the sensing matrix separately and simultaneously are presented to decrease the coherence between different target responses. Simulation results show that our optimized results lead to smaller coherence, with higher sparsity and better recovery accuracy observed in the CSR system than the nonoptimized transmit waveform and sensing matrix.


Introduction
Compared with the whole scene observed by radar systems, the target scene is sparse therein in the majority of cases. Classical radars do not take advantages of this sparsity and lead to complicated and expensive radar receiver consisting of high-rate analog-to-digital (AD) converters, large memories, and fast computing systems.
Recently with an emerging theory of 'compressive sensing' (CS) [1][2][3][4], a radically new concept of compressive sensing radar (CSR) has been proposed [5]. According to CS theory, CSR can recover the target scene from far fewer samples or measurements than traditional methods. To make this possible, CSR relies on two principles: sparsity, which restricts the number of targets of interest, and incoherence, which says the dissimilarity between targets of interest. Obvious characteristics of the CSR system can be summarized as follows [5,6]: • Eliminating the need for the pulse compression matched filter at the receiver • Reducing the required receiver AD conversion bandwidth so that it need operate only at the low 'information rate' rather than at the high Nyquist rate • Providing the potential to achieve higher resolution between targets than traditional radars whose resolution is limited by the uncertainty principles Two different tasks of CSR have been investigated by only a few papers. The first radar task is to detect and estimate targets in distinct range, Doppler and angle cells [6,7]. The second is imaging, including range profiling, synthetic aperture radar (SAR) and inverse synthetic aperture radar (ISAR) [8][9][10]. In both cases, CSR can work in the situation of sparse targets/scene. CSR was demonstrated to be capable of successfully working with an AD converter operating at a sampling frequency lower than the Nyquist rate. An exact recovery of target scene can be implemented with four times undersampling for CSR SAR imaging [5]. CSR was considered to transmit a sufficiently incoherent pulse and reconstruct the sparse target scene by the greedy algorithm. Better resolution in the timefrequency plane over traditional radar can be provided [6]. In [8][9][10], CS technique was applied to range profiling, azimuth domain focusing, and (ω, k) domain focusing in SAR imaging. http://asp.eurasipjournals.com/content/2014/1/94 CSR waveform design was also investigated by Chen [11] and Subotic [7]. To effectively reconstruct the target scene, it is required that the correlations between target responses must be small. A multiple-input multipleoutput (MIMO) radar waveform design method has been proposed based on simulated annealing (SA) algorithm [11]. For a distributed radar system, waveform and position impacts have been examined by considerations of sparsity of the target scene and the restricted isometry property (RIP) [7], [12].
In [5], Baraniuk has suggested an interesting choice for the sensing matrix of CSR system. He pointed out that in majority of cases we can use a causal, quasi-Toeplitz matrix where each row is a right shift of the row immediately above it. The sensing matrix based on random filtering is plotted in Figure 1. The measurement process in CSR can be seen as a linear time-invariant filter followed by decimation. When we choose a pseudo-noise (PN) sequence as the initial row, this approach can be named as random filtering. Obvious benefits of random filtering structure-based CSR can be summarized as [13] • The sensing matrix is stored and applied efficiently.
• Fast fourier transformation (FFT) can be used to replace convolution for long filters. • It is easily implementable in software or hardware.
The constitutions of the two kinds of random filtering structure-based compressive sensing radar is illustrated in Figure 2. As we can see, the random filter can be implemented both in time and frequency domain. In Figure 2, we also note that two Toeplitz matrices have been used and concatenated. The first appears due to the convolution of the transmitted waveform with the sparse target scene, spreading the energy of each target in the scene. The second is caused by the convolution of a random filter with the received signal. Actually, the input for the second Toeplitz matrix is no longer sparse in the canonical basis because of the concatenation of the two matrices.
CS theory suggests that the sensing matrix and sparse representation matrix be as incoherent (orthogonal) as possible. A measure of coherence between and can be expressed as where , denotes the inner product, i is the ith row of , and j is the jth columns of . ρ( , ) plays an important role in the successful recovery of basis pursuit (BP) and orthogonal matching pursuit (OMP) algorithms [14]. Low coherence between and means small ρ( , ). We note that a key notion of the RIP in CS theory also requires the orthogonality between and [15]. The RIP can be expressed as [16] ∀|T| ≤ S : where ∀|T| ≤ S mean for any sparsity T, which is less than S, | · | is the cardinality operator, · 2 represents the l 2 norm being equivalent to the square root of the sum of squares of all the elements, D = is the equivalent dictionary, D T is a subset extracted from D, θ T are the coefficients corresponding to the T selected columns, and 0 < δ S < 1 is the S-restricted isometry constant (RIC). If the RIP holds, any subset of columns of D are nearly orthogonal and the incoherence between and is ensured. However, the RIP is difficult for us to verify [17]. Therefore, some matrices have been proved to be incoherent enough with any fixed sparsifying basis with overwhelming probability, such as Gaussians or ±1 random matrices [16].
Elad has proposed an alternative framework towards the incoherence required by CS [18]. This alternative, which has been shown to be computationally more efficient and produce significantly better results, can be described as [19] where | · | is the absolute value, μ(D) is often called the mutual coherence of the matrix D, and d i = i denotes the ith column of D. The mutual coherence is known to be a sub-optimal metric to quantify CS matrices as compared to RIC. Notably, μ ≥ 1/ √ M for a M × N Gaussian matrix. Thus, using the mutual coherence metric, we have a sub-optimal quadratic scaling of M with the sparsity S. In comparison, a linear scaling of M with S is achieved with the RIC.
If D is designed such that μ(D) is as small as possible, the orthogonality between and can be guaranteed and successful recovery will be implemented in CS process. Here, the transposition (·) T was applied for image processing in real number domain. http://asp.eurasipjournals.com/content/2014/1/94 In this paper, we are concerned with incoherence between the sensing and sparse representation matrices in random filtering structure-based CSR. With the thought that the sensing matrix and transmit waveform in CSR can be changed in mind, we will investigate the problem of how to design the sensing matrix and transmit waveform to guarantee incoherence in the CSR system.
A similar thought has appeared in image processing and can be traced back to Elad's work [18]. Elad first attempted to decrease the average mutual coherence by optimizing the sensing matrix. His work showed that designing a sensing matrix is a better choice than a random matrix, and it indeed leads to better CS performance. Abolghasemi proposed a gradient descent method to optimize the sensing matrix [20]. Duarte-Carvajalino extended Elad's work and proposed to optimize the sparse representation and sensing matrices simultaneously [15]. Due to more freedom degrees introduced in CS, this new CS framework can offer better performance than only optimizing the sensing matrix. Overall, the results of these methods show enhancement in terms of both reconstruction accuracy and the maximum allowable sparsity CS can recover.
The remainder of this paper is organized as follows. First, we study the theory of random filtering structurebased compressive sensing radar in Section 2. Then, we introduce our proposed algorithms to design the transmit waveform and sensing matrix in Section 3. In Section 4, we present detailed experimental results demonstrating the superiority of our framework. Finally, concluding remarks and directions for future research are presented.

Sparse representation dictionary
Consider a target scene in time-frequency plane is discretized into an L × M grid, and we define time and frequency shift matrices as where L and M denote the numbers of range and Doppler bin CSR measures, N is the length of transmit waveform vector, ω M = e √ −12π/M is the Mth root of unity. The (l, m)th basis element in time-frequency plane can be defined as Assuming the transmit signal x of length N, the receiver will receive a signal h l,m = p l,m · x of length L to observe L range bins. Note that there are necessarily LM grid points in timefrequency plane; we concatenate these received signals and obtain the received signal basis dictionary where is the product we define as A and B are matrices of the same size m × n and s is a vector of size n × 1. The sparse representation dictionary contains all the possible signal reflected from the target in any grid of time-frequency plane.

Random filtering measurement
In CS, the sensing matrix measures and encodes P < L linear projections of the signal. By random filtering measurement, this process can be seen as the convolution of the received signal and the FIR filter f of length B, which approximates the analog filtering in the digital domain. To take P measurements of the signal, downsampling of the FIR filter output is then carried out. This process can be represented by a matrix , where is a P × L matrix. This matrix is banded and quasi-Toeplitz: each row has B nonzero elements, and each row of is is a shifted copy of the first row. http://asp.eurasipjournals.com/content/2014/1/94 By the fixed filter f which is defined on the region [0, B − 1], the (p, l)th entry of the sensing matrix can be obtained by where K = L/P is the decimation factor and · is the floor function. The pth row of can be expressed as

Scene recovery
The received signal of the target scene can be expressed as where θ = [θ 1,1 , θ 1,2 , . . . , θ L,M ] T is the scattering coefficient vector of the scene and n is the additive Gaussian white noise vector of length L. The measured signal is written as Then, the target scene can be recovered by reconstruction methods including using algorithms such as orthogonal matching pursuit [21,22] and basis pursuit [23,24]. The latter program is solving the following convex problem: where · 1 denotes the l 1 norm of a vector or matrix which is equal to the sum of absolute value of all the elements and > 0 takes into account the possibility of noise in the linear measurements and of nonexact sparsity. Regularized orthogonal matching pursuit (ROMP) has been proposed to take advantage of OMP and BP algorithms [25]. Equation 3 defines the maximum absolute value of normalized inner product between all columns in the equivalent dictionary D. Suppose the sparsity of the target scene θ 0 satisfies the following inequality

The FIR filter and transmit waveform design for compressive sensing radar
where · 0 denotes the l 0 norm counting the number of nonzeros in a vector or matrix. θ is necessarily the sparsest solution (min θ 0 ) such that y = Dθ.
A fast greedy algorithm such as OMP is guaranteed to succeed in finding the correct solution in the presence of noise n. The root mean square error (RMSE) of the solution θ obeys [14] θ − θ 2 where = n 2 and δ ≥ = y − Dθ 2 . The mutual coherence μ(D) that affects both the recoverable sparsity of target scene and the recovery accuracy is demonstrated in (12) and (13). We note that in the equivalent dictionary D, the FIR filter f in the sensing matrix and the transmit waveform x are variables in the CSR system. Therefore, the FIR filter f and transmit waveform x design problem in the CSR system can be described as where (·) H denotes the conjugate transposition. Here, we replace the transposition (·) T in Equation 3 by the conjugate transposition (·) H to process the columns of the complex equivalent dictionary D in the CSR system.

The transmit waveform optimization
With the fixed FIR filter f, the transmit waveform x optimization problem can be given by where We use the square of μ i,j and transform this fraction in (15) into a weighted summation. This expression can be obtained by where λ is the weighted coefficient, This transformation is valid by taking the log in (15) (due to the monotonicity of the log function). The optimization problem in (15) can be transformed to This minimax problem is difficult for us to solve, and we replace by minimum summation of γ i,j and can be described as This transformation is also valid because it is noted that for any vector x ∈ R N , we have 1/ √ N x 2 ≤ x ∞ ≤ x 2 . Thus, the optimization problem is simplified to the minimization problem defined in (18), in which the transmit waveform x has constant energy constraint When the transmit vector x is equal to the eigenvector of corresponding to the smallest eigenvalue, the minimization of x H x in (18) is achieved subject to the energy constraint of x H x = 1. However, the matrix which depends on x lead to indirect solution. Therefore, an iterative procedure must be applied. The specific steps involved in this iterative procedure are described below: • Step A1: Set the x with random generated values or use some existing sequence (i.e., Frank sequence or Golomb sequence), k = 0. • Step A2: Compute the matrix k+1 in terms of x k . • Step A3: Find the smallest eigenvalue and the corresponding normalized eigenvector v k+1 of the matrix k+1 . • Step A4: Repeat the above steps until the convergence criteria is satisfied, e.g., x k − x k+1 2 < 1 , where x k+1 = v k+1 is the waveform obtained at the k th iteration and 1 is a predefined threshold.

The FIR filter optimization
We define the complex Gram matrix G = D H D whose entry at the ith row and jth column is g i,j . Unlike the Gram matrix definition in [18], we do not compute G using the matrix D after normalizing each of its columns. In practice, the small absolute off-diagonal elements in G is desired. In the ideal case, minimum possible coherence occurs when g i,j = 0, i = j, and we have where G is a LM×LM diagonal matrix with self-coherence of each columns. We might be able to design the sensing matrix making the Gram matrix G close to G as possible. This process can be written as follows: where · F is the Frobenius norm for a matrix. By (10) and (11), the Gram matrix can be rewritten as According to (20), the sensing matrix optimization problem with the fixed sparse representation matrix can be described as With the successful case in [26], simpler criteria can be given to replace (22) and written as where U is a P×N semiunitary matrix (i.e., U H U = I), g = Diag( √ g 11 , √ g 22 , . . . , √ g NN ), Diag(·) denotes the diagonal matrix with diagonal elements as indicated.
Considering the sensing matrix also has energy constraint, the above criteria can be rewritten as and U are both unknown variables in (24). Our strategy to solve this minimization problem is calculating one variable while the other is fixed and iterating this process until convergence appears.
First, we consider the minimization problem with the known sensing matrix in (24) has the following solution [27] where U 1 and U 2 can be obtained by the following singular value decomposition (SVD) expression here U 1 is a P × P unitary matrix, U 2 is an P × LM semiunitary matrix, and is a P × P diagonal matrix. Then, we try to find f which determines the sensing matrix with given U. We decompose the sensing matrix by rows and have After transforming the above expression, we obtain The FIR filter f can be obtained by the least square (LS) estimator The FIR filter f optimization method can be summarized as follows: • Step B1: Generate the FIR filter f with random complex values, then compute the initial sensing matrix , set k = 0.
• Step B2: Compute the SVD of k g −1 k and the unitary matrix U.
• Step B3: Compute the FIR filter f k+1 that minimizes (24) by (31); under the constraint f k+1 • Step B4: Repeat the above steps until the convergence criteria is satisfied, e.g., f k − f k+1 2 < 2 , where f k+1 is the FIR filter obtained at the k + 1th iteration and 2 is a predefined threshold.
Because P LM, the SVD of the P×LM matrix g −1 in Step B2 requires a large computation amount for large values of L and M.

Joint optimization
With the above discussion, now we turn to the transmit waveform x and FIR filter f joint optimization problem. The method will be considered to combine the introduced iterative approaches for optimizing the transmit waveform and FIR filter. Considering these two variables cannot be optimized simultaneously in an iteration, we split each iteration into two parts, which optimize one variable while the other is fixed. With the transmit waveform and FIR filter optimization approaches in Sections 3.1 and 3.2, the joint optimization method can be summarized as

•
Step C1: k = 0, generate the transmit waveform x k with random complex values and constant energy, and set the FIR filter f k with random complex values; compute the corresponding sensing matrix k , the sparse representation matrix k , the equivalent dictionary D k , and the Gram matrix Step C2: Assume the deterministic sensing matrix = k , optimize the transmit waveform using Step A2 to Step A4, and obtain x.
• Step C3: With the deterministic waveform x = x, optimize the FIR filter using Step B2 to Step B4 and obtain f. • Step C4: k = k + 1, Step C5: Compute the Gram matrix Repeat the above steps until the convergence criteria is satisfied, e.g., G k − G k+1 2 < 3 , where 3 is a predefined threshold.
The proposed three algorithms in this section is stopped whenever the innovations is less than a certain value ( 1 , 2 , and 3 , respectively). The order of the magnitude of these values will be given in the simulation section. Similarly, the number of iterations will also be tested.

Simulation
In this section, we will complete computer simulations with three aspects. First, simulation examples will be given to demonstrate the effectiveness of our proposed methods for decreasing the coherence between the sparsifying representation and sensing matrices. Second, in CS theory, RIP is an important rule. Simulation results will show that our designed result can ensure this rule finely. Third, the target scene recovery experiment will be given to show the improved recovery accuracy by our methods.

Transmit waveform and sensing matrix optimization results
The CSR system transmits a waveform of length N = 19 and measures a target scene with L = 80 range and M = 1 Doppler bin. The sensing matrix compresses the received signal with the FIR filter f of length B = 40 and obtains the measured data of length P = 40. We optimize the transmit waveform x and FIR filter f of the CSR system separately and simultaneously, and results are compared for these different approaches. The parameters 1 , 2 , and 3 are set to be 10 −8 , 10 −5 , and 10 −5 , respectively. These algorithm will stop after hundreds of iterations in our simulations. The recovery algorithm used here is OMP. In [6], Herman has proved that the Alltop sequence has nearly ideal incoherence properties for the dictionary . This sequence is defined as where n ≥ 5 is a prime. Considering the important property of the Alltop sequence, we use it as a  standard sequence for comparing with our optimization result. In the process of the sensing matrix optimization, we also use this sequence as the fixed transmit waveform. Figure 3 shows the histogram of the cross-correlations between different target responses with random sensing matrix. The optimized waveform is compared with the Frank, LFM, and Alltop sequences. As a standard sequence, the maximum and average values of the histogram of the Alltop sequence obviously outperform Frank and LFM sequences. However, our proposed method can provide a sequence with better performance than the Alltop sequence. The average and maximum cross-correlations are decreased about 0.07 and 0.27, respectively.
In Figure 4a, the result for the optimized sensing matrix is given. The Alltop sequence is used as the fixed transmit waveform. As we can see, the cross-correlations of the optimized sensing matrix are obviously smaller than that of the random matrix. The average and maximum crosscorrelations are lowered about 0.17 and 0.32, respectively. The result for the optimized transmit waveform and sensing matrix are also shown in Figure 4b. Compared with that of the Alltop sequence with random sensing matrix, the average and maximum cross-correlations are decreased about 0.09 and 0.44, respectively. The detailed information of the cross-correlations is summarized in Table 1. From the above results, we may conclude that the sensing matrix optimization offers greater performance improvement than the optimized transmit waveform because it can provide more optimization variables.

RIP verification
With normalized θ T , i.e., θ T 2 2 = 1, the RIP in (2) can be simplified to We note that where λ min and λ max are the smallest and largest eigenvalues of D H T D T . The S-restricted isometry constant can be obtained by With the constraint 0 < δ S < 1 [15], the eigenvalues λ min and λ max should satisfy Thus, to verify if the equivalent dictionary D T satisfies the RIP, we should compute the eigenvalues λ min and λ max of all the possible subset matrix D T . This operation needs great amount of computation. We use 1,000 times Monte Carlo simulations instead. The smallest and 10   largest eigenvalues (λ min , λ max ) are plotted in Figure 5. As we can see, (35) holds with the sparsity K ≤ 3 for Alltop sequence and random filtering measurement matrix. However, with optimized transmit waveform and FIR filter, (35) can be satisfied until the sparsity is increased to 8. We note that this simulation does not verify the RIP but instead provides to proxy to it.

Target scene recovery
To verify the effectiveness of the transmit waveform and sensing matrix optimization in random filtering structurebased CSR system, we assume that a sparse target scene has L = 25 range and M = 25 Doppler bins. The CSR system transmits a waveform of length N = 31 and measures the received signal with a sensing matrix of P × L, where P = 25. There are two target scenes in Figure 6. The scene is recovered by the Alltop sequence and random sensing matrix or optimized waveform and sensing matrix, and this process is carried out by 1,00 times. Figure 7 shows the root mean squared (RMS) errors of the recovery results of the two target scenes versus SNR. As we can see, recovery errors with optimized waveform and sensing matrix are much smaller than the result with Alltop sequence and random sensing matrix. With increased SNR, the recovery error is obviously decreased.

Conclusions
A new notion of random filtering structure-based compressive sensing radar was proposed in this paper. To decrease the coherence between the sparse representation and sensing matrices, a computational framework for optimizing the transmit waveform and sensing matrix separately and simultaneously was introduced. We showed that optimized transmit waveform and sensing matrices lead to smaller mutual coherence between different target responses. We also use Monte Carlo simulations to verify whether our optimized results satisfy RIP. Simulation results demonstrate that our optimized results can obey RIP with much higher sparsity than nonoptimized waveform and sensing matrix. As we can see, the reconstruction accuracy was significantly improved by the optimized transmit waveform and sensing matrices for a given target scene.