Bearing estimation of low probability of intercept sources via polynomial matrices and sparse linear arrays

Recent years have seen a steady convergence of Radar and Communications band Radio Frequency (RF) transceiver systems. Not only have Communications systems colonised large swathes of previously allocated Radar bands but there is also a convergence of technologies driven by the relatively low cost of software-defined transceivers and solid-sate RF sources. Thus, where conventional radar transmissions are characterised by short narrowband pulses with high peak power, new classes of 'pulse-compression' radar are being developed to exploit this new technology. The resulting Low Probability of Intercept waveforms are designed to spread energy in both time and frequency, yielding a very low instantaneous power spectral density. Methods to detect, analyse and distinguish such sources require longer acquisition periods to collect more energy from the sources. Here, a novel solution is provided for detection and separation based on direction-finding utilising polynomial matrix methods in conjunction with sparse array geometries. This approach provides enhanced detection, separation and direction finding while using relatively few antenna elements.


| INTRODUCTION
The advent of spread-spectrum, low probability of intercept (LPI) radar systems and the increasing overlap of Radar and Communication system bands is providing increasing challenges for passive receiver systems. Such radio frequency (RF) systems exhibit much lower peak power, making even detection more difficult. Moreover, even when detection is possible, the convergence of the underlying technologies around solidstate sources and software-defined signal structure makes distinguishing Radar from Communications signals much more difficult [1,2]. This technology convergence is likely to be an increasing issue that is compounded by the ever denser occupation of RF bands leading to cognitive systems concepts characterised by non-stationary signal characteristics [3].
Spread-spectrum RF waveforms provide an efficient means of exploiting the available bandwidth, with pulsed waveforms being replaced by phase or frequency coding superimposed on Continuous Wave (CW) waveforms. The key characteristic of such waveforms is the spreading of energy in both time and frequency yielding a low power spectral density (PSD) [4]. This results in the need both for much longer signal integration times and new processing methods that no longer rely on fast Fourier transform (FFT) processing.
There are many different schemes for producing spreadspectrum waveforms. Perhaps the oldest and simplest are frequency-modulated continuous waveforms (FMCWs) [5], which are still used in many Radar systems. Here, a very broadbandwidth signal can be generated with low hardware and software costs since the receivers need only simple homodyne approaches whereby the received signal is mixed with the transmitted signal, and thus only beat frequencies need to be digitised [6]. In more modern systems modulating waveforms are synthesised purely digitally, promoting flexibility. Thus in direct sequence spread spectrum radar, the approach mimics that of digital communication; using a method originally devised to improve interference and multipath robustness [7]. This technique spreads the signal energy across a wide range of frequencies via some modulation with pseudo-noise codes.
To detect and locate such sources, broadband antennas and subsequent analogue and digital signal processing circuitry are required. However, these are expensive in terms of both monetary costs, and power consumption, the latter leading to further issues such as heat dissipation. Although the obvious solution is to use only a small number of antenna elements, this compromises performance. With array-element spacing small enough to provide unambiguous angle resolution, the array aperture is small and the angular resolution is low. Moreover, close element spacing will also tend to exacerbate mutual coupling between adjacent antennas. While these issues might be addressable for narrow band sensors, they are made much worse on broadband systems where avoiding ambiguity at the highest frequencies will seriously compromise angular resolution for lower frequency signals.
In this paper, the aperture problem is addressed by using sparse linear array geometries. While such sparse array designs such as co-prime [8], nested [9], super-nested [10] and minimum redundancy arrays (MRAs) [11] are much exploited in the context of narrowband signals, they have not previously received much attention in broadband systems. Traditionally, for an array of N elements, subspace-based methods such as MU-SIC can resolve a maximum of N − 1 sources. The virtual array methods presented in this paper are not bound by this limitation and have the ability to resolve more sources than sensors, due to increased degrees of freedom gain from this processing [12].
To deal with the broadband signals considered here, this paper exploits relatively recently developed polynomial matrix methods that have demonstrated promising solutions for a range of broadband beamforming and direction of arrival (DoA) estimation problems [13][14][15][16]. Such methods have demonstrated improved performance over independent frequency bin methods [14], and the spatio-spectrum generated provides a joint frequency-angle estimate, rather than simply an angle-only estimate as is the case in other well-known broadband DoA methods such as coherent signal subspace processing [17,18]. In applying these methods we extend the usual analysis, based on directly sampled signals, to allow application to downconverted signals that are characteristic of the RF world.
Presenting a novel approach to spread spectrum radar direction-finding by exploiting polynomial matrix methods, in conjunction with sparse arrays, this paper specifically addresses the following: � Application to realistic scenarios where it is necessary to include a frequency down-conversion stage in the signal model. � Practical analysis of virtual broadband uniform linear sensor arrays using polynomial matrix methods. � Clarification of the spatial smoothing stage that is necessary for constructing non-singular covariances for sparse array analysis. � Redundancy averaging in the covariance matrix construction to improve spatio-spectral analysis accuracy. � Polynomial MUSIC implementation without the need to use any fractional delay filter approximation. � Statistical analysis of the performance in terms of accuracy, and probability of separating close sources.
The remainder of the paper is organised as follows. Section 2 introduces the history and theory behind select sparse arrays, as well as the reasons for improved performance. Section 3 discusses the model for broadband signals illuminating an arbitrary array, with complex down-conversion. The motivation and theory of polynomial space-time covariance matrices are discussed in Section 4, as well as the polynomial eigenvalue decomposition (PEVD). The spatio-spectral polynomial-MU-SIC (PMUSIC), and Root-MUSIC algorithms are discussed in Section 5. The performance of the proposed methods is analysed in Section 6, and finally a summary of conclusions and future works in this area is discussed in Section 7.
To conform with standard notation, vectors and matrices are denoted via lower-and uppercase bold letter, for example a and A. The convolution operation is denoted by the ⊗ symbol, and the Kronecker product is denoted via the ⊙ symbol. The complex conjugate is denoted by a superscript * , for example a*, the Hermitian transpose of a vector or matrix is represented with a a H ; the para-Hermitian transpose of a polynomial vector or matrix is represented as a P (z), and is defined as a P (z) = a H (z −1 ).

| SPARSE LINEAR ARRAYS
The linear array has played a crucial role in many areas of sensing, including radar [5] and electronic warfare [19]. The uniform linear array (ULA) is perhaps the simplest linear array. However, since accuracy and resolution is dependent on the electrical length of array, using a ULA to cover a wide range of frequencies makes little sense since many antenna elements are required to have an array with sufficient electrical length for adequate accuracy and resolution for lower frequency sources, while still maintaining ambiguity-free localisation of higher frequency sources. Sparse arrays such as co-prime, nested, super nested and MRAs can be implemented to cover a wide aperture, and thus sufficiently large electrical length for all frequencies of interest, with a relativity low number of physical antennas. Figure 1 shows the physical array geometries for an eight element uniform linear, nested, super-nested, and restricted MRA geometries.
Consider a set of N sensors placed on a uniform grid of spacing d = λ min /2, where λ min is the wavelength of the highest frequency source anticipatedensuring unambiguous bearing estimation (i.e. 25 mm for an X-band radar). The sensor x n is physically located at nd from the first sensor in the array. The difference set of the array is defined as The significance of this difference set is that it defines the spatial lag for which second order statistics can be estimated under wide sense stationary conditions, which arises naturally when estimating a spatial covariance matrix: where n is the discrete time index, and τ is the discrete lag parameter.
The weight function w(u) of a sparse array is defined as the number of times that the difference u occurs within its difference set. It is to be noted that this weight function is not to be confused with the weight vector commonly used in beamforming algorithms. An example shown in Figure 2(a) demonstrates the weight function for an eightelement ULA.
If the same difference, u, occurs more than once, that is w(u) > 1 then it said to be redundant. In Section 4.1, it is shown that the redundancies can be used to enhance the estimate of the virtual space time covariance matrix.

| Minimum redundancy array
A class of linear arrays was shown in [11] which achieves maximum spatial resolution for a given number of elements by minimising the redundancies in the difference set. Whilst the original scope of [11] is for radio astronomy, the aim of detecting and locating faint distant sources is very similar to an electronic surveillance problem.
It is worth noting that the case of the zero redundancy array is the goal of a MRA. However, such an array does not exist for more than four elements [20]. Calculating the sensor locations of a MRA is not a trivial task as closed form expressions for estimating sensor locations do not exist, however, look up tables and examples are provided in [11] and some methods can be found in [21]. In this paper, we only consider the case of a restricted array, in which the difference set fills the holes in the entire array aperture. An example of an eight element MRA is provided in Figure 1 and its weight function in Figure 2(b). There is also no closed form expression for the maximum number of sources such an array can resolve with subspace methods. However, in the example provided, it is shown that the 8-element restricted MRA can resolve a maximum of 23 sources since a 24 element virtual ULA is generated in later processing stages.

| Nested, and super nested array
In comparison with the MRA, it is trivial to calculate the sensor positions for a nested array. It is the union of two ULAs; a Nyquist spaced (d ) array with N 1 elements, and a sparse ULA of N 2 sensors spaced at (N 1 + 1)d, over an aperture of N 2 (N 1 + 1)d, yielding a sparse array with N 1 + N 2 sensors overall, with the ability to resolve N 2 (N 1 + 1) − 1 sources. An 8-element (N 1 = N 2 = 4) Nested array is shown in Figure 1, and its weight function in Figure 2(c). Unlike a co-prime array [12], this array yields a contiguous difference set, akin to the restricted class of MRAs [22]. While such an F I G U R E 1 Physical Array Geometries for the uniform linear-, nested-, super nested-and minimum redundancy arrays array is easy to design, in general, its weight function contains more redundancies than an MRA, and another major drawback to the conventional nested array is the issue of mutual coupling between sensors of the dense Nyquist portion of the array.
Recently, the super-nested array was proposed in [10,23] which aims to reduce the effect of mutual coupling in the array by redistributing the dense portion of the array across the entire aperture. The goal of the second-order supernested array is to minimise the pairs spaced at Nyquist, that is the w(1) = 1, which will reduce the effect of mutual coupling substantially, and a third-order super-nested array also minimises the pairs spaced at 2 Nyquist. Figure 1 shows an 8-element second-order super-nested array. Its weight function in Figure 2(d) demonstrates that while there is the same redundancy as the conventional nested array, the number of pairs spaced at Nyquist (u = ±1) is reduced from 4 to 2.

| DATA MODEL
As a concrete example of a typical broadband source, we consider the case of signals with an FMCW waveform. We also imagine that this signal has been down-converted to some lower frequency to allow unambiguous digitisation with only modest receiver bandwidth. Nonetheless, for the wide frequency modulation of LPI radar, the narrow-band array approximation is no longer valid. Thus, it is necessary to consider the signals at the separate array elements in terms of time delay and not simply phase shift.

| FMCW radar source model
A broadband frequency-modulated continuous wave radar signal is a common example of a LPI signal as it results in a low instantaneous PSD. To the radar, the waveform is known and integration times can be chosen so that a high enough total energy can be obtained to allow detection. Such waveforms are indeed used in many radar applications [5]. The source model for a triangular linear frequency modulated signal is where α is the signal amplitude, f c is the carrier frequency, b is the signal bandwidth and τ c is the chirp duration.

| Antenna signal model
Consider such a source, s(t), illuminating a linear array. Assuming the source is in the far field, and the array elements are identical and isotropic, then the signal at the antenna array is where τ m is the inter-element time delay from the reference antenna, (position zero) and is τ m = S m d sin θ/c, where S m denotes the integer valued position from the reference, d is the Nyquist rate spacing (λ min /2), λ min being the smallest wavelength anticipated, θ the DoA, c being the propagation speed and ν is the additive noise vector. In many areas of array signal processing, the signal is assumed to be narrowband, that is the time period of the complex envelope is considerably greater than the time it takes the wavefront to traverse the entire array aperture. Owing to this narrowband approximation, it can be modelled simply as a time shifted sinusoid, that is a phase shift. The signal seen at the antenna can thus be modelled as an instantaneous mixture of a steering vector at the signal's centre frequency, and the narrowband source signal, s α (t) Assuming perfectly calibrated hardware, the baseband signal sampled at t = nT s can be modelled as b where ω d is the frequency of the down-converting oscillator and ν(n) is the receiver noise. Now, consider a broadband source, s β (t). Clearly, this narrowband approximation is invalid and the time delays in (4) cannot be modelled as a simple phase shift. It can still be factorised in terms of a steering vector multiplied by the source signal, provided this multiplication is interpreted as a convolution. Thus, in terms of the continuous time Dirac delta function δ(t), and the following convolutive model, and thus the model of the digitised down-converted signal b where b s β ðnÞ is the down-converted source signal, e δðn − e τÞ represents an ideal fractional delay FIR filter, and e τ is the normalised delay in samples. In reality, such a filter would be infinite in length and thus non-causal and non-realisable. There are multiple methods to approximating such a filter, including simple windowed sinc functions [24], modified Farrow structures [25] and filter bank methods [26].

| POLYNOMIAL SPACE-TIME COVARIANCE MATRICES
In the case of a narrowband scenario in (6), owing to the instantaneous mixture model, only instantaneous temporal correlations are of interest, thus the common step in many narrowband array signal processing problems is to form a data covariance matrix at the lag of k = 0.
However, in the case of broadband signals illuminating the array, it is now necessary to consider a range of temporal correlations due to the convolutive mixture model in (9). Thus, the polynomial space-time covariance is defined as where W should be sufficiently large that R xx (±W) is approximately a zero matrix. Assuming a good estimate of the polynomial covariance matrix, this may also be expressed as [15]: where A(z) = [a 1 (z), a 2 (z), …, a L (z)] is the steering matrix, R ss (z) is z-transform of the source covariance matrix R ss (k) = E [s(n)s H (n − k)], and σ 2 ν is the noise variance. Without loss of generality, we make the simplifying assumption that all sources illuminating the array are uncorrelated, thus the source covariance matrix will be diagonal, where 〈R ss ðzÞ〉 ii ¼ σ 2 i ðzÞ is the z-transform of the autocorrelation function of the i th source. This assumption is not strictly necessary for the main results of this paper, but it is useful in providing a more easily understood description of the sparse array analysis. In the case of highly correlated or even coherent sources, there will be significant energy in the off-diagonal elements of this matrix. Thus, the assumption that the source covariance matrix is diagonal is invalid in the case of highly correlated sources. We can rewrite (12) as This model is applicable, whatever the sensor configuration, but to exploit the features of particular sensor configurations, we now need to consider the structure of the steering vectors in more detail. Now, if we consider a linear sensor array; then, from (9) we note that in the time domain, the rows of the steering vectors are convolutional products of the primitive (Nyquist) ideal, unitary, fractional delay filter, characterised by a delay τ l , say, for source l. Thus, in the zrepresentation: � The rows of the steering vectors are just powers of the primitive (Nyquist) steering element. that is a nl (z), the n th row of the l th source steering vector is: � The para-Hermitian conjugate of the steering vector is just the inverse of the steering vector: Of course for the sparse arrays described in Section 2, the set of M integers, {m n : m n ∈{1: M}} is sparse but still provides a complete set of integer sensor spatial displacements in the set where M >> N, this is exploited below in Section 4.1 where a virtual covariance matrix is constructed based on this set M of sensor pairs. COVENTRY ET AL. -5

| Formation of virtual covariance matrices
The definition of the difference set was given in (1). The elements of the polynomial covariance matrix contain all spatial differences from this set, thus the component 〈R xx ðzÞ〉 ij represents the spatial auto-(for i = j) and cross-(for i ≠ j) correlations for spatial difference i − j. By vectorising this matrix, a virtual array is generated with positions contained in the difference set.

γðzÞ ¼ vecðR xx ðzÞÞ ð17Þ
Using the simplification from (13), this can be rewritten as Here, e 1 is the vectorised identity matrix that is a vector of zeros with ones at N + 1 intervals, corresponding to the diagonal, auto-correlation components, of R xx (z). Clearly the first term can be rewritten in matrix form to give something analogous to the original sensor model of (9): where σ 2 (z) is a vector of source autocorrelation functions, and represents the steering vector for the much larger virtual array of sensor data pairs, covering all the sensor displacement in the difference set M. Clearly, this representation contains redundancies as spelt out in the weight functions of the difference set, as introduced in Section 2. In particular, for all the sparse arrays of Section 2 there is an N-fold redundancy of the zero displacement auto-correlation terms. We can truncate the arrays in (19) to eliminate the redundancies. If we assume each element is corrupted by some zero mean independent and identically distributed noise process, then rather than simply removing the redundant elements they can be averaged to enhance the estimate of the virtual ULA and in the process reduce the N 2 size of γ(z) to the (2M − 1) element ULA: Here, e e is a vector of zeros, except a one at the array position zero. If we re-order the rows of B(z) so that the differences are in order of increasing sensor displacement then, using the known steering vector structure summarised in (14) and (15), we are able to obtain the virtual array steering vector in the rationalised form characterised by the Vandermonde structure: With this condensation and reordering of (19), A 1 (z) now has precisely the structure of a conventional ULA steering vector and the model expressed in (21) is almost identical in structure to that for the ULA summarised in (9).
With that in mind, the obvious next step is to construct the quadratic form: Not surprisingly, this has just the form of the covariance (12) and just as in the case of (12) the rank of γ 1 ðzÞγ P 1 ðzÞ will be equal to the rank of the source matrix σ 2 (z)σ 2 P (z). However, we can immediately see that: 1. Since the source matrix is simply the outer product of two vectors, the source matrix is seriously rank deficient; being rank 1, irrespective of the number of sources. 2. Moreover, the noise term is also of rank 1, since e e and e e P are also just column and row vectors, respectively, each being zero except for a single 1 at the zero displacement element.
Thus, despite the familiar form, (17) is not a representation that enables the usual forms of quadratic signal analysis. Nonetheless, the situation can be recovered as described below.

| Rank Recovery by spatial smoothing
Spatial smoothing [27] is a well-known method of restoring the rank of this source covariance matrix; an approach that was recently extended to the case of broadband sources [28]. In this approach, applicable to ULAs, a new array is constructed by taking a shortened sub-array and forming the average of this array with translations of it. In general, averaging over just L such translations is sufficient to increase the rank to L, allowing L sources to be resolved.
In [8] it was demonstrated that by forming M overlapping subarrays of M elements could restore the rank of (23). It was suggested this step was required to restore the rank of the source covariance matrix. While this is true, the spatial smoothing theory stated above requires only L subarrays to restore the rank of this matrix, rather than M. In our findings, it was noticed that M subarrays were required to achieve a valid space time covariance matrix. Consequently, not only does (23) exhibit a singular source matrix but the noise is also singular. Since we assume L < M, to remove both singularities it is necessary to form an M-dimensional sub-array obtained as the average of its M possible translations.
Starting with the shortened array defined by the M-dimensional steering vector, A 0 (z), given as the lowest M rows of A 1 (z) in (22), the M translations of this ULA are: where Q(z) is the diagonal matrix, and e e i is a vector of zeros apart from a single 1 at array index i.
We are now in a position to form the covariance of this spatially smoothed, M-dimensional, virtual ULA: Inserting (20), we obtain the full structure of this virtual array covariance: The additional phase terms in the vector outer-product terms ensure incoherent summation that restores the source covariance matrix rank, as it can be proved, just as in [28]. It is also easy to see that the translational shift, in the summation, steps the single 1 along the noise matrix diagonal; restoring it also to full rank. Thus I in (27) is just the M-dimensional identity matrix.

| POLYNOMIAL MUSIC
The MUSIC algorithm [29] has been a powerful tool for super resolution DoA estimation of narrowband sources for decades, and has recently been extended to broadband scenarios [14] due to the advent of the PEVD, [30][31][32].
Consider the decomposition of the space-time covariance matrix from (27) is spectrally majorised such that λ 1 (e jω ) > λ 2 (e jω ) > … > λ M (e jω ), where where K is the polynomial length and k represents the polynomial index. Thus, the eigenvalue PSD can be estimated through the use of a Fourier transform. From (27), it is clear that the number of significant eigenvalues, and thus the dimensions of the signal subspace will be the rank of e R ss , which will be equal to the number of sources illuminating the array. This can thus be separated into signal, and noise subspaces It is clear that the eigenvectors related to the signal subspace, U s , contain information of the steering matrix, A 0 . However, the columns of U s will be orthonormal, whereas the columns of the true steering matrix may not be. The columns of the true steering matrix, however, will be orthogonal to the noise subspace.

| Spatio-spectrum polynomial MUSIC
Consider the polynomial where a(θ, z) is the test polynomial steering vector for DoA θ. This test steering vector is a vector of ideal fractional delay FIR filters, as described in (9). The ideal filters are, however, non-causal, and thus non-realisable. Realisable COVENTRY ET AL.
approximations do exist, such as the methods used in [26]. However, since the goal is to estimate a discretely sampled spatio-spectrum, the above polynomial can be evaluated in discrete points for z = e jω , Γðθ; e jω Þ ¼ a H ðθ; e jω ÞU ν ðe jω ÞU H ν ðe jω Þaðθ; e jω Þ ð34Þ and now the test steering vector becomes the same as the narrowband case across a range of frequencies.
aðθ; e jω Þ ¼ This test steering vector is then scanned across the noise subspace in discrete spatial (θ) and spectral (e jω ) steps, and the spatio-spectrum is calculated as and the peaks denote the spatio-spectral location of the sources.

| Root spatio-spectrum polynomial MUSIC
The scanning portion of thePMUSIC algorithm is computationally expensive. It was shown in [28] that it was possible to extend the computationally cheaper Root-MUSIC algorithm to broadband scenarios by reformulating the problem. Similar to the above example, the goal is to estimate the polynomial steering vector, a(z) to solve Using the identities from (14) and (15), this can be rewritten as Exploiting the Vandermonde structure of the steering vector for the ULA, we can express Γ(z) as a Laurent polynomial of a 1 (z) with 2M − 1 coefficients, where the coefficients are simply the sum of the diagonal elements of C(z).
where b i (z) is the sum of C(z) down its i th diagonal, with 0 th being the main diagonal. Note that since Γ(a 1 [z]) is a polynomial of a 1 (z), and a 1 (z) itself is a polynomial of z. As such, conventional polynomial root-finding methods cannot directly be applied to estimate the filter solutions for a 1 (z) and the resulting directions of arrival. Instead, similar to (34), the polynomials of z can be evaluated for z = e jω Γða 1 ðe jω ÞÞ ¼ Thus, conventional polynomial root-finding methods can be used at each discrete frequency point of interest. Since this is in discrete frequencies, the steering vector will be just as in (35). Thus the directions of arrival for the L sources will be given by the L roots lying closest to the unit circle. From the roots, the actual directions of arrival is calculated as where q l (e jω ) is the root relating to the l th source.

| PERFORMANCE ANALYSIS
In this section, the performance of the above methods is demonstrated, and analysed. In all these simulation scenarios, we consider the case of uncorrelated sources and a perfectly calibrated array, each with eight elements. The PEVD algorithm used is the SMD method [31]. Section 6.1 demonstrates the powerful ability of estimating a spatio-spectrum of several spectrally overlapping, spatially close broadband sources using a uniform linear-, super-nested, and MRAs. Section 6.2 analyses the accuracy of the Root spatio-spectrum estimator when used in conjunction with the array geometries discussed and demonstrates the benefit of redundancy averaging. Section 6.3 shows the empirical resolution of the three array geometries.

| Detection & Spatio-Spectrum Estimation
Six temporally overlapping FMCW waveforms of varying carrier frequencies and bandwidths, linearly spaced between −35°a nd 35°are present, with an acquisition period of N s = 20,000 samples. Table 1 shows the ground truth for all simulated sources, and Table 2 shows the receiver properties. Since the polynomial eigenvalues are representative of a PSD, this can be used to determine the dimensions of the signal and noise subspaces. Figure 3 shows the eigenvalue PSDs for this simulation, note only the first 8-eigenvalues are displayed for the MRA and super-nested array. From these figures, it is easy to see that there are 6 significant eigenvalues, and thus 6 uncorrelated spectrally overlapping sources present. Separating the noise subspace, the spatio-spectrum polynomial music algorithm from Section 5 is used to estimate the spatio spectrum. Figure 4 shows the estimated spatio spectrum via the uniform linear, super nested, and MRAs, and a normalised 1D projection (spatial) of these can be seen in Figure 4(d). It is easy to see that all array geometries correctly estimate both the frequency, bandwidth and DoA of all sources present. However, results are clearer in the case of the sparse arrays due to the increased resolution offered, given by the much larger aperture. While these examples give a vague idea on performance, the subsequent simulations analyse the statistical performance in terms of both estimator variance, and resolution.

| Accuracy
In this simulation, we compare the DoA estimation variance across a range of signal to noise ratios (SNRs). A Monte-Carlo simulation of I = 500 runs for each SNR, with a fixed bandwidth and randomised DoAs for each run. The Root Polynomial MUSIC algorithm from Section 5.2 is used for improved computation times; in [28], it was noted that the statistical performance was largely the same as the heuristic searching method in Section 6.1. Since this is an unbiased estimator, the estimator variance, σ 2 E , is simply calculated as the mean squared error, that is where e θ i is the estimated DoA, and θ is the ground truth. Figure 5 shows the simulated performance of root polynomial MUSIC algorithms for the three array geometries. Results demonstrate that for a fixed number of antenna elements, the sparse array geometries have an improved performance over a ULA. It is to be expected that the minimum redundancy array should perform better than the super-nested array due to the greater degrees of freedom offered, and this is confirmed in Figure 5(a). The effect of redundancy averaging on the variance can be seen in Figure 5(b). This averaging has a clear improvement in estimator accuracy, and intuitively, this has a greater effect on the super-nested array as this array geometry has more redundancies in its weight function.

| Probability of resolution
Since the methods proposed in Section 6.1 are of superresolution, the resolution is not bound by the Rayleigh criterion. In this section, a Monte-Carlo simulation of 100 runs is performed to analyse the empirical probability of resolving two sources spaced at Δθ = θ 1 − θ 2 by using the spatial-only function of the Polynomial MUSIC algorithm [33].
where N ωL and N ωH are the indexes of the frequency bins containing sufficient energy. The following inequality is used to determine whether a pair of sources are resolved At the point P mu ð θ 1 þθ 2 2 Þ ¼ 0, the peaks merge into one, and the DoAs will not be resolved. The probability of a single run is calculated as Resolution analysis was performed with two FMCW sources centred at 10 GHz, with 2 GHz bandwidth. Two scenarios are simulated, one where the received signals have equal power, and another where there is a 10 dB difference. The received SNR was set to 15 dB and the redundancy averaging method was used for the formation of the space-time covariance matrices. Figure 6  Signal to noise ratio (SNR) 10 dB COVENTRY ET AL. -9 shows the simulated probability of resolution over a range of angular separations. It is easy to see that the sparse methods perform considerably better in both scenarios. In the case of equal power sources, both the minimum redundancy, and supernested array demonstrate sub 1°resolution, as opposed to the 2°f rom the ULA. In the case of 10 dB power difference between the sources, the resolution improvement using sparse arrays is further improved over the ULA. More interestingly, there is now a greater difference between the minimum redundancy and super-nested arrays, demonstrating how the increased degrees of freedom translate to a greater angular resolving power. While all arrays have equal number of elements, it is expected that the F I G U R E 3 Frequency domain representation of the polynomial eigenvalues in the case of the three array geometries F I G U R E 4 Spatio-spectrum estimation of six broadband frequency-modulated continuous waveform sources 10sparse methods improve upon performance not only due to the increased aperture, but also due to increased degrees of freedom.
Since the MRA has greater freedoms than the super-nested array, this is also expected to perform better.

| Remarks on computational complexity
The results in this section have demonstrated the clear performance advantages of having increased degrees of freedom from sparse arrays and virtual array generation, even though the number of array elements remain constant. The computational cost of the multiple shift maximum element sequential matrix diagonalisation polynomial eigenvalue decomposition algorithm scales as O(M 3 W), where M is number of sensors (or number of virtual sensors in the case of the sparse arrays) and W is the polynomial length determined by (11) [34]. For the analysis in this paper, the polynomial matrix length was constant for all array geometries, and M is 8, 20 and 24 for the ULA, super-nested and MRAs. Thus, this improvement in performance comes with an increase in computational cost compared with a ULA with the same number of elements.

| CONCLUSION
In this paper, we have presented a novel solution to the direction-finding aspect of locating FMCW LPI radar sources F I G U R E 5 Root PMUSIC estimator variance for the three array geometries in the case of the three array geometries F I G U R E 6 Probability of resolution for the three array geometries COVENTRY ET AL. that provide impressive accuracy and resolution, even at low SNRs. Results demonstrate that sparse arrays used in conjunction with polynomial matrix methods provide an elegant and promising solution. Results also demonstrate that the older method of generating a restricted sparse array, the MRA, generally out performs the more modern super-nested array due to the larger electrical length of the array and thus greater degrees of freedom generated. The effect of redundancy averaging further improves accuracy of the estimators by enhancing the estimate of the polynomial space-time covariance matrices. Future works involve expanding upon this framework to include detection capability and analysis when used in conjunction with other LPI radar detection techniques such as time frequency analysis for a more complete electronic support system.