DIRECTION OF ARRIVAL ESTIMATION

Abstract

Iterative methods for direction of arrival estimation of a signal at a receiver with a plurality of spatially separated sensor elements are described. A quantized estimate of the angle of arrival is obtained from a compressive sensing solution of a set of equations. The estimate is refined in a subsequent iteration by a computed error based a quantized estimate of the direction of arrival in relation to quantization points offset from the quantization points for the first quantized estimate of the angle of arrival. The iterations converge on an estimated direction of arrival.

Claims

1. A method for estimating a direction of arrival of a signal, comprising: receiving, at a processor, data indicative of a set of complex envelope voltage outputs from an antenna element array comprising a plurality of spatially separated antenna elements; applying, by the processor, a computational technique for solving an under-determined set of linear equations to the outputs and a plurality of grids of candidate directions of arrival to generate a plurality of solutions, the plurality of solutions indicating a first interim direction of arrival estimate in relation to a first grid of the plurality of grids of candidate directions of arrival and a second interim direction of arrival estimate in relation to a second grid of the plurality of grids of candidate directions of arrival, wherein the first grid and the second grid include a number of grid points that is higher than the than outputs and are rotated relative to each other; generating, by the processor, a first angular discriminant based on the first interim direction of arrival estimate and the second interim direction of arrival estimate; and generating, by the processor, a direction of arrival estimation between the first interim direction of arrival estimate and the second interim direction of arrival estimate based on the first angular discriminant.

2. The method of claim 1, wherein the plurality of grids of candidate directions of arrival includes a third grid, rotated relative to the first grid and the second grid and wherein plurality of solutions indicate a third interim direction of arrival estimate, and wherein the method includes, prior to generating the first angular discriminant: generating a second angular discriminant based on the third interim direction of arrival estimate and the first interim direction of arrival estimate; and determining that the second angular discriminant does not satisfy a threshold condition, wherein generating the first angular discriminant is responsive to the determination that the second angular discriminant does not satisfy the threshold condition.

3. The method of claim 2, wherein the first grid and the second grid are rotated relative to each other by an amount determined, by the processor, based on the second angular discriminant.

4. The method of claim 1, wherein the plurality of grids of candidate directions of arrival comprise at least three grids of candidate directions of arrival and wherein the method includes iteratively applying, by the processor, the computational technique to the outputs and pairs of candidate directions of arrival to generate solutions and generating an angular discriminant for each pair of candidate directions of arrival until the generated angular discriminant satisfies a threshold condition, wherein the first angular discriminant satisfies the threshold condition.

5. The method of claim 1, wherein the angular discriminant is based on relative magnitudes of the solutions indicating the first interim direction of arrival estimate and the second interim direction of arrival estimate.

6. The method of claim 1, further comprising producing a null in a received antenna pattern, wherein the null is located in the antenna pattern based on the generated direction of arrival estimation.

7. The method of claim 1, further comprising determining, based on the generated direction of arrival estimation, that a transmission associated with the set of complex envelope voltage outputs is from an authorized user.

8. The method of claim 1, further comprising determining, based on the generated direction of arrival estimation, that a transmission associated with the set of complex envelope voltage outputs is from an unauthorized user.

9. The method of claim 1, wherein the plurality of spatially separated antenna elements includes a plurality of antenna elements spatially separated across two orthogonal dimensions and wherein the generated direction of arrival estimation is a first dimensional estimation, in relation to one of the two orthogonal dimensions.

10. The method of claim 9, further including repeating the method for the other of the two orthogonal dimensions to generate a second dimensional estimation.

11. The method of claim 10, further comprising generating a three-dimensional direction of arrival estimation based on the first dimensional estimation and the second dimensional estimation.

12. The method of claim 11, wherein the first dimensional estimation is an estimation of azimuth and wherein the second dimensional estimation is an estimation of elevation.

13. A method for estimating a direction of arrival of a signal comprising: receiving, at a processor, input signals representative of detection of one or more signals received at a respective plurality of spatially separated sensor elements; in an initial determination and in at least a first iteration determining, by a processor, based on phase information in the input signals and a known geometry of the spatially separated sensor elements, at least one sparse solution indicating one or more estimated directions of arrival amongst a set of candidate directions of arrival, wherein for each iteration the set of candidate directions of arrival are rotated, and wherein each iteration is performed based on preceding sparse solutions to cause the iterations to display convergence in the sparse solutions; and determining by the processor, based on preceding sparse solutions, that a convergence threshold for the sparse solutions has been satisfied and in response generating data indicating a direction of arrival estimation for the signal based on the immediately preceding sparse solutions.

14. The method of claim 13, further comprising determining that the convergence threshold for the sparse solutions has been satisfied based on relative magnitudes of the preceding sparse solutions.

15. The method of claim 13, wherein the spatially separated sensor elements comprise an array of antenna elements of an antenna and the method further comprises producing a null in a receive antenna pattern of the antenna, wherein the null is located in the antenna pattern based on the generated direction of arrival estimation.

16. The method of claim 13, further comprising determining, based on the generated direction of arrival estimation, that a transmission associated with the set of complex envelope voltage outputs is either from an authorized user or from an unauthorized user.

17. The method of claim 13, wherein the plurality of spatially separated sensor elements is a first set of sensor elements and the direction of arrival estimation is in relation to azimuth or elevation and comprising repeating the method for a second set of spatially separated sensor elements to generate another direction of arrival estimation in relation to elevation or azimuth respectively.

18. An iterative method for direction of arrival estimation of a signal at a receiver with a plurality of spatially separated sensor elements, the method comprising: generating a first quantized estimate of the angle of arrival from a compressive sensing solution of a set of equations relating sensor output signals from the sensor elements to direction of arrival; and refining the first quantized estimate in at least one subsequent iteration by a computed error based a quantized estimate of the direction of arrival in relation to quantization points offset from the quantization points for the first quantized estimate of the angle of arrival, wherein the offset is selected to cause the iterations to converge on an estimated direction of arrival.

19. The method of claim 18, wherein the plurality of spatially separated sensor elements comprise an array of antenna elements of an antenna and the method further comprises producing a null in a receive antenna pattern of the antenna, wherein the null is located in the antenna pattern based on the generated estimated direction of arrival.

20. The method of claim 18, further comprising determining, based on the generated estimated direction of arrival, that a transmission associated with the set of complex envelope voltage outputs is either from an authorized user or from an unauthorized user.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0097] FIG. 1 illustrates a block diagram of an embodiment of an iterative process for DOA estimation.

[0098] FIG. 2 illustrates a block diagram of another embodiment of an iterative process for DOA estimation.

[0099] FIG. 3 illustrates a contour plot of mutual coherence for UCA as an antenna geometry.

[0100] FIG. 4 illustrates a line plot illustrating the impact of the radius of an UCA on the mutual coherence and condition number of the measurement matrix Φ.

[0101] FIG. 5 illustrates a line plot of the impact of external noise in DOA estimations.

[0102] FIG. 6 illustrates a line plot of the convergence of the algorithm through a first DOA estimate and through two iterations. The Cramer-Rao Lower Bound (CRB) is also plotted.

[0103] FIG. 7 shows a flow diagram of a process to perform DOA estimation.

[0104] FIG. 8 shows a flow diagram of another process to perform DOA estimation.

[0105] FIG. 9 shows a block diagram of a radio receiver for DOA estimation.

[0106] FIG. 10 shows an MSE performance comparison of an embodiment of the disclosed method using UCA as an antenna geometry against other DOA techniques, Root-MUSIC and Beam-Forming, for antenna geometry of a uniform circular array. The inter-element spacing is

[00009]$\frac{\lambda }{2}.$

[0107] FIG. 11 shows an MSE performance comparison of an embodiment of the disclosed method using ULA as an antenna geometry against other DOA techniques, ESPRIT, Root-MUSIC and Beam-Forming, for antenna geometry of a uniform linear array. The inter-element spacing is

[00010]$\frac{\lambda }{2}.$

DETAILED DESCRIPTION

[0108] A method for direction of arrival (DOA) estimation involves the determination of the DOA of a signal, for example a radio signal, from a measured characteristic of the signal, for example from signal complex voltages at the outputs of antenna elements configured in an array. In other examples the signal is a sound signal detected by an array of sound detectors, for example an array of microphones. In other examples the signal is a phase coherent light signal detected by an array of photodetectors.

[0109] The array may be 2-dimensional or 3-dimensional. The range of potential angles of arrival is quantized into a grid, in which the number of grid points in the quantization is greater than the number of antenna elements, whereby a sparse recovery problem is created. The grid may have points that are uniformly spaced or non-uniformly spaced. For clarity of illustration the description herein is given with reference to a uniform grid, which in at least some embodiments provides an advantage of requiring reduced computational resources. The grid may extend around all directions or occupy only a subset of directions.

[0110] In some embodiments the array of antenna elements is a circular array, which in one embodiment is a uniform circular array (UCA). In other embodiments, the array is a linear array, for example a uniform linear array. The disclosed methods may be used with other forms of uniform and non-uniform arrays. The described technique is applicable to any antenna array geometry.

[0111] In some embodiments the measurements are transformed using a transformation function that increases sparsity. For example, in one embodiment the measurements are transformed by a decorrelating transform, which may be an orthogonal transform. In other embodiments the transformation step or process is omitted.

[0112] The signals are processed using a computational technique for solving an under-determined set of linear equations. In some embodiments compressive sensing (CS) is used to determine a DOA estimate. Embodiments of compressive sensing are described, for example, in United States Patent publication number 2006/0029279 A1 (Donoho), which is hereby incorporated herein by reference. In some embodiments CS algorithms utilise basis pursuit or another greedy algorithm such as CoSaMP. CS, in particular CoSaMP is the basis of the examples provided herein, but as noted above other techniques for solving an under-determined set of linear equations may be used. The solution indicates the grid point with the greatest signal magnitude and this grid point corresponds to at least an initial estimate of the angle of arrival.

[0113] In some embodiments CS is applied with respect to antenna array output and both a first grid and a second grid, different from the first grid. In one implementation the second grid is or can be viewed as an electronic rotation of the first grid, with grid points rotationally offset (frame of reference being potential DOAs) from the grid points in the first grid. CS for the second grid may use the same measurements, e.g., the same set of complex envelope voltage outputs, as were used for the first grid.

[0114] In some embodiments CS is performed incorporating a fixed or a variable phase shift from a previous CS determination, for example as described with reference to equations (8) herein. A fixed phase shift may be to rotate the grid points to a mid-point between points in the preceding determination or in a preceding two determinations. A variable phase shift may be determined based on solutions in preceding determinations. In one embodiment a fixed phase shift is used between an initial determination and a second determination and a variable phase shift is used for a third determination.

[0115] A measure of relative magnitude between a first identified grid point and a second identified grid point from CS for the first and second grids respectively is used to determine an angular discriminant. The identified grid points indicate a resultant DOA estimate for the first and second grids, for example by being the grid point identified by a maximum value in an output vector from CS. For example, disparate magnitudes indicate an estimated DOA off-centre to the first and second grid points, whereas close or equal magnitudes indicate an estimated DOA at or near the centre of the identified grid points. The angular discriminant is used to determine an estimated direction of arrival. In some embodiments the estimated direction of arrival is determined directly from the angular discriminant and a preceding CS determination. In some embodiments the angular discriminant is used to identify a third grid incorporating a phase shift from the grid used for the preceding CS determination. In some embodiments, information identifying or related to the angular discriminant is output as an indication of a measure of error. In some embodiments the angular discriminant is used as a stop condition for an iterative process to arrive at a DOA estimate.

[0116] In some embodiments CS based on the first and second grids, as described above follow CS based on a third grid, whereby the first grid is rotationally offset from the third grid in one direction and the second grid is rotational offset from the third grid in the opposite direction (e.g., clockwise and anti-clockwise respectively). In one embodiment, the grid points of the first and second grids are offset to a mid-point between the grid points of the third grid. In some embodiments rotational offset is iteratively performed until a stop condition is met. In some embodiments the stop condition is the reaching of a threshold. The threshold selected may depend on the application. Examples of a threshold that may be selected include a change in DOA estimate between iterations is less than a predetermined value, a certain number of iterations have been completed, a certain amount of time has elapsed, or the angular discriminant is below a certain level. The threshold may be complex, for example continue until one of a plurality of conditions are met, continue until all of a plurality of conditions are met, or continue until a plurality of conditions are met or one or more other conditions are met.

[0117] In some embodiments the iterative procedure is applied to determine the maximum signal amplitude at each of the new grid point sets. These maximum amplitudes are then used in an angular discriminant (e.g., as described herein with reference to Equation (8)), which creates an estimate of the angular estimation error. This estimated angular error is then used to control the magnitude of the offset of the grid points.

[0118] In some embodiments the antenna elements are configured to optimise the results of CS. Optimisation variables may include one or both of the number of antenna elements and antenna element placement.

[0119] For example, in some embodiments the distance between the antenna elements is configured to optimise one or more parameters that affect at least one of the accuracy and convergence rate of CS and iterative CS. In the case of a UCA, a measure of distance between the antenna elements is the radius (equivalently the diameter) of the array. In some embodiments the one or more parameters include or consist of one or both of mutual coherence and condition number of a matrix of actual or simulated measurements from the array of antenna elements for one or more representative signals (e.g., signals of a particular narrow-band wavelength) of interest for DOA estimation, whereby optimisation includes minimising or at least reducing the mutual coherence and/or condition number. In some embodiments, the optimisation is constrained, for example having regard to physical size or shape constraints for the antenna array. In some embodiments, optimisation comprises locating the local minimum of mutual coherence or condition number that is associated with the smallest value for the condition number or mutual coherence respectively.

[0120] In some embodiments the antenna elements are odd in number, which at least in some implementations results in a performance improvement relative to an even number of antenna elements. The methods described are however also applicable to an even number of antenna elements. In particular implementations there are an odd number of elements of at least five elements, or at least nine elements. In some embodiments an odd number of antenna elements, particularly at least five, at least seven or at least nine, are used in combination with an optimised distance between antenna elements, for example an optimised radius in a UCA, as described above. In general, a lower number of antenna elements results in reduced computational complexity, but with a potential cost in accuracy.

[0121] In some embodiments the measured characteristics of the signal, for example from signal complex voltages represented as a set of complex envelope voltage outputs, obtained at the outputs of antenna elements configured in an array, is a single measurement value (e.g., a single vector of dimension M, representing the voltages from each antenna element). For example, in one embodiment measurements are obtained from all antenna elements at substantially the same time, so as to obtain a single time sample for use in CS. In another embodiment the sampling is repeated a plurality of times and the results are combined, for example by averaging, into a single measurement value. Using an average may statistically improve the accuracy of the DOA estimate. Using a single measurement value reduces the computational complexity for CS.

[0122] From the foregoing description, it can be appreciated that the problem of determining the angle of arrival is represented by a relationship, represented by an under-determined set of equations, which map signals origination at each grid point to the set of complex envelope voltage outputs of the antenna elements. Accordingly, the relations from the received signal from the possible grid points to the voltages at the antenna outputs can be represented by a system of equations. An example discussion of this issue is now provided.

Problem Formulation

[0123] A typical realization of the angle of arrival problem is with respect to a planner array of M isotropic omnidirectional elements equally distributed along a circular ring of Uniform Circular Array (UCA) with radius r and angular separation of

[00011]$\frac{2\pi }{M}$

radians. The inter-element spacing

[00012]$d=2r\sin \left(\frac{\pi }{M}\right)$

is the length of the straight line between two adjacent antenna elements. The angular positions of the antenna elements in UCA are represented by {γ.sub.m} where

[00013]${\gamma }_m=\frac{2\pi \left(m-1\right)}{M}.$

An electromagnetic plane wave impinges on the antenna elements from some unknown DOA θ.sub.p. The incident signal is considered to be narrow-band and characterized by the same frequency content. The narrow-band assumption states that all frequencies in the observed band have the same phase shift. This simplifies the construction of steering vector given in equation (2). Under the following assumption, the output of m.sup.th antenna array can be written as,

[00014]$\begin{array}{cc}{v}_m=\underset{i=1}{\overset{P}{.Math.}}{s}_i^{\mathrm{inc}}{\tau }_m\left({\theta }_i\right)\mathrm{where}{t}_m\left({\theta }_i\right)=\text{?}\text{?}\text{indicates text missing or illegible when filed}& \left(1A\right)\end{array}$

[0124] and [0125] s.sub.i.sup.inc represents the ith impinging wave, [0126] b is the angular wavenumber (2π/λ), [0127] r is the radius of UCA, [0128] θ.sub.p in the angle of arrival of pth incident wave, [0129] γ.sub.m is the angular position of mth element, [0130] λ is the wavelength of the wave,

[00015]$\left(2A\right)$

[0131] and, P is the number of impinging waves.

Or in the case of three dimensions:

[00016]$\begin{array}{cc}{v}_m=\underset{i=1}{\overset{P}{.Math.}}{s}_i^{\mathrm{inc}}{\tau }_m\left({\theta }_i,{\psi }_i\right)& \left(1\right)\end{array}$

where,

[00017]$\begin{array}{cc}{\tau }_m\left({\theta }_i,{\psi }_i\right)=e^{-\mathrm{jbrcos}\left({\theta }_i-{\gamma }_m\right)\cos \left({\psi }_i\right)}& \left(1\right)\end{array}$

and,

[0132] s.sub.i.sup.inc represents the i.sup.th impinging wave

[0133] b is the angular wavenumber

[00018]$\left(\frac{2\pi }{\lambda }\right)$

[0134] r is the radius of the UCA

[0135] θ.sub.i is the aximuth angle of the i.sup.th impinging wave

[0136] Ψ.sub.i is the elevation angle of the i.sup.th impinging wave

[0137] γ.sub.m is the angular position of the m.sup.th element

[00019]$\begin{array}{cc}{v}_m=\underset{i=1}{\overset{P}{.Math.}}{s}_i^{\mathrm{inc}}{\tau }_m\left({\theta }_i,{\psi }_i\right)\mathrm{where},& \left(1\right)\\ {\tau }_m\left({\theta }_i,{\psi }_i\right)=e^{-\mathrm{jbrcos}\left({\theta }_i-{\gamma }_m\right)\cos \left({\psi }_i\right)}& \left(2\right)\end{array}$

[0138] and, [0139] s.sub.i.sup.inc represents the i.sup.th impinging wave [0140] b is the angular wavenumber

[00020]$\left(\frac{2\pi }{\lambda }\right)$ [0141] r is the radius of the UCA [0142] θ.sub.i is the azimuth angle of the i.sup.th impinging wave [0143] Ψ.sub.i is the elevation angle of the i.sup.th impinging wave [0144] γ.sub.m is the angular position of the m.sup.th element [0145] λ is the wavelength of the wave

[0146] The phase of τ.sub.m is the phase shift due to the increased travel distance of the incoming signal in reference to the first element, while it is being received by the m.sup.th element of UCA.

DOA Estimation Using Compressive Sensing

[0147] Received open-circuit voltage information at each antenna element is combined to formulate a sparse matrix problem, which may be solved using CS techniques to identify the DOA of an unknown target. To incorporate the architecture of CS into the system model, angular space is being quantized into N discrete regions, each region having a representation value, or grid point. SF discretises the entire 2π radian angular domain into N possible DOAs, Θ={{circumflex over (θ)}.sub.n, 1≤n≤N} where N denotes the number of grid points. The incoming DOA, θ can be anywhere in the range [−π,π). Using the relationship between output of each antenna elements and DOA of the target in equation (1) and rewriting it in matrix form gives:

[00021]$\begin{array}{cc}V=\Phi \left(\Theta \right)S\mathrm{where},& \left(3\right)\\ \Phi \left(\Theta \right)=\left(\begin{array}{cccc}{T}_1\left({\hat{\theta }}_1\right)& T_1\left({\hat{\theta }}_2\right)& \mathrm{.Math.}& T_1\left({\hat{\theta }}_N\right)\\ T_2\left({\hat{\theta }}_1\right)& T_2\left({\hat{\theta }}_2\right)& \mathrm{.Math.}& T_2\left({\hat{\theta }}_N\right)\\ \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ T_M\left({\hat{\theta }}_1\right)& T_M\left({\hat{\theta }}_2\right)& \mathrm{.Math.}& T_M\left({\hat{\theta }}_N\right)\end{array}\right)& \left(4\right)\end{array}$

[0148] and, V={ν.sub.m, 1≤m≤M} is a one-dimensional column vector representing the complex output at each antenna elements of UCA. Φ(Θ) is the dictionary matrix, where τ.sub.m ({circumflex over (θ)}.sub.n) is calculated using equation (2). In general, Φ(Θ) can be computed for any antenna array geometry, (in 2 or 3-dimensional space) and any set of grid points (in 2 or 3-dimensional space). The iterative algorithm is in this sense universally applicable.

[0149] The column Φ(Θ) represents an M-element array response vector, for an incoming plane wave arriving from the direction {circumflex over (θ)}.sub.n. The vector S={s.sub.n.sup.inc, 1≤n≤N} is a one-dimensional column vector of size N, where s.sub.n.sup.inc represents the incoming plane wave from the direction {circumflex over (θ)}.sub.n. The outputs of antenna elements in equation (3) will be corrupted with a noise vector η.sub.M×1. The entries of η.sub.M×1 are statistically independent and are extracted from a complex Gaussian distribution with zero mean and variance σ.sup.2. The effect of noise on the output observation can be expressed as.

V.sub.n=Φ(Θ)S+η  (5)

[0150] The system defined in (5) is an under-determined set of equations, where M<N, and can be formulated as a CS problem to recover an estimate Ŝ of the original sparse vector S via convex optimization as shown in (1). Therefore

[00022]$\begin{array}{cc}\hat{S}=\underset{SϵC^N}{\min }{.Math.S.Math.}_{0}s.t.{.Math.V_n-\Phi \left(\Theta \right)S.Math.}_2<\mathrm{.Math.}& \left(6\right)\end{array}$

[0151] where is the ∥•∥.sub.p is the l.sub.p norm and ϵ is the regularization parameter that is being determined by the noise or quantization level. Since the model assumes a single transmitting source among the N possible DOAs, the recovered sparse vector will have only one nonzero element. The index n of the non-zero element refers to the angular grid) ({circumflex over (θ)}.sub.n) corresponding to the source DOA. Additionally, ϵ can be increased to cater for optimized designed for small antenna geometries, with close spacing between the antenna elements.

Multi-Resolution Approach

[0152] An assumption that the source is located on one of the angular grid points may not provide sufficiently accurate DOA estimation for at least some applications. In an ideal situation under this assumption, when the DOA of the source is on the grid, i.e. θ∈: {{circumflex over (θ)}.sub.n, 1≤n≤N}, the sparse vector solution discussed above enables detection of the DOA of the source. However, in a typical scenario when the source DOA is off the grid, i.e., θ={circumflex over (θ)}.sub.n+Δθ, where

[00023]$-\frac{\pi }{N}\le \Delta \theta \le \frac{\pi }{N},$

the DOA estimate includes an error. The dictionary mismatch between processed observation V.sub.n and measurement matrix Φ, forces the optimized solution vector S to generate several peaks at neighbouring grid points. In such cases, for an off-grid DOAθ, CS generates peaks at {circumflex over (θ)}.sub.k and {circumflex over (θ)}.sub.k+1, which are the neighbouring angular grids closest to the original off-grid DOAθ. To address this, the estimation process includes a two-stage strategy, wherein at the first stage, an index corresponding to the maximum amplitude is chosen as a coarse estimation k . The coarse estimate of θ, θ̆.sub.0, may be obtained from,

[00024]$\begin{array}{cc}{n}_{\max }={\max }^{-1}\left\{.Math.{\hat{S}}_0\left[n\right].Math.\right\}\mathrm{such}\mathrm{that},{\hat{S}}_0=\underset{Sϵ{•}^N}{\min }{.Math.S.Math.}_1\mathrm{st}{.Math.V_n-\Phi \left({\Theta }_0\right)S.Math.}_2<\mathrm{.Math.}\mathrm{and}{\stackrel{.Math.}{\theta }}_0={\Theta }_0\left(n_{\max }\right)& \left(7\right)\end{array}$

[0153] where, Ŝ.sub.0[n] is the n.sup.th element of the recovered sparse vector after CS processing, and Θ.sub.0 represents the set of N discrete azimuth angular grid points.

[0154] Assuming the Signal to Noise Ratio (SNR) is relatively high, there is a high probability that

[00025]$\theta \in \left[{\stackrel{.Math.}{\theta }}_0-\frac{\omega }{2},{\stackrel{.Math.}{\theta }}_0+\frac{\omega }{2}\right)$$\mathrm{where}$$\omega =\frac{2\pi }{N}$

is the angular grid separation. The determination of the course estimate in (7) is followed by an iterative process in the second stage. In the second stage, new grid points are determined, which usefully may be located at ½ the distance between the grid points used for the previous estimated angle of arrival. In one embodiment, two modified sparse vectors Ŝ.sub.t,1 and Ŝ.sub.t,2 are recovered by introducing a grid shift of

[00026]$-\frac{\omega }{2}$$\mathrm{and}$$\frac{\omega }{2},$

on the N possible DOAs Θ. In another embodiment, a single new modified vector Ŝ.sub.t is obtained due to the angular rotation of the grid points by

[00027]$\frac{\omega }{2}$

in one direction. t=1,2,3, . . . is the iteration number. In each iteration the magnitudes of the recovered sparse vector are used to determine a correction factor, which enables the algorithm to converge to an accurate estimate.

[0155] The coarse estimate is obtained by obtaining a compressive sensing solution of:

[0156] V.sub.n=Φ(Θ.sub.0)Ŝ.sub.o.

[0157] In particular, using compressive sensing, Ŝ.sub.0 is obtained by:

[0158] Ŝ.sub.0=min∥S.sub.0∥.sub.1, such that ∥V.sub.n−Φ(Θ.sub.o) S.sub.o∥.sub.2<ϵ,

[0159] where ∥•∥.sub.p is the l.sub.p norm of a vector.

[0160] Then, n.sub.max=max.sup.−1[|Ŝ.sub.o[n]|: 1≤n≤N] and θ̆.sub.o=θ(n.sub.max),

[0161] where θ̆ is a coarse (quantized) estimate of the angle of arrival of the signal.

[0162] More accurate estimates of the angle of arrival than the initial coarse estimate are obtainable through creating additive correction terms and applying these correction terms to the coarse estimate of the angle of arrival. The determination of the correction terms and the more accurate estimate of the angle of arrival are done in an iterative computational procedure.

[0163] For the first step of the iteration, define,

[00028]$\Delta {\Theta }_1=D\left({\alpha }_1,{\beta }_1\right)=\left(\frac{{\alpha }_1-{\beta }_1}{{\alpha }_1+{\beta }_1}\right)\left(\frac{\omega }{2}\right),\mathrm{where}$$V_n=\Phi \left[{\left({\Theta }_0+\frac{\omega }{2}\stackrel{\_}{u}\right)}_{2\pi }\right]{\hat{S}}_1$

[0164] where,

[0165] (Q).sub.2π=modulo(Q+πū, 2πū)−πū for any vector Q,

[0166] ū is the N×1 vector for which every element is 1.

[0167] ΔΘ.sub.1 is the correction factor for the first iteration,

[0168] D(α.sub.1,β.sub.1) is an angle of arrival error discriminant based on parameter defined in the following.

[0169] Ŝ.sub.1′ is solved by using compressive sensing.

Therefore Ŝ.sub.1=min∥S.sub.1∥.sub.1 such that

[00029]${.Math.V_n=\Phi \left[{\left({\Theta }_0+\frac{\omega }{2}\stackrel{\_}{u}\right)}_{2\pi }\right]{\hat{S}}_1.Math.}_2<\mathrm{.Math.}$$n_{\max }={\max }^{-1}\left[.Math.{\hat{S}}_0\left[n\right].Math.\text{:}1\le n\le N\right],{\alpha }_1=.Math.{\hat{S}}_1\left(n_{\max }\right).Math.,{\beta }_1=.Math.{\hat{S}}_1\left(k_{\max }\right).Math.,\mathrm{where},k_{\max }=\{\begin{array}{c}{n}_{\max }-1,\mathrm{for}2\le n_{\max }\le N,\\ N,\mathrm{for}{n}_{\max }=1.\end{array}$

[0170] A correction factor for the phase estimate is computed using the error discriminant,

[00030]$\Delta {\Theta }_1=D\left({\alpha }_1,{\beta }_1\right)=\left(\frac{{\alpha }_1-{\beta }_1}{{\alpha }_1+{\beta }_1}\right)\left(\frac{\omega }{2}\right)$

All of the grid points are rotated by ΔΘ.sub.1 which is represented as, Θ.sub.1=(Θ.sub.0+ΔΘ.sub.1ū).sub.2π
and therefore the estimate of the angle of arrival after the t=1 iteration is, θ̆.sub.1=Θ.sub.1 (n.sub.max)=Θ.sub.0(n.sub.max)+ΔΘ.sub.1=θ̆.sub.0+ΔΘ.sub.1
where ΔΘ.sub.1 is a correction factor for the coarse estimate θ̆.sub.0.
For the second iteration, if required, define

[00031]$\Delta {\Theta }_2=D\left({\alpha }_2,{\beta }_2\right)=\left(\frac{{\alpha }_2-{\beta }_2}{{\alpha }_2+{\beta }_2}\right)\left(\frac{\omega }{2}\right)$$\mathrm{where}$$V_n=\Phi \left[{\left({\Theta }_1+\frac{\omega }{2}\stackrel{\_}{u}\right)}_{2\pi }\right]{\hat{S}}_2$

Ŝ.sub.1′ is solved using compressive sensing, Ŝ.sub.2=min∥S.sub.1∥.sub.1 such that

[00032]${.Math.V_n-{\Phi \left({\Theta }_1+\frac{\omega }{2}\stackrel{\_}{u}\right)}_{2\pi }{S}_2.Math.}_2<\mathrm{.Math.}$$n_{\mathrm{ma}x}={\max }^{-1}\left[.Math.{\hat{S}}_o\left[n\right].Math.\text{:}1\le n\le N\right],{\alpha }_2=.Math.{\hat{S}}_2\left(n_{\mathrm{ma}x}\right).Math.,{\beta }_2=.Math.{\hat{S}}_2\left(k_{\mathrm{ma}x}\right).Math.,\mathrm{where},k_{\mathrm{ma}x}=\{\begin{array}{cc}{n}_{\mathrm{ma}x}-1,& \mathrm{for}2\le n_{\mathrm{ma}x}\le N,\\ N,& \mathrm{for}{n}_{\mathrm{ma}x}=1\end{array}.$

A correction factor for the phase estimate is computed using the error discriminant:

[00033]${\mathrm{\Delta \Theta }}_2=D\left({\alpha }_2,{\beta }_2\right)=\left(\frac{{\alpha }_2-{\beta }_2}{{\alpha }_2+{\beta }_2}\right)\left(\frac{\omega }{2}\right).$

All of the grid points are rotated by ΔΘ.sub.2, which is represented as: Θ.sub.2=(Θ.sub.1+ΔΘ.sub.2ū).sub.2π=(Θ.sub.0+ΔΘ.sub.1ū+ΔΘ.sub.2ū).sub.2π.
Therefore, the estimate of the angle of arrival after the t=2 iteration is: .sub.2=Θ.sub.2(n.sub.max)=Θ.sub.1(n.sub.max)+ΔΘ.sub.2=.sub.0+ΔΘ.sub.1+ΔΘ.sub.2,
where ΔΘ.sub.2 is the additional correction factor on the coarse estimate.

[0171] The iterative algorithm almost always converges after the second (t=2) iteration. Accordingly, in some embodiments the number of iterations is fixed at 2. In other embodiments, more than two iterations may be used.

[0172] For iteration t:

[00034]${\mathrm{\Delta \Theta }}_t=D\left({\alpha }_t,{\beta }_t\right)=\left(\frac{{\alpha }_t-{\beta }_t}{{\alpha }_t+{\beta }_t}\right)\left(\frac{\omega }{2}\right),\mathrm{where}{V}_n=\Phi \left[{\left({\Theta }_{t-1}+\frac{\omega }{2}\stackrel{\_}{u}\right)}_{2\pi }\right]{\hat{S}}_t.$

[0173] Ŝ.sub.t′ is solved by using compressive sensing,

[0174] Ŝ.sub.t=min∥S.sub.t∥.sub.1 such that

[00035]${.Math.V_n-{\Phi \left({\Theta }_{t-1}+\frac{\omega }{2}\stackrel{\_}{u}\right)}_{2\pi }{S}_t.Math.}_2<\mathrm{.Math.}$$n_{\mathrm{ma}x}={\max }^{-1}\left[.Math.{\hat{S}}_o\left[n\right].Math.\text{:}1\le n\le N\right],{\alpha }_t=.Math.{\hat{S}}_t\left(n_{\mathrm{ma}x}\right).Math.,{\beta }_t=.Math.{\hat{S}}_t\left(k_{\mathrm{ma}x}\right).Math.,\mathrm{where},k_{\mathrm{ma}x}=\{\begin{array}{cc}{n}_{\mathrm{ma}x}-1,& \mathrm{for}2\le n_{\mathrm{ma}x}\le N,\\ N,& \mathrm{for}{n}_{\mathrm{ma}x}=1\end{array}.$

[0175] A correction factor for the phase estimate is computed using the error discriminant,

[00036]${\mathrm{\Delta \Theta }}_t=D\left({\alpha }_t,{\beta }_t\right)=\left(\frac{{\alpha }_t-{\beta }_t}{{\alpha }_t+{\beta }_t}\right)\left(\frac{\omega }{2}\right)$

[0176] All of the grid points are rotated by ΔΘ.sub.t which is represented as,

[00037]$\begin{array}{cc}{\Theta }_t={\left({\Theta }_{t-1}+{\mathrm{\Delta \Theta }}_t\stackrel{\_}{u}\right)}_{2\pi }={\left({\Theta }_0+\left(\underset{k=1}{\overset{t}{.Math.}}{\mathrm{\Delta \Theta }}_k\right)\stackrel{\_}{u}\right)}_{2\pi }& \left(8\right)\end{array}$

[0177] and therefore the estimate of the angle of arrival after

[0178] the iteration t is,

[00038]${\stackrel{(}{\theta }}_t={\Theta }_t\left(n_{\mathrm{ma}x}\right)={\Theta }_{t-1}\left(n_{\mathrm{ma}x}\right)+{\mathrm{\Delta \Theta }}_t={\stackrel{(}{\theta }}_0+\left(\underset{k=1}{\overset{t}{.Math.}}{\mathrm{\Delta \Theta }}_k\right),$

[0179] where, ΔΘ.sub.t is the correction factor for the coarse estimate after iteration t.

The variables in the description above can be described as follows: [0180] Ŝ.sub.t is the N×1 minimum sparsity solution angle of arrival indicator vector for iteration t, [0181] V.sub.n is the M×1 vector of complex envelope voltages at each of the M antenna elements, [0182] Φ is the M×N observation matrix which depends upon the geometry of the antenna locations in the array and the positions of the angular grid points, [0183] α.sub.t is the magnitude of the solution vector element corresponding to a counter-clockwise shift of the previous angle of arrival estimation by ½ a quantization interval, [0184] β.sub.t is the magnitude of the solution vector element corresponding to a clockwise shift of the previous angle of arrival estimation by ½ a quantization interval,

[00039]$\omega =\frac{2\pi }{N}$

is the angular quantization step size,

[00040]$D\left({\alpha }_t,{\beta }_t\right)=\left[\frac{{\alpha }_t-{\beta }_t}{{\alpha }_t+{\beta }_t}\right]\frac{\omega }{2}$

is the error discriminant for iteration t,

[00041]${\mathrm{\Delta \Theta }}_t=\left[\frac{{\alpha }_t-{\beta }_t}{{\alpha }_t+{\beta }_t}\right]\frac{\omega }{2}$

is the scalar output of the error discriminant for iteration t [0185] Θ.sub.t is the vector of N grid points for iteration t, [0186] Θ.sub.t (n.sub.max) is the angle of arrival estimate for iteration t.
In the preceding description, a vector Q, (Q).sub.2π is defined as

(Q).sub.2π=modulo(Q+πū, 2πū)−πū.

[0187] This equation describes the 2π modulo operation on each of the elements in Q. By adopting this, a rotation in only one direction (either direction) is required to determine the angular discriminant. In alternative embodiments a rotation of

[00042]$-\frac{\omega }{2}\stackrel{\_}{u}$

is used in addition to the rotation

[00043]$+\frac{\omega }{2}\stackrel{\_}{u},$

with a corresponding increase in computational resources. In the case of a non-uniform grid of points, ω may differ for each direction, to maintain the rotated position at a half interval to either side of the course estimate.

[0188] In an alternative embodiment, an angular discriminant, for example the angular discriminant D(α.sub.t,β.sub.t) described above, is determined based on a first or initial compressive sensing, for example V.sub.n=Φ(Θ.sub.0)Ŝ.sub.o as described above. In the solution vector Ŝ.sub.o a maximum term is identified, and the maximum adjacent term is also identified. The angular discriminant is then computed based on the identified maximum term and maximum adjacent term. The grid points are then rotated by the angular discriminant for a first iteration. Subsequent iterations, if any, are performed as described above. If the adjacent terms on both sides of the maximum term are equal, then the estimated angle of arrival is the angle corresponding to the grid point for the maximum.

[0189] If there are two equal sized adjacent maxima in a compressive sensing solution, then the estimated angle of arrival is determined as the mid-point of the grid points corresponding to the two maxima.

[0190] The initial estimate of the angle of arrival is the angle corresponding to the mid-point between the angles corresponding to the two largest magnitude values of This midpoint is the initial angle of arrival estimate θ̆.sub.o [0191] α.sub.1 and β.sub.1 are the two largest magnitudes of adjacent components of the Ŝ.sub.o vector. [0192] α.sub.1 and β.sub.1 are then the input to D(α.sub.1,β.sub.1) and the output of this discriminant is the correction factor ΔΘ.sub.1. The angle of arrival estimate after the first iteration is θ̆.sub.1=θ̆.sub.0+ΔΘ.sub.1.

[0193] The grid points are then rotated by Θ.sub.1=(Θ.sub.0+ΔΘ.sub.1ū).sub.2π

[0194] The alternative embodiment of the algorithm has reduced computational complexity. One compressive sensing solution and one computation of the Φ matrix are eliminated by the alternative embodiment.

[0195] In some embodiments, the stopping criterion of the algorithm is determined by a user-defined threshold, for example for a fixed number of iterations or for D(α.sub.t,β.sub.t)<Ω.

[0196] The algorithm terminates at iteration t and .sub.t=Θ.sub.t(n.sub.max).

[0197] In some embodiments, a compressive sensing algorithm is used that accommodates pre-determination of the sparsity of the solution vector. For one transmitting source, the sparsity is set to 2, because there will be two significant adjacent elements in the solution vector. In the case of P signal sources, the sparsity is set to 2P . The CoSaMP algorithm is an example algorithm that accommodates predetermination of the sparsity of the solution vector and as such has particular application to the disclosed methods of direction of arrival estimation.

Compressive Sensing

[0198] An example discussion of CS is now provided.

[0199] Compressive sensing is a mathematical framework that deals with the recovery of a sparse vector x.sub.n×1, from an observation vector y.sub.n×1 with M<<N. The measurement paradigm consists of linear projection of the signal vector via a known projection matrix Ψ.sub.M×N. As M<<N, the recovery of sparse vector x from the measurement vector y becomes an underdetermined problem with an infinite number of solutions. In a CS framework, an accurate estimation of a sparse signal x can be obtained in the following reconstruction problem:

min∥x∥.sub.1 s.t. ∥y−Ψx∥.sub.2≤ζ,   (9)

[0200] where ∥•∥.sub.P is the l.sub.p-norm and ζ bounds the amount of noise in the observation data. A vector x is said to be K-sparse, if ∥x∥.sub.0=K. A matrix Ψ is said to have satisfied the RIP (Restricted Isometry Property) of order K, if there exists a δ.sub.k∈(0, 1) such that

(1−δ.sub.k)∥x∥.sub.2.sup.2≤∥Ψx∥.sub.2.sup.2≤(1+δ.sub.k)∥x∥.sub.2.sup.2   (10)

[0201] If ψ satisfied the above condition, there is a high probability of successfully recovering a sparse signal from noisy measurements, as long as the spark(ψ)>2K. The spark of a matrix is the smallest number of columns in matrix Ψ that are linearly independent. The larger the spark, the bigger the signal space, allowing CS to guarantee exact recovery. Although the spark and the RIP provides guarantees for the recovery of a k-sparse vector, verifying that a matrix satisfies any of the above properties has a combinatorial computation complexity, since each time one must consider

[00044]$\left(\begin{array}{c}n\\ k\end{array}\right)\hspace{1em}$

submatrices. Therefore, it is preferable to use a property of a matrix which is easily computable and provides guarantees of recovery.

[0202] The mutual coherence of a matrix ψ, μ(Ψ) is the largest absolute inner product between two columns Ψ.sub.i and Ψ.sub.j where

Ψ.sub.i is the i th column of Ψ and Ψ.sub.j is the j th column of Ψ.

[00045]$\begin{array}{cc}\mu \left(\Psi \right)=\underset{1<i\ne j<n}{\max }\frac{.Math.\left({\psi }_i,{\psi }_j\right).Math.}{{.Math.{\psi }_i.Math.}_2{.Math.{\psi }_j.Math.}_2}& \left(11\right)\end{array}$

[0203] The mutual coherence of a matrix Ψ is always bounded in the range

[00046]$\mu \left(\Psi \right)\in \left[\sqrt{\frac{N-M}{M\left(N-1\right)}},1\right],$

where the lower bound is known as the Welch Bound. Note that when N>>M, the lower bound is approximately equal to

[00047]$\frac{1}{\sqrt{M}}.$

If the original signal x satisfies the following requirements,

[00048]$\begin{array}{cc}{.Math.x.Math.}_0=K<\frac{1}{2}\left(1+\frac{1}{\mu \left(\Psi \right)}\right),& \left(12\right)\end{array}$

[0204] then, CS algorithms such as basis pursuit or other greedy algorithms such as COSAMP can be used to guarantee the recovery of x from under-determined set of equations.

[0205] A rectangular matrix such as Ψ.sub.M×N does not possess quantifiable parameters such as eigenvalues to determine the structure of the matrix. However, Q=Ψ.sup.TΨ can be considered as a square matrix and the eigenvalues of Q can be related back to quantify the property of Ψ. The singular values ρ.sub.1, . . . ρ.sub.m of a m×n matrix Ψ are the positive square roots, ρ.sub.i=√Λ.sub.i>0, of the non-zero eigenvalues of the associated Gram matrix Q=Ψ.sup.TΨ. Singular values of Ψ can be used to introduce another quantifiable parameter known as condition number, expressed as,

[00049]$\begin{array}{cc}{\rm Y}\left(\Phi \right)=\frac{{\rho }_{\max }}{{\rho }_{\min }}& \left(13\right)\end{array}$

[0206] where ρ.sub.min and ρ.sub.max are the smallest and largest singular value associated with the matrix Ψ. The condition number plays a vital role in providing a geometric interpretation of the action of the matrix. A matrix with lower condition number suggests strong convergence to an accurate and unique solution.

Array Geometry Optimization

[0207] In some embodiments the array of antenna elements is configured to optimise CS. Additionally or alternatively, the number of elements in the array may vary from a conventional approach of an inter-element separation,

[00050]$d\in \left[\frac{\lambda }{2},\lambda \right]$

between the antenna elements to avoid ambiguity between the steering vectors of distinct DOAs.

[0208] FIG. 3 shows a contour plot on the effect of antenna elements on μ(Φ). In the contour plot μ(Φ) is measured when constructed by varying M and the radius r in the range of [7,16] and [λ,10λ] respectively. The figure suggests that for an increase in odd M (7,9,11, . . . ) the reduction in mutual coherence is much sharper compared to even M (8,10,12, . . . ). Especially for M=9, Φ has a lower μ compared to M∈[8,10,12] and μ(Φ) achieves the minimum point for M=16. When the minimum μ(Φ) for M=11 is compared with M=12, it can be seen that, in case of UCA constructed with M=11 antenna elements, μ(Φ) is reduced by a factor of 10. The result clearly shows that the ambiguity between two distinct DOA can be significantly reduced in case of odd number of M>9.

[0209] Although

[00051]$d=\frac{\lambda }{2}$

has been used as an optimum separation to perform trade-off between mutual coupling and grating lobes, a geometry with

[00052]$d>\frac{\lambda }{2}$

or d>λ or d>2λ or d>3λ or d>4λ or d>5λ may be used with the DOA methods described herein. Also, a geometry with or d>1.5λ or d>2.5λ or d>3.5λ or d>4.5λ or d>5.5λ may be used with the DOA methods described herein. Referring to the example of UCA, FIG. 4 shows the variation in γ(φ) and μ(φ) with respect to varying radius of UCA for M=13 when r is varied between

[00053]$\frac{\lambda }{2}$

and 10λ with an increment of

[00054]$\frac{\lambda }{2}.$

In the plot the arrow marked r* refers to a radius of the array such that the inter-element spacing between the antenna elements is

[00055]$\frac{\lambda }{2}.$

It can be seen that at r* both γ(Φ) and μ(Φ) are significantly higher than at other points and hence is not optimal for an accurate recovery of using CS. An antenna array may be optimised by having a radius such that both γ(Φ) and μ(Φ) are minimized. An optimum radius r.sup.opt that contributes to the minimization of γ(Φ) and μ(Φ) may maximize the incoherence between the columns of Φ and efficient utilization of the vector space for CS operation. FIG. 3 indicates that at optimum radius r.sup.opt=6λ, γ(Φ) and μ(Φ) are reduced by a factor of 12 and 20 respectively compared to r*. FIG. 3 also indicates that the performance remains relatively steady between around 5λ to 9λ and at above around 9λ the performance starts to degrade. Accordingly, an antenna may be configured with d<9λ or d<10λ. A similar approach can be used to obtain the optimum radius for any UCA with M antenna elements.

[0210] In general, an optimum separation of antenna elements in the array is dependent on the number of antenna elements in the array. For example, for a UCA, a radius of about 6λ may be suited to about 9 to 17 antenna elements. A radius of about 8λ may be suited to 19 to 21 antenna elements. In some embodiments, the radius may be selected so that the number of antenna elements is within about 1.5 to 3 times the radius, or within 2 to 3 times the number of antenna elements.

[0211] In some embodiments, the distance between the antenna elements may be about 0.5λ. In other embodiments the distance between the antenna elements may be between 1λ and 2λ or between 1λ and 1.5λ or between 1.1λ and 1.5λ.

Phase Determination

[0212] In some embodiments, the phase information in V.sub.n described above, which is utilised for DOA estimation, is directly indicative of the relative phase between the complex envelopes received at the antenna elements. In other embodiments the phase information in V.sub.n is indicative of the phase of the complex envelope relative to a local oscillator. Using a local oscillator facilitates embodiments with larger signal to noise ratio.

Example Process Flow

[0213] FIG. 1 shows an example process flow for determining an estimated direction of arrival of a signal at an array of M antenna elements. An initial grid of N candidate directions of arrival is defined, with N greater than M to result in a sparse problem. The grid is rotated in both directions and an angular discriminant determined based on the rotated grid. The angular discriminant controls the extent of rotation of any subsequent iterations. An estimated direction of arrival is output when a threshold condition is met.

[0214] FIG. 2 shows another example process flow for determining an estimated direction of arrival of a signal at an array of M antenna elements. An initial grid of N candidate directions of arrival is defined, with N greater than M to result in a sparse problem. After an initial compressive sensing determination, the grid is rotated and the compressive sensing determination repeated. An angular discriminant is determined based on the solution vectors. The angular discriminant controls the extent of rotation of any subsequent iterations. An estimated direction of arrival is output when a threshold condition is met.

[0215] The process flows may be modified to enable DOA estimation in three-dimensional space.

[0216] In some embodiments, an estimated angle of arrival in three-dimensional space is determined based on individual determinations for transverse planes. For example, in some embodiments, the candidate directions of arrival are located in a first plane and the estimated direction of arrival is an estimated direction of arrival for that plane. To determine the DOA in three-dimensional space, the method further includes repeating the determinations in respect of candidate angles of arrival located in a second plane having at least a component substantially transverse to the first plane, to determine an estimated direction of arrival for the second plane. An estimated direction of arrival is then based on the estimated direction of arrival for the first and second planes.

[0217] In some embodiments, the second plane is perpendicular to the first plane. In some embodiments, the second plane intersects the first plane along a line having a direction corresponding to the first estimated direction of arrival. In some embodiments the method further comprises repeating the determinations in respect of a grid of candidate angles of arrival located in a third plane, the third plane intersecting points on a line in three-dimensional space corresponding to the second estimated direction of arrival. In some embodiments, the method comprises iteratively determining estimated directions of arrival in planes with substantial components transverse to the preceding plane until a threshold minimum variation in estimated direction of arrival is reached.

[0218] In some embodiments, the N potential angles of arrival are spatially separated in three-dimensional space, whereby Ŝ.sub.t for t=0 has solution vector elements for both azimuth and elevation. The method may then comprise applying at least the iterations t=1 and t=2 to determine the azimuth in relation to the largest absolute value adjacent pair of elements with constant elevation and applying at least the iterations t=1 and t=2 to determining the elevation in relation to the largest absolute value adjacent pair of elements with constant azimuth.

Simulation

[0219] A simulation was carried out on N=180 angular grid points, with

[00056]$\omega =\frac{2\pi }{N}.$

The scanning angle ranges between [−π, π) radians. The signal is considered to be transmitted at centre frequency of ƒ.sub.c MHz, and the wavelength is λ. A simulated UCA consists of 13 isotropic antenna elements distributed evenly on a circular ring with r=r.sup.opt=6λ. The inter-element distance d between the antenna elements is approximately 3λ. The simulation scenario has one source, transmitting from any angle in the range between [−π,π) radians. The signals have been supposed to be arriving on the antenna with equal strength in order to perform an unbiased analysis of the accuracy of the method with respect to the angles of arrival.

[0220] In order to determine the robustness of the system model, the following noise sensitivity test was considered. The Signal-to-Noise-Ratio (SNR) is calculated at the receiver as the ratio of the sum of the power received from M antenna elements to σ.sup.2 where, σ.sup.2 is the variance of the complex Gaussian noise. The measured data are characterized by SNR.sub.db=[−10, −5, 0, 5, 10, 15, 20, 25], defined as,

[00057]$\begin{array}{cc}{\mathrm{SNR}}_{\mathrm{dB}}=10{\log }_{10}\left[\frac{{.Math.}_{m=1}^M{.Math.v_m.Math.}^2}{M{\sigma }^2}\right]& \left(14\right)\end{array}$

[0221] where, ν.sub.m, m=1, . . . M, is the noiseless complex voltage observation at each antenna element. Since the actual DOA can be placed anywhere in the range [−π,π), T=1000 different scenarios were considered, to give a consistent statistical validation. Compressive Sampling Matching Pursuit (CoSaMP) performs the CS operation. The performance parameter of the algorithm is characterized as Mean Square Error (MSE), where MSE is defined as,

[00058]$\begin{array}{cc}\stackrel{32}{\mathrm{MSE}}=\frac{\underset{a=1}{\overset{T}{.Math.}}.Math.{\stackrel{\_}{\theta }}_{\mathrm{org},a}-{\stackrel{\_}{\theta }}_{\mathrm{est},a}.Math.}{T}& \left(15\right)\end{array}$

where, θ.sub.org,a is the original DOA of the source and θ.sub.est,a is the estimated DOA using the algorithm.

[0222] The MSE of the proposed algorithm is compared with the Cramer-Rao Lower Bound (CRB or CRLB), as

[00059]$\begin{array}{cc}\mathrm{CRLB}\ge \frac{-1}{\left(-\frac{2}{{\sigma }^2}\right)\left(M\frac{b^2{r}^2}{2}\right)}=\frac{{\sigma }^2}{{\mathrm{Mb}}^2{r}^2}& \left(16\right)\end{array}$

[0223] A set of results are presented in FIG. 5 to show the impact of external noise in estimating the original DOA of the source when an UCA is constructed with r=r* and r=r.sup.opt respectively (see FIG. 4). The results show that for a UCA geometry with r.sup.opt, the algorithm under the simulation conditions achieves the CRB for SNR.sub.dB>10 and remains on the bound for higher SNR. On the other hand, in case r*, the graph approaches the bound for SNR.sub.dB=3, but deviates away for higher SNR. The plot provides a clear indication that, UCA constructed with optimum radius (r.sup.opt) is more efficient in detecting the DOA of a source with minimum error.

[0224] Another set of simulations were carried out to examine convergence of the recursive algorithm in achieving the CRLB of DOA estimation. FIG. 6 shows plots that graphically display the convergence of the algorithm through a first estimate and through two iterations. The CRB is also plotted.

[0225] For SNR of 5 dB or less, the estimates have the same mean square error angle. Above 5 dB the first (course) estimate remains constant at a MSE about 10.sup.−4 whereas the first and second iterations perform closer to the CRB, the second iteration converging on the bound for about SNR>10.

[0226] COSAMP has a complexity of O(MN) in determining the solution of a sparse vector. The proposed method converges to the bound using just 2 iterations. Compared to Eigen-Value Decomposition (EVD) based DOA estimation such as (MUSIC and Root-MUSIC), the proposed algorithm therefore has much lower computational complexity.

[0227] Although the simulation was performed with N=180 for M=13, N may be increased or decreased. A reduction in N reduces computational complexity. For example, N may be reduced to approximately 100, or approximately 50, or approximately 40, or increased to approximately 250, 360 or more. In general, a minimum for N may be determined by the maintenance of sufficient sparsity for CS, which for some implementations may be between about two to three times M, whereas a maximum may be determined computational cost.

[0228] FIGS. 7 and 8 show flow diagrams of embodiments of a process to perform DOA. The process may be performed, for example, by a computational processor configured with instructions, for example instructions held in non-volatile memory. In some embodiments the computational processor is in communication with an antenna array, for example a UCA as described herein above. In such embodiments DOA estimation may be performed substantially in real-time.

[0229] In steps 100, 100A a set of complex envelope voltage outputs are received from the antenna element array. These may be stored in transient or non-transient memory for further processing. In step 101, 101A the set of complex envelope voltage outputs are transformed by an orthogonal transform to increase sparsity. Step 101, 101A is omitted in other embodiments. In step 102, 102A CS is applied to the transformed outputs and a grid including a higher number of grid points than measured outputs, to reveal a first DOA estimate, specified by one of the grid points (the DOA grid point). In step 103, 103A two new grids or one new grid is defined, rotated with respect to the first grid. For example, the new grids include grid points that are rotated a half grid quantization interval. In step 104, 104A CS is applied to the new grid(s) and an angular discriminant is determined based on preceding CS solutions. In step 105, 105A a decision is made whether a threshold condition, for example based on the angular discriminant, has been met. If so, the process ends and the latest DOA estimate is used as the final DOA estimate. If not, the process returns to step 103, 103A.

[0230] Accordingly, the solutions of the CS operations yield the magnitudes derived from the two shifted sets of grid points, respectively. The magnitudes of the shifted grid points closest in angle to the previous direction of arrival estimate are used as the input to a phase error discriminant. The output of the phase error discriminate is then used to adjust the estimate of the angle of arrival. This process is continued in an iterative manner until the output of the phase error discriminant is below an acceptable, user-defined threshold. On each iteration, the estimate of the angle of arrival improves, until there is negligible discriminate output.

[0231] FIG. 9 shows a block diagram of a radio receiver 200 for DOA estimation. The radio receiver includes an antenna element array 201, filter and amplifier circuitry 202, sampling circuity 203, which may include a local oscillator as a phase reference for the received signals, to sample the filtered and amplified signals from filter and amplifier circuitry 202 and a processor 204 for processing the signals, for example according to an embodiment described herein. The processes described herein above are suited for samples of a narrow-band signal. Where a DOA estimate is required for a wide-band signal, the filter and amplifier circuitry 202 includes a filter to compensate for phase shift in the received signal as a function of frequency. The processor 204 may hold the samples in transient memory 205 or process the signals according to instructions held in non-transient memory 206. The processor 204 provides the DOA estimate via an output 207, for example a display and/or a data communication interface to another device. The radio receiver 200 may receive a radio signal from a radio source 300, a plurality of radio signals from the radio source 300 and/or one or more radio signals from a plurality of radio sources 300 and determine DOA estimate(s) based on each of the received signals, or where appropriate based on a combined estimate for a plurality of signals.

[0232] FIGS. 10 and 11 show simulation results obtained for an embodiment of the disclosed methods. The simulations provide insight into the MSE performance of each of two example antenna array geometries under the influence of varying SNR. The MSEs of the estimation are compared to the theoretical CRLBs of the respective antenna geometries. The DOA of the incoming signal is determined for the disclosed embodiment (labelled “Proposed” in FIG. 10 and “ICSDOA” in FIG. 11). DOAs are randomly chosen, with an assumption that the DOAs are the in the range [−π, π). The number of Monte Carlo runs for each DOAs is set to T=5000. Two antenna geometries are considered, where both UCA and ULA are constructed with M=9 antenna elements with an inter-element separation of

[00060]$d_{\mathrm{UCA}}=d_{\mathrm{UCA}}=\frac{\lambda }{2}.$

The number of angular grid points for UCA and ULA is set to be N.sub.UCA=N.sub.ULA=180, with grid interval

[00061]${\omega }_{\mathrm{UCA}}=\frac{2\pi }{N_{\mathrm{UCA}}}\mathrm{and}{\omega }_{\mathrm{ULA}}=\frac{2\pi }{N_{\mathrm{ULA}}}$

respectively.

[0233] The CRLB for both UCA and ULA are shown. The CRLB of the ULA is lower than that of the UCA. The MSE plots for both the antenna geometries behave in a similar fashion, dipping off at an approximate SNR=6 dB and continuing to be on the CRLBs for higher SNR. For SNR<5 dB, the MSEs are relatively higher than the CRLBs with ULA having a lower MSE than UCA. The high MSE at low SNR regions can be associated with the inaccurate grid estimation of the disclosed algorithm, where the underlying CS operation fails to detect the angular grid on which the source is located.

[0234] Embodiments of the iterative compressive sensing direction of arrival estimation algorithm (ICSDOA) described above have significantly less computational complexity than previous algorithms that obtain estimates of the angle of arrival.

[0235] MUSIC is the Multiple Signal Classification Algorithm.

[0236] Root MUSIC is the Root Multiple Signal Classification Algorithm.

[0237] ESPIRIT is the Estimation Signal Parameter via a Rotation Invariant Technique.

[0238] Certain embodiments of the MUSIC Algorithm have computational complexity of O(PM.sup.2N+M.sup.2), where [0239] P is the number of time samples (or snapshots) of the signals at the antenna outputs, [0240] M is the number of antennas in the array, [0241] N is the number of elements in the quantized grid

[0242] Certain embodiments of the Root-MUSIC algorithm also have computational complexity O(PM.sup.2N+M.sup.2).

[0243] Certain embodiments of the ESPIRIT Algorithm have computational complexity of O(PM.sup.2+M.sup.3).

[0244] Certain embodiments of MUSIC, Root MUSIC, and ESPIRIT require P to be much greater than 1 for successful operation.

[0245] Certain embodiments of the ICSDOA iterative algorithm for the estimation of the angle of arrival has computational complexity of O(3MN), where O(3MN) is for compressive sensing and O(3MN) is for the evaluation of the dictionary matrix Φ. An alternative implementation of the iterative algorithm has computational complexity of O(4MN).

[0246] The ICSDOA iterative algorithm has much less computational complexity than at least certain embodiments of MUSIC, Root Music, or ESPIRIT. The iterative algorithm obtains an estimate with only one time sample from the antenna elements. A sequence of output estimates may be further operated upon, if required, with signal processing to produce a reduced error estimate of the angle of arrival.

[0247] Further aspects and embodiments of the present disclosure will be apparent from the following description, given by of example to a radio signal. In other example the embodiments are applied to a sound signal. In other examples the embodiments are applied to a phase coherent light signal.

[0248] A method for use in a direction of arrival estimation for a radio signal, includes: receiving, at a computational processor, a set of measurements of a radio signal from a radio source taken by an array of antenna elements; generating first and second measures of a direction of arrival estimate, the generating based on first and second grids of potential direction of arrivals respectively, the first and second grids offset from each other; generating an angular discriminant based on the first and second measures; generating third and fourth measures of a direction of arrival estimate, the generating based on third and fourth grids of potential direction of arrivals respectively, the third and fourth grids offset from the first and second measures by an amount based on the angular discriminant.

[0249] A method for use in a direction of arrival estimation for a radio signal, includes: receiving, at a computational processor, a set of measurements of a radio signal from a radio source taken by an array of antenna elements; generating, by the computational processor based on the received measurements, a first measure associated with a first direction of arrival estimate for the radio signal, based on a first grid with a plurality of grid points corresponding to potential directions of arrival, the grid comprising a larger number of grid points than antenna elements in the array of antenna elements and a lower resolution of grid points than required to achieve a target accuracy for the direction of arrival estimation; generating, by the computational processor, a second measure associated with a second direction of arrival estimate for the radio signal, based on a second grid comprising grid points around the first direction of arrival estimate that are offset to grid points in the first grid; and determining, by the computational processor, an angular discriminant based on the first measure and the second measure, wherein the measures associated with the direction of arrival estimates are based on a solution to a sparse problem defined by the received set of measurements and the respective grid points.

[0250] In some embodiments the set of measurements comprise a single measurement value.

[0251] In some embodiments the set of measurements comprises measurements from a circular array. In some implementations the circular array is a uniform circular array.

[0252] In some embodiments the set of measurements comprises measurements corresponding to an odd number of antenna elements.

[0253] In some embodiments the set of measurements comprises at least 5 measurements, or between 7 and 25 measurements, or between 9 and 23 measurements, or between 9 and 21 measurements, or between 9 and 19 measurements, or between 9 and 17 measurements, or between 9 and 15 measurements.

[0254] An antenna array for direction of arrival estimation includes a plurality of antenna elements arranged in a substantially uniform array, each antenna element configured to provide a measurement signal of a radio signal having a base wavelength, wherein a distance between pairs of antenna elements is substantially equal to a distance that minimises at least one of or a combined measure of a mutual coherence and a condition number of a matrix of said measurement signals of a radio signal by the antenna array at the base wavelength.

[0255] An antenna array for direction of arrival estimation includes a plurality of antenna elements arranged in a substantially uniform array, each antenna element configured to provide a measurement signal of a radio signal having a base wavelength, wherein a distance between pairs of antenna elements is greater than a distance corresponding to one wavelength at the base wavelength.

[0256] In some embodiments the distance between pairs of antenna elements in the antenna array is greater than a distance corresponding to two wavelengths at the base wavelength.

[0257] In some embodiments the distance between pairs of antenna element in the array is one half wavelength.

[0258] In some embodiment the distance between pairs of antenna elements is less than one wavelength.

[0259] In some embodiment the distance between pairs of antenna elements is less than one half wavelength.

[0260] In some embodiments the distance between pairs of antenna elements in the antenna array is less than a distance corresponding to ten wavelengths at the base wavelength. In some implementations the distance between pairs of antenna elements is about a distance corresponding to six wavelengths at the base wavelength.

[0261] In some embodiments the array comprises an odd number of antenna elements.

[0262] In some embodiments the number of antenna elements is at least 5 or at least 7 or at least 9.

[0263] In some embodiments the number of antenna elements less than or equal to 15. In other embodiments the number of antenna elements is more than 15.

[0264] In some embodiments the antenna array is substantially a uniform circular array.

[0265] In some embodiments the antenna array is substantially a uniform linear array.

[0266] In some embodiment the antenna array is of a geometry that is neither a uniform linear array nor a uniform circular array.

[0267] A radio receiver for direction of arrival estimation includes: an antenna array according to any embodiment described in the preceding paragraphs; and a computational processor configured to receive measurement signals for a radio signal from the antenna array and generate a direction of arrival estimation based on the measurement signals, the direction of arrival estimation utilising compressive sensing.

[0268] In some embodiments of the radio receiver, the computational processor is configured to perform the method as described in the preceding paragraphs.

[0269] The radio receiver of claim 17 or claim 18, configured to provide a direction of arrival estimate anywhere within the range 0 to 2π.

[0270] A method of direction of arrival estimation for a radio signal includes: receiving, at a computational processor, a set of measurements of a radio signal from a radio source taken by an array of antenna elements; generating, by the computational processor based on the received measurements, a direction of arrival estimate for the radio signal, wherein the direction of arrival estimate is based on compressive sensing of a sparse problem defined by a de-correlating transform of the received set of measurements and a grid with a plurality of grid points corresponding to potential directions of arrival, the grid comprising a larger number of grid points than antenna elements in the array of antenna elements.

[0271] In some embodiments the set of measurements of a radio signal from a radio source taken by an array of antenna elements, is a set of measurements from a circular array, which may be a uniform circular array.

[0272] In some embodiments the set of measurements consists of measurements corresponding to an odd number of antenna elements, for example between 5, 7 or 9 elements and 15 elements.

[0273] A method of direction of arrival estimation for a radio signal includes generating, by a computational processor based on received measurements of a radio signal from an antenna element array, a direction of arrival estimate for the radio signal from within a possible range of 0 to 2π, wherein the direction of arrival estimate is based on compressive sensing of a sparse problem defined by the received set of measurements and a grid with a plurality of grid points corresponding to potential directions of arrival, the grid comprising a larger number of grid points than antenna elements in the array of antenna elements.

[0274] It will be understood that the invention disclosed and defined in this specification extends to all alternative combinations of two or more of the individual features mentioned or evident from the text or drawings. All of these different combinations constitute various alternative aspects of the invention.