SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

Open Access Research

Cooperative DoA-only localization of primary users in cognitive radio networks

Federico Penna1* and Danijela Cabric2

Author Affiliations

1 Fraunhofer Heinrich Hertz Institute, Berlin D-10587, Germany

2 Electrical Engineering Department, University of California Los Angeles (UCLA), Los Angeles, CA 90095-1594, USA

For all author emails, please log on.

EURASIP Journal on Wireless Communications and Networking 2013, 2013:107  doi:10.1186/1687-1499-2013-107

The electronic version of this article is the complete one and can be found online at: http://jwcn.eurasipjournals.com/content/2013/1/107


Received:2 July 2012
Accepted:26 March 2013
Published:19 April 2013

© 2013 Penna and Cabric; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

In this paper, we investigate localization techniques based on direction-of-arrival (DoA) measurements to estimate the position of primary users in cognitive radio networks. In the proposed approach, different multi-antenna sensors estimate the DoA of a primary signal by subspace-based techniques and the target position is estimated from the available DoA measurements using maximum-likelihood, least-squares, or Stansfield estimators. The resulting localization performance is evaluated numerically and compared to the Cramr-Rao bound derived under a considered problem setting. The impact of several system parameters (in particular, number of sensors and number of antennas per sensors) is thoroughly analyzed and discussed.

Keywords:
Cognitive radio; Spectrum sensing; Localization; DoA estimation; Cramr-Rao bound

1 Introduction

After more than 10 years of research on cognitive radio (CR) systems [1,2], a large variety of spectrum sensing techniques have been proposed (see, e.g., [3-5]). Even the most sophisticated sensing techniques, however, fail to provide full protection to primary users (PUs) due to unavoidable detection delays (depending on the number of samples needed before a decision is made) and vulnerability to adverse radio conditions (e.g., fading, physical obstacles). This conclusion has been recently acknowledged by the Federal Communications Commission, whose latest memorandum [6] replaced spectrum sensing by a database-oriented approach such that PU positions and times of activity are known in advance.

Knowledge of PU position, in particular, provides important information for planning and management of a CR network: it can be exploited to identify ‘spectrum holes in space’ [7], to implement power control mechanisms, to design location-aware routing protocols, and to handle mobility of secondary users (SUs). Also, performance and reliability of a CR network can be greatly improved if SUs are able to localize PU signal sources besides detecting only their presence, like in traditional spectrum sensing. PU localization capabilities are particularly beneficial for cognitive networks operating in spatially diverse radio environments without databases.

The problem of PU localization in CR networks has attracted some attention over the past few years. For example, in [8], a Bayesian approach based on range and phase measurements is proposed, while a semi-range method is introduced in [9]. However, range-based approaches are difficult to apply in practice because the transmission power of PUs is typically unknown as well as the propagation channel between PUs and SUs. In addition, theoretical bounds and approximate maximum-likelihood algorithms based on received signal strength measurements were derived in [10], showing that the achievable estimation accuracy is relatively poor if all PU transmission parameters are unknown and SUs can rely only on energy measurements. Therefore, different approaches need to be explored for application in CR networks. Time-difference-of-arrival (TDoA) techniques could be adopted in principle, but they require perfect synchronization among different SU receivers, which is not always possible in CR systems. A weighted centroid approach was proposed in [11], which is suitable for dense CR networks with many sensors close to the PU position.

In this paper, we consider PU localization based on direction-of-arrival (DoA) measurements only. As an advantage compared to range-based localization, the DoA approach does not require knowledge of PU transmission power nor of the propagation channel. This localization technique, often referred to as bearings-only target location, was originally proposed and analyzed in the navigation literaturea (see, for example, [12-16]).

DoA estimates can be obtained in different ways, for example, by directional antennas or antenna arrays [17]. This work concentrates on the latter case. Note that arrays may be multiple antennas of one receiver or multiple users cooperating with each other (‘virtual arrays’). In the following, we will use the general term sensor to denote a sensing unit (single device or virtual array) composed of a number m > 1 of antennas. Each sensor can estimate the DoA of the PU signal using classic array processing techniques [18] such as the multiple signal classification (MUSIC) algorithm [19].

The contributions of this paper are as follows:

• Derivation of Cramér-Rao bound (CRB) of the achievable PU localization accuracy as a function of number of sensors in the network, number of antennas per sensor, distance of the sensors from the target, and array orientation with respect to the target.

• Comparison of the above bound with practical algorithms. It is shown that a combination of MUSIC and Stansfield algorithm [20] can nearly reach the CRB in conditions of practical interest for CR networks.

• Study of the impact of number of sensors and number of antennas on localization accuracy and numerical evaluation of the trade-off between these two factors.

We remark that results on the CRB of multi-sensor localization were already derived in the navigation literature (e.g., [13-16]) and later applied to wireless sensor networks (see, e.g., the survey paper [17]). However, our work is, to the best of our knowledge, the first to consider the joint effect of the number of sensors and the number of antennas per sensor in a multi-sensor, multi-antenna target localization system. The results of this paper show the impact of several parameters (number of sensor elements, choice of DoA fusion algorithm, array orientation, signal-to-noise ratio) on the overall localization performance. Applied to cognitive radio networks, these results provide important design guidelines for implementing practical PU localization systems. In this paper, as common in the localization literature, we adopt the position root mean square error (RMSE) as performance metric to evaluate the proposed algorithms. The RMSE has a theoretical justification, in that it allows for direct comparison with the CRB. In practical CR applications, an error in the estimation of PU position (expressed by the RMSE) leads to an erroneous identification of the boundaries of a spectrum hole in space and, ultimately, to a reduction of the throughput for the secondary network or to an increase of interference for the primary network. While the relation between PU localization accuracy and PU/SU communication performance is strongly application-dependent, minimizing the PU position RMSE translates in general into maximizing the efficiency of dynamic spectrum usage.

The paper is organized as follows: the problem model is introduced in section 2. In section 3, we briefly state preliminary results on the DoA estimation variance achieved by the MUSIC algorithm and its corresponding CRB. In section 4, we address the problem of localization by combining DoA estimates of multiple sensors: we first derive the CRB of the achievable localization error, highlighting its dependence on distances and angles between sensors and target, and then we introduce practical algorithms: maximum-likelihood (ML), least-squares (LS), and Stansfield estimators. In section 5, we evaluate by numerical simulations the performance of the proposed methods for DoA estimation and localization; we study the impact of various system parameters and investigate the trade-off between number of sensors and number of antennas per sensor. In section 6, we discuss implementation issues of the proposed method. Section 7 concludes the paper and summarizes the key results.

2 Mathematical model

Consider a CR network consisting of M sensors with m antennas each. Let n be the number of PUsb. The position of the ith PU (target) is denoted by <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M1">View MathML</a>, and the position of jth SU (assumed known) is <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M2">View MathML</a>. We then define the distance between PU i and sensor j as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M3">View MathML</a>

(1)

where ∥ · ∥ denotes Euclidean norm, and the DoA of the ith signal with respect to the jth sensor is given by

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M4">View MathML</a>

(2)

Angles θi,j can be estimated by each sensor j using subspace-based methods (such as MUSIC), which exploit the eigendecomposition of the m × m covariance matrix constructed from N signal samples received at the m antennas. A generic received signal sample <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M5">View MathML</a> at sensor j and time instant t can be expressed as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M6">View MathML</a>

(3)

where <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M7">View MathML</a> is a Gaussian noise vector with zero mean and covariance σ2Im, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M8">View MathML</a> is a vector of transmitted signal samples with zero mean and covariance P, and

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M9">View MathML</a>

(4)

is the array response matrix. If primary signals are uncorrelated, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M10">View MathML</a>. We then define the signal-to-noise ratio (SNR) of signal i as <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M11">View MathML</a>.

Remark 1. By convention, the DoA support is limited to one semi-plane, e.g.,θi,j∈ [− π /2, π /2] in Equation 2. This choice is derived from theπ-ambiguity of MUSIC and other subspace-based methods using uniform linear arrays with isotropic antennas. Such ambiguity can be resolved if some prior information about PU location is available (e.g., a coarse, initial estimate obtained by directional antennas) so that the reference semi-plane is known in advance.

Remark 2. Angles denoted asθi,j, e.g., in Equation 2, are defined with respect to the positive direction of thexaxis. The same angles, when denoted as<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M12">View MathML</a>, e.g, in Equation 3, are defined with respect to the specific reference system of thejth sensor. An example is shown in Figure1. It is assumed that each SU knows the relative orientation of its own array; therefore, angles<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M13">View MathML</a>can be converted intoθi,j.

thumbnailFigure 1. Example: relative positions of a sensor (<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M14">View MathML</a>)and of a primary user(<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M15">View MathML</a>). The angle θi,j is defined with respect to a common reference system (horizontal axis in the figure), while <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M16">View MathML</a> is defined with respect to the orthogonal to the positive direction of the sensor array axis.

3 DoA estimation error

3.1 Cramér-Rao bound

Given an unbiased estimator of <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M17">View MathML</a>, the estimation variance is lower-bounded by the CRB, denoted by <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M18">View MathML</a>. In [21], the CRB of DoA estimation is derived as a function of the SNR ρi, the number of antennas m of sensor j, and the number of received signal samples N. Under the assumption of uncorrelated signals (diagonal P), the result is as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M19">View MathML</a>

(5)

where

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M20">View MathML</a>

(6)

and <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M21">View MathML</a>. Note that A is a short-hand notation for <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M22">View MathML</a>.

Let us now consider the special case of uniform linear array (ULA) with single signal (n = 1). The ULA array response is given as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M23">View MathML</a>

(7)

where ι is the imaginary unit, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M24">View MathML</a> (with λ = signal wavelength), D is the spacing between array elements (which should satisfy D ≤ λ /2 to avoid spatial aliasing). By convention, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M25">View MathML</a> for ULAs is defined as the angle between the direction orthogonal to the array line and the impinging wave. In the presence of a single signal, the matrix A consists of a single column equal to <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M26">View MathML</a> (where index i = 1 can be dropped with no ambiguity). Therefore, Equation 6 reduces to the following:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M27">View MathML</a>

(8)

After simple algebraic manipulations (see also Sec. IV-A of [21]), we find that

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M28">View MathML</a>

(9)

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M29">View MathML</a>

(10)

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M30">View MathML</a>

(11)

We thus obtain the following:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M31">View MathML</a>

(12)

The factor <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M32">View MathML</a> can be interpreted as the effect of the orientation of the ULA with respect to the line between target signal and sensor. The DoA estimation CRB is minimized when the ULA is perpendicular to such line (<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M33">View MathML</a>), whereas it tends to when the ULA is aligned to the same direction of the received signal (<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M34">View MathML</a>).

3.2 Practical algorithms

Multi-antenna DoA estimation can be performed by classic array processing techniques such as MUSIC [19] or one of its numerous variants (see [18] for a survey).

In [21], it was shown that the signal parameter estimates obtained through MUSIC are, asymptotically in the sample size, unbiased and Gaussian distributed. Furthermore, the estimation variance achieved by MUSIC (denoted <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M35">View MathML</a>) can be expressed as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M36">View MathML</a>

(13)

which converges to the corresponding CRB <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M37">View MathML</a> (Equation 5) as the SNR ρi increases. The symbol [X]ij represents the element (i,j) of a matrix X.

Considering again the special case of ULA with single PU signal, we obtain the following:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M38">View MathML</a>

(14)

4 Localization error

4.1 Cramér-Rao bound

Assume now that one of the nodes in the network (the fusion center) gathers M DoA estimates referred to PU i, obtained by M different sensors in the network. The set of available observations, properly converted according to a common reference system (see section 2, Remark 2), is expressed by the following vector:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M39">View MathML</a>

(15)

DoA estimates can be written as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M40">View MathML</a>

(16)

where, recalling that the MUSIC estimator is asymptotically unbiased and Gaussian distributed (see section 3.2), error variables wi,j are modeled as zero-mean Gaussian random variables: <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M41">View MathML</a>.

The fusion center, then, computes an estimate of the position of the ith PU based on observations <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M42">View MathML</a>. If such estimate <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M43">View MathML</a> is unbiased, the estimation error covariance matrix is lower-bounded by the CRB:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M44">View MathML</a>

(17)

where F is the Fisher information matrix (FIM) given by

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M45">View MathML</a>

(18)

Thanks to the Gaussian distribution of wi,j, the FIM can be expressed as follows (cf. [16]):

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M46">View MathML</a>

(19)

where <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M47">View MathML</a> and <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M48">View MathML</a>. We now observe that Δxi,j = Ri,j cos(θi,j) and Δyi,j = Ri,j sin(θi,j), hence

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M49">View MathML</a>

(20)

with

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M50">View MathML</a>

By applying the CRB inequality (Equation 17) to the diagonal elements of F−1, we find that the estimation error on the x-dimension is lower-bounded by

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M51">View MathML</a>

(21)

and similarly, the error on the y-dimension by

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M52">View MathML</a>

(22)

If all sensors are located at equal distance Ri from the target and have the same DoA estimation variance <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M53">View MathML</a>, the above expressions are simplified to the following:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M54">View MathML</a>

(23)

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M55">View MathML</a>

(24)

These expressions show that the estimation error variance scales linearly with the DoA estimation variance and quadratically with the target distance.

Finally, we express the quantity <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M56">View MathML</a> that bounds the MSE norm of any unbiased estimator <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M57">View MathML</a>. From Equations 21 to 22 and defining <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M58">View MathML</a>, we obtain the following:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M59">View MathML</a>

(25)

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M60">View MathML</a>

(26)

The above expression is useful to identify non-observability conditions, i.e., configurations where <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M61">View MathML</a>. Consider for example the case M=2. Equation 26 becomes

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M62">View MathML</a>

(27)

Hence, the CRB is infinite if

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M63">View MathML</a>

(28)

for k = 0,1,2. This result is consistent with [14], where it was shown from geometric considerations that non-observability occurs if the sensors are placed on a straight line passing through the target. For M=2, this condition corresponds exactly to θi,1 − θi,2 = π.

Note that there is an infinite number of configurations satisfying Equation 28. Thus, if we assume angles θi,1 and θi,2 to be random and uniformly distributed in [0,2 π), a configuration with M = 2 sensors does not provide reliable PU localization. Mathematically, the reason is that if both θi,1 and θi,2 are uniformly distributed, so is θ = θi,1 − θi,2, and the expected value of Equation 27 with respect to θ tends to because

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M64">View MathML</a>

A physical interpretation is that the ambiguity of π in DoA estimation makes it impossible to distinguish between angles whose difference is π.

If we increase the number of sensors (M ≥ 3), configurations such that the denominator of Equation 26 is 0 (i.e., all sensors lying on the same straight line) become events with zero probability. We conclude that the minimum number of sensors should be M = 3.

4.2 Practical algorithms

In this section, we introduce different DoA-only localization methods (ML, LS, and Stansfield estimators) and discuss their application in CR networks. Before continuing, the following preliminary observations need to be made:

• For simplicity of presentation, we consider M DoA measurements to be collected by the fusion center at the same time. In practice, such measurements may also be obtained during different time slots. In this case, the considered estimators are implemented through iterative numerical methods (which, in fact, can be used in the case of simultaneous measurements as well). However, multiple measurements per sensor are beneficial only if sensors move; otherwise, they may result in matrix singularities which ultimately degrade the estimation quality.

• We do not consider the extended Kalman filter because this method provides poor performance when applied to DoA-only localization, as proven in [12].

• While CRB results presented in section 4.1 hold for an arbitrary number of coexisting signals (as long as n<m), in practice, data association of multiple DoAs observed by multiple sensors is a non-trivial problem (precisely, it is NP-hard. See, e.g., [22] and references therein). For this reason, in this paper, we assume ideal identification of PUs by SUs, i.e., either there is a single PU observed by multiple sensors or there are multiple PUs, but each of them can be unambiguously identified by SUs by exploiting other information, e.g., detection of carrier frequency, bandwidth, modulation type. We thus remove index i in the following (e.g., θj≡ θi,j).

4.2.1 Maximum-likelihood estimator

The ML estimator (see, e.g., [13,15,16]) is given as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M65">View MathML</a>

(29)

with

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M66">View MathML</a>

(30)

with <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M67">View MathML</a>, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M68">View MathML</a>, and <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M69">View MathML</a>. The symbol ⊖ denotes difference modulo 2π; we then definec:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M70">View MathML</a>

(31)

The minimization of Equation 29 is a non-linear least-squares problem and can be performed by iterative methods. The most popular of these methods is the Gauss-Newton algorithm [15,16], but it is very sensitive to the initial value. More robust iterative methods are the Levenberg-Marquardt algorithm [23] or the trust-region-reflective (TRR) algorithm [24].

As shown in [15], the ML algorithm, asymptotically in the number of measurements M, is unbiased and achieves the CRB. However, we will show in section 5.2 that the theoretical performance of the ML estimator is degraded when adopting iterative implementations as cases of non-convergence often occur.

4.2.2 Stansfield estimator

The Stansfield estimator, introduced in [20], is based on an approximation of the ML problem Equation 30 by letting Δθj≈ sin(Δθj), which holds for Δθj→0. This formulation results in a linearization of the originally non-linear ML problem so that it can be solved in closed form as follows (see [15] for more details):

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M71">View MathML</a>

(32)

where

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M72">View MathML</a>

()

and <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M73">View MathML</a>. Note that the Stansfield estimator involves the distances Rj between sensors and target; however, as discussed in [16] and confirmed by simulation results (see section 5.2), the dependency of the estimator on R is quite weak. For this reason, if distances are unknown (as it is the case in CR networks), R can be roughly estimated from the received signal strength or simply replaced by an identity matrix. In the rest of the paper, we will consider the latter case (worst-case scenario).

The Stansfield estimator was shown [15] to provide excellent performance in terms of mean square error (which can be even smaller than that of the ML estimator when M is small), but it was also shown to be biased (hence, not directly comparable to the CRB).

4.2.3 Least-squares estimator

Finally, a simple linear LS position estimate can be obtained by solving the system of equations:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M74">View MathML</a>

(33)

following directly from Equation 2 where DoA θj is replaced by its estimate. Writing in matrix form Equation 33 for j from 1 to M, we obtain the following:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M75">View MathML</a>

(34)

with

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M76">View MathML</a>

The solution is then given as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M77">View MathML</a>

(35)

Similar to ML and Stansfield estimators, a weighted version of the LS estimator can also be obtained:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M78">View MathML</a>

(36)

The linear LS estimator is not asymptotically unbiased: as shown in simulation results, its estimation accuracy tends to saturate as M grows.

5 Numerical results and analysis

In this section, we investigate the performance of the methods presented in the previous sections as a function of several systems parameters. We first consider DoA estimation (section 3): we compare the performance of MUSIC against the CRB and show how DoA estimation accuracy scales with the number of antennas m. We then consider the DoA-based localization algorithms presented in section 4: we compare the position error achieved by LS, ML, and Stansfield estimators and analyze the impact of the number of sensors M as well as the number of antennas per sensor m. We finally investigate the trade-off between the above parameters, in case of ideal conditions (uniform distance of all sensors from target, optimal orientation of sensor arrays) or realistic conditions (random distance and orientation).

5.1 DoA estimation

The DoA estimation variance achieved by the MUSIC algorithm using a ULA is given by Equation 14, and the corresponding CRB is given by Equation 12. In order to compare MUSIC performance against the CRB, we evaluate the DoA RMSE, i.e., the square root of the average (Δθ)2 over several simulations, with Δθ defined as in Equation 30. Without loss of generality, we assume ideal orientation of the array, such that <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M79">View MathML</a> in Equations 12 and 14. Different array orientations would result in a constant (in logarithmic scale) penalty term.

In Figure 2, the DoA RMSE, averaged over 104 iterations, is plotted as a function of the number of antennas m for different values of SNR (ρ = 0 dB and ρ = −5 dB). The number of received signal samples is set to N = 100. The following observations can be made:

Validation of MUSIC performance expressions. The formula of MUSIC estimation variance Equation 14 matches with the simulation results.

MUSIC vs. CRB. The gap between MUSIC estimation error and CRB is small and tends to reduce as the SNR increases (see section 3.2).

Impact of number of antennas. The DoA estimation accuracy (for both MUSIC and CRB) monotonically improves with m, following an almost linear trend on a logarithmic scale. No saturation is observed for practical values of m (i.e., m ≤ 10), which suggests that the more antennas are used by each sensor, the more accurate DoA estimates (and, therefore, PU localization) would be obtained.

thumbnailFigure 2. RMSE of the DoA estimation error vs. number of antennas.N = 100 received signal samples, 104 Monte Carlo iterations.

5.2 Localization

We now analyze the performance of the ML, LS, and Stansfield localization methods introduced in section 4.

5.2.1 Error distributions

In Figure 3, we evaluate the distribution of the position error defined as follows:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M80">View MathML</a>

(37)

thumbnailFigure 3. CDFs of position error<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M81">View MathML</a>(x is in distance units) for ML, LS, and Stansfield estimators. Random DoA estimation variance <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M82">View MathML</a> at different sensors.

computed over 1,000 simulations. The true PU position is set, without loss of generality, to xP = [0,0]T.

At every simulation, M sensor positions are selected randomly on a unit-radius circle centered in xP so that all sensors are at equal distance (Rj = 1 unit) from the target. The DoA estimation variance of each sensor (<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M83">View MathML</a>) is random, uniformly distributed between 10° and 30°. For each of the considered estimators, we test two versions: weighted (W) estimators (W-ML, W-Stansfield, W-LS) that assume knowledge of the values <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M84">View MathML</a>, i.e., of the matrix W in Equations 30, 32, and 36, respectively, and non-weighted estimators (ML, Stansfield, LS) that replace W by an identity matrix, thus assigning equal weight to all sensors.

The ML method is implemented iteratively, starting from an initial guess <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M85">View MathML</a> selected randomly at each Monte Carlo simulation, with σ0=1.5 distance units. The TRR implementation is adopted as it proves to be more robust than the Gauss-Newton algorithm (which does not converge if the initial guess is far from the true value).

Results, shown for three cases M = 3, M = 5, and M = 8, reveal the following:

• The Stansfield method provides the best performance among the three considered estimators. It reaches convergence in 100% of cases, while the ML CDF saturates around 90/95% due to cases of non-convergence (which become more frequent as M increases). Notice that suboptimality of the ML estimator is due exclusively to iterative implementation, which is sensitive to the initial estimate and may converge to local minima. The Gauss-Newton implementation provides even worse results, i.e., it fails to converge in 20%/25% of cases in the considered setting. The LS solution is close to the Stansfield solution for M = 3, but it becomes much more suboptimal as M grows.

• Weighted estimators provide in general a modest performance improvement over non-weighted ones, for the considered DoA variance range (<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M86">View MathML</a>). The gap increases for large M when more different DoA measurements (hence more weights) are available.

5.2.2 Position RMSE and CRB vs. number of sensors

We now compare the performance of ML, LS, and Stansfield estimators against the CRB derived in section 4.1, and we study their evolution with the sensor number M. To this purpose, we consider the position RMSE computed as the square root of the average square error <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M87">View MathML</a> over multiple iterations and random sensor positions. The lower bound of the position RMSE for unbiased estimators is given by the square root of the average of <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M88">View MathML</a> (see section 4.1) over the same random sensor positions.

In our simulation, we assume sensors to be equidistant from the target PU (Rj = 1 distance unit ∀ j) and, at every new random configuration, we generate M sensor positions as <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M89">View MathML</a> with ϕjU [0,2π). Angles ϕj are defined as the angles of sensors j with respect to target, i.e., ϕj = (θj− π) mod 2 π; hence, angles θj are also random with the same uniform distribution as ϕj. We generate 100 random sensor positions (i.e., 100 random vectors [ϕ1,…,ϕM]), where we compute the CRB and test the performance of LS, ML, and Stansfield methods by running 1,000 Monte Carlo simulations per configuration.

The simulation is repeated for values of M from 2 to 10. Note that LS and Stansfield estimators are not directly comparable to the CRB because they are not unbiased (so their RMSE could be, in principle, lower than the CRB); however, the CRB still represents an important benchmark to evaluate the quality of the position RMSE achieved by practical algorithms. For M = 2, when generating [ϕ1,ϕ2], we impose the constraints 5° <|ϕ1ϕ2|<175° so as to avoid configurations close to non-observability conditions (see section 4.1) that would result in RMSE and CRB → . In this case, the DoA estimation variance is constant (<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M90">View MathML</a>) so that weighted and non-weighted algorithms are identical.

Results are shown in Figure 4. The Stansfield method provides the best performance for all values of M, and it achieves the CRB for M > 3. The LS performance is almost equivalent to that of the Stansfield method for low numbers of sensors (M = 2 and M = 3), but it does not further improve as M grows; therefore, if the linear LS method is adopted, there is no advantage in having more than three or four sensors. The ML-TRR algorithm, finally, turns out to be suboptimal due to non-convergence issues. For a more fair comparison, we have plotted the performance curve of ML, ignoring cases of evident non-convergence (i.e., discarding values of RMSE >1,000). Still, this method is strongly suboptimal compared to both Stansfield and LS.

thumbnailFigure 4. Position RMSE vs. number of sensors, averaged over 100 random sensor positions, 1,000 iterations per sensor configuration.Rj = 1, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M91">View MathML</a>j.

5.2.3 Position RMSE and CRB vs. DoA (for M = 2)

To give more insight into the case M = 2, in Figure 5, we compare the CRB and the RMSE of the Stansfield method as a function of the DoA difference:

<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M92">View MathML</a>

(38)

thumbnailFigure 5. Position CRB and RMSE of the Stansfield estimator for M  = 2 as a function ofθΔ = |θ2θ1|.

Values of θΔ from 0 to 2π are considered, with a step of π/20. It can be observed that, as predicted by the theory, θΔ = 0 and θΔ = π are non-observability conditions; therefore, both the CRB and the Stansfield error tend to diverge as θΔ approaches such critical values. The divergence of the CRB is sharp, while the Stansfield method exhibits a smoother behavior. The Stansfield RMSE nearly attains the CRB for values of θΔ far from 0 and π, and it becomes suboptimal in the regions close to non-observability.

5.2.4 Position RMSE vs. number of antennas

Similar to section 5.1 where we have evaluated the impact of m on DoA estimation accuracy, we now investigate the resulting impact of m on the position RMSE. Results shown in Figure 6 refer to a scenario with M = 5 sensors and assume the same number of antennas m, ranging from 2 to 10, for all sensors. The SNR is set to ρ=0 dB; the number of received signal samples used by MUSIC is N = 100. The CRB reported in the plot is computed with <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M93">View MathML</a> given by Equation 12.

thumbnailFigure 6. Average position RMSE (100 sensor configurations×1,000 iterations) vs. number of antennas per sensor.M = 5, N = 100, ρ = 0 dB.

We observe that, both from the point of view of the CRB and of practical algorithms, the position RMSE is reduced by increasing the number of antennas. Thus, the localization accuracy as a function of the number of antennas is not affected by saturation phenomena that occur when increasing the number of sensors (see Figure 4). We also note that the RMSE of the Stansfield algorithm tends to achieve the CRB as m.

5.3 Trade-off between number of sensors and number of antennas

So far, we have analyzed the position accuracy as a function of different system parameters, in particular, the number of sensors M and the number of antennas m. A natural question may be asked: which one of these two factors is more important with regard to PU localization accuracy? To address this question, let us consider a CR network with a total of K antennas; assume that such antennas can be arbitrarily organized into arrays (for simplicity, ULAs) to localize a PU transmitter; then, if every array has m elements, the number of arrays (hence, of DoA measurements) is M = ⌊K/m⌋.

We then formulate the problem as follows: given a fixed number of antennas K, we want to determine the optimal trade-off [m,M] between number of antennas per sensor and number of sensors, with m×MK. It has been shown in Figure 6 that increasing m constantly improves the localization accuracy; on the contrary, increasing the number of sensors besides a certain value (around M = 4) does not provide significant advantages, as shown in Figure 4. For this reason, we expect that the best trade-off will be obtained for m>M.

We next consider three cases: (1) uniform sensor-target distance and ideal array orientation for all sensors (Rj=1 unit, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M94">View MathML</a>j); (2) random distance: RjU[0.5,1.5] units ∀j; and (3) random distance and random array orientation: <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M95">View MathML</a>. Results are shown, respectively, in Figures 7, 8, 9. The total number of antennas is set as K = 36, a number that allows for several exact combinations [m,M] such that m × M = 36. The pair [m = 18,M = 2] is excluded a priori, because M = 2 gives rise to non-observability conditions, as discussed in section 4.1. The performance of LS and Stansfield estimators is plotted along with two CRBs: one computed with <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M96">View MathML</a> given by Equation 12 and the other one with <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M97">View MathML</a> given by Equation 14. Again, the number of samples used in MUSIC is N = 100 and the SNR is ρ = 0 dB. All simulations are averaged over 100 random sensor configurations and 1,000 iterations per configuration.

thumbnailFigure 7. Antennas-sensors trade-off: uniform sensor-target distance, ideal array orientation.

thumbnailFigure 8. Antennas-sensors trade-off: random sensor-target distance, ideal array orientation.

thumbnailFigure 9. Antennas-sensors trade-off: random sensor-target distance, random array orientation.

In ideal conditions (Figure 7), the optimal trade-off is located, as expected, in the region m > M for all curves. In particular, according to the CRB and to the Stansfield method, the optimal combination is [m = 9,M = 4], while using the LS method, the optimum is at [m = 12, M = 3]. Such discrepancy is explained by the results of Figure 4: for a given DoA estimation variance, the LS localization RMSE does not significantly improve from M = 3 to M = 4; therefore, it is more advantageous to keep M = 3 and increase the number of antennas from 9 to 12 so as to improve <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M98">View MathML</a>. We note, then, that the gap between the two CRBs is very moderate (i.e., MUSIC by itself provides excellent DoA estimation performance, as already observed in Figure 2), and that the Stansfield RMSE equals CRB(<a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M99">View MathML</a>) except for M = 3 (cf. Figure 4).

When introducing variable distances (Figure 8), a similar behavior is observed. The only difference is that the Stansfield estimator becomes slightly suboptimal compared to the CRB because we assume distances to be unknown, so R is still considered an identity matrix in Equation 32. The conclusion is that different distances of sensors from the target only have a minor impact on the resulting localization accuracy,

Very different results are obtained when considering non-ideal array orientation (Figure 9). Note that the assumption of random uniformly distributed <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M100">View MathML</a> may correspond either to a scenario where each sensor has a different, random orientation or to a scenario where all sensors have equal orientation (in this case, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M101">View MathML</a> + a constant offset; therefore, <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M102">View MathML</a> is random like θj). We consider here two variants of the LS and Stansfield estimators: one where the orientation is known and coefficients <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M103">View MathML</a> are incorporated in the weight matrix W, and one where the orientation is unknown, so W = I. Figure 9 shows that ‘weighted’ estimators, W-LS and W-Stansfield, achieve good performance (close to the CRB for W-Stansfield), while standard estimators are strongly suboptimal and provide nearly constant RMSE for all values of m and M. The best trade-off is located at values of M lower than in the previous scenarios: [m = 6, M = 6] for Stansfield and CRB, and [m = 9,M = 4] for LS, i.e., one step less than in the other cases (Figures 7 and 8). This behavior can be explained by observing that, when each sensor has a random orientation, more ‘sensor diversity’ is necessary to achieve good performance. In terms of absolute values, the best RMSE achievable with random orientation is approximately 0.02, while it is approximately 0.008 in the case of ideal orientation. The gap is approximately 4 dB.

The above results show that it is important that the array orientation is taken into account by SUs. We remark that <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M104">View MathML</a> is simply the DoA of the impinging wave, which is estimated through MUSIC; therefore, it is sufficient to use the MUSIC DoA estimates <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M105">View MathML</a> to fill the weight matrix W. Then, if sensors are able to rotate, another possibility is to rotate the array towards <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M106">View MathML</a> so as to optimize the orientation factor. This can be implemented by an iterative procedure. First the sensor, using MUSIC, estimates the DoA <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M107">View MathML</a>, and then, it rotates by <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M108">View MathML</a>, so that the new orientation is <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M109">View MathML</a>. At this point, if a new estimate of <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M110">View MathML</a> is taken using MUSIC, the DoA estimation variance gets reduced by a factor <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M111">View MathML</a> as followed from Equation 14. Then, if needed, the array rotation can be readjusted based on the new estimate. Iteration can be stopped when the estimated DoA stabilizes around <a onClick="popup('http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://jwcn.eurasipjournals.com/content/2013/1/107/mathml/M112">View MathML</a>, i.e., when the array is approximately orthogonal to the signal direction. The results of the above array orientation procedure are illustrated in Figure 10. The plots show the convergence of the position RMSE as a function of the number of iterations, averaged over 104 Monte Carlo simulations for a given sensor configuration, with ρ = 0 dB, N = 100 samples, m=6 antennas per sensor, and M = 6 sensors. The iterative algorithm converges very quickly, with a steep error reduction after one iteration. The resulting position RMSE nearly achieves the CRB, like in the case of ideal array orientations (Figure 8). Similar results are obtained with different values of m and M as well.

thumbnailFigure 10. Convergence of position RMSE vs. iterations of the array orientation procedure.N=100 received signal samples, SNR ρ=0 dB, 104 Monte Carlo iterations.

6 Implementation issues

The key requirements for practical implementation of the proposed localization method in a cognitive radio network can be listed as follows:

1. Availability or multi-antenna SU devices (possibly with rotation capabilities),

2. Existence of a wired or wireless link for sending DoA measurements from distributed sensors to a fusion center, and

3. Sufficient processing power and energy supply to support DoA estimation at the sensor nodes and DoA fusion at the fusion center.

Regarding the first point, multi-antenna cognitive radios are nowadays a common assumption in the spectrum sensing literature (see, for example, [5] and references therein) and can be therefore exploited to combine DoA estimation with traditional signal detection with little additional complexity. Furthermore, the number of antennas we considered in our numerical examples is within practical limits of current base stations and Wi-Fi access points. We also remark that, if antenna arrays are not available, directional antennas or electronically steerable antennas [25] can be used as well. A more challenging requirement is the availability of rotating arrays, which may be feasible only in certain scenarios (e.g., vehicle terminals, dedicated spectrum monitors, base stations). However, the availability of such type of devices is not a prerequisite of the algorithm but only a possible way to optimize performance (see Figure 10). In the absence of rotating arrays, localization performance can be still improved by increasing the number of sensors or using proper weighting matrices (see Figure 9).

Coming to the second point, let us first consider the case of physical antenna arrays (SUs equipped with multiple antennas). The proposed algorithm involves exchange of DoA measurements from M−1 distributed sensors to a fusion center (it is assumed that the fusion center itself is one of the M cooperating sensors). In addition, if the array orientation of sensor nodes is unknown to the fusion center, M−1 additional angle measurements should be sent from the distributed sensors in order to establish a common coordinate system. Hence, if a scalar angle measurement is encoded by one ‘information unit’ (e.g., 1 byte), the signaling overhead of the algorithm amounts to 2(M − 1) = O(M) information units, i.e., it scales linearly with the number of sensors. Given the moderate bandwidth requirement, the algorithm lends itself well to both wired and wireless implementations. In the latter case, the sensors can communicate to the fusion center via a dedicated channel (‘common control channel’ [26,27]), which may be either predefined or selected dynamically after a preliminary spectrum sensing phase. In the case of virtual arrays (m single-antenna users performing MUSIC cooperatively), there is an additional signaling overhead within each sensor array. Assuming that one user acts as a local fusion center, each of the remaining M − 1 users must transmit N received signal samples over the wireless channel, hence the signaling overhead per sensor is N(m − 1) complex numbers. In this case, the number of samples (N) should not be too large, and synchronization between the users within a virtual array is crucial.

Finally, numerical complexity (and, consequently, energy consumption) is dominated by the eigenvalue or singular value decomposition (EVD/SVD) necessary to implement MUSIC at each sensor node. The number of floating-point operations required for EVD (applied to the m × m covariance matrix) is ≈ 10m3, while for SVD (applied to the m × N data matrix) it is ≈k1Nm2+k2m3, where k1,k2 are constants [28,29]. Therefore, in general, the complexity of EVD/SVD scales as O(m3). Clearly, this task may be critical in terms of processing power and battery consumption for sensor nodes. We note, however, that (1) according to our simulation results, accurate localization can be achieved with a limited number of antennas per sensor, resulting in an affordable complexity for typical state-of-the-art embedded processors (for a practical DSP implementation of MUSIC, we refer to [30]); (2) although we have considered here the classic MUSIC algorithm, low-complexity variants exist in the literature (see, for example, [31-33]). The DoA fusion step (using LS or Stansfield estimation) then involves solving a system of linear equations (Equations 32 and 34). This task can be performed by LU factorization and Gaussian elimination [29] with complexity O(m3). Since it is reasonable to assume for the fusion center higher computation power than the other sensor nodes and wired connection to the network backbone, the computational burden of DoA fusion is generally less critical than the one involving distributed sensors.

Signaling overhead and numerical complexity, with respect to the distributed sensors and the fusion center, are summarized in Table 1. Note that in the first case (multi-antenna sensors), there is no signaling towards the distributed sensors but only from the sensors to the fusion center. Also note that the fusion center performs first (local) DoA estimation and then fusion of DoA estimates from all sensors; since both tasks have cubic complexity (respectively, in the number of antennas and in the number of sensors), the resulting complexity is O(max {m,M}3).

Table 1. Signaling overhead and numerical complexity

We finally remark that the proposed method is robust to possible failures of sensor nodes (due to low battery or other reasons), as long as at least three sensors remain operational.

7 Conclusions

In this paper, we have investigated PU localization in CR networks by DoA-based techniques, where DoA estimates are obtained by multi-antenna sensors. The key advantage of the proposed approach is that it does not rely on range measurements (difficult to obtain in a CR network where PUs do not collaborate with SUs) nor on precise synchronization of different sensors (needed for TDoA techniques).

The following conclusions can be drawn from the results presented in the paper:

• Very accurate PU localization can be obtained using the MUSIC algorithm for DoA estimation at each sensor and the Stansfield method to combine DoA measurements together. In this way, the PU position RMSE nearly achieves the CRB.

• Other DoA fusion methods, namely the linear LS estimator and the ML estimator implemented iteratively, are inferior to the Stansfield estimator.

• In order to maximize the localization accuracy, the effect of the number of antennas per sensor prevails over the effect of the number of sensors.

• If sensor arrays have random orientation, the performance is affected by a major penalty factor, which should be compensated by properly weighting DoA measurements from different sensors or rotating sensors (when possible) so as to optimize their orientation.

Endnotes

aIn the navigation context, it typically assumed that multiple observations are obtained by a unique sensor at different positions (i.e., a maneuvering ship), whereas in CR networks we assume that observations are taken by multiple sensors at the same time, and then gathered by one node acting as a fusion center.

bIf subspace methods (e.g., MUSIC) are used, the maximum number of identifiable primary signals is M − 1. Thus, we assume n<m.

cThe adopted formulation ensures that the difference between two arbitrary angles is always computed correctly (e.g., |1° ⊖351° | = 10°. On the contrary, |1° − 351 ° | = 1 ° ⊖351° =  350° is formally correct too, but leads to overestimate an error that is in fact small, which is a problem for iterative implementation of the ML method).

dTRR is implemented using the Optimization Toolbox of Matlab v7.8.

eAn exact ML solution would give at least the same performance of the Stansfield estimator, which is in fact an approximation of the ML estimator for Δ θ →0. However, an ML implementation through exhaustive search (e.g., on a grid of points) would be prohibitive in terms of complexity and still suboptimal due to discretization of the search space.

fThe main source of complexity of DoA estimation is given by the eigenvalue decomposition of the covariance matrix, which is also required for most of the existing multi-antenna detection techniques.

gThe sensors are assumed to be aware of their own position as well as that of the fusion center so that they are able to estimate the relative angle. This step is not necessary if the positions of sensor nodes are fixed and/or known in advance at the fusion center.

Competing interests

Both authors declare that they have no competing interests.

Authors’ information

FP received his M.S. degree (with Hons.) and Ph.D. degree in communications engineering from the Politecnico di Torino, Italy in 2008 and 2012, respectively. He is currently a research associate with the Fraunhofer Heinrich Hertz Institute in Berlin, Germany. His research interests include cognitive radio networks, distributed signal processing, multi-sensor signal detection, and cooperative localization. DC received her Dipl. Ing. degree from the University of Belgrade, Serbia in 1998 and her M.Sc. degree in electrical engineering from the University of California, Los Angeles (UCLA) in 2001. She received her Ph.D. degree in electrical engineering from the University of California, Berkeley in 2007. She was a member of the Berkeley Wireless Research Center, University of California, Berkeley. In 2008, she joined the faculty of Electrical Engineering, UCLA, as an assistant professor. Her key contributions involve the novel radio architecture, signal processing, and networking techniques to implement spectrum sensing functionality in cognitive radios. DC was awarded a Samueli Fellowship in 2008 and an Okawa Foundation research grant in 2009.

Acknowledgements

This research was supported by the National Science Foundation under CNS grant 1149981 and by the German Federal Ministry of Education and Research under grant 01BU1224. Some preliminary results of this paper were presented at the IEEE Global Communications Conference (GLOBECOM), Houston, TX, USA, December 2011 [34].

References

  1. J Mitola, GQ Maguire, Cognitive radios: making software radios more personal. IEEE Personal Commun 6(4), 13–18 (1999). Publisher Full Text OpenURL

  2. S Haykin, Cognitive radio: brain-empowered wireless communications. IEEE Trans. Commun 23(2), 201–220 (2005)

  3. D Cabric, SM Mishra, RW Brodersen, Implementation issues in spectrum sensing for cognitive radios. Proc. IEEE Asilomar Conf. Signals, Syst. Comput 1, 772–776 (2004)

  4. A Ghasemi, ES Sousa, Spectrum sensing in cognitive radio networks: requirements, challenges and design trade-offs. IEEE Commun. Mag 46(4), 32–39 (2008)

  5. Y Zeng, YC Liang, AT Hoang, R Zhang, A review on spectrum sensing for cognitive radio: challenges and solutions. EURASIP J. Adv. Signal Process 2010, 1–15 (2010)

  6. FCC, Second memorandum of order:. September 23, 2013 (2013) (http://hraunfoss, 2013), . fcc.gov/edocs public/attachmatch/FCC-10-174A1.pdf webcite. Accessed 17 April 2013 OpenURL

  7. R Tandra, SM Mishra, A Sahai, What is a spectrum hole and what does it take to recognize one? Proc of the IEEE 97(5), 824–848 (2009)

  8. S Kandeepan, S Reisenfeld, TC Aysal, D Lowe, R Piesiewicz, Bayesian tracking in cooperative localization for cognitive radio networks. Proceedings of IEEE Vehicular Technology Confernece (VTC 2009 Spring) Barcelona April 2009 (New York: IEEE, 2009), pp. 1–5

  9. Z Ma, W Chen, KB Letaief, Z Cao, A semi range-based iterative localization algorithm for cognitive radio networks. IEEE Trans. Vehicular Technol 59(2), 704–717 (2010)

  10. R Martin, Algorithms and bounds for estimating location, directionality, and environmental parameters of primary spectrum users. IEEE Trans. Wireless Commun 8(11), 5692–5701 (2009)

  11. J Wang, P Urriza, Y Han, D Cabric, Performance analysis of weighted centroid algorithm for primary user localization in cognitive radio networks. Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, November 2010 (New York: IEEE, 2010), pp. 966–970

  12. VJ Aidala, Kalman filter behavior in bearings-only tracking applications. IEEE Trans. on Aerospace and Electron. Syst 15(1), 29–39 (1979)

  13. SC Nardone, AG Lindgren, KF Gong, Fundamental properties and performance of conventional bearings-only target motion analysis, IEEE Trans. Automatic Control 29(9), 775–787 (1984). Publisher Full Text OpenURL

  14. M Gavish, E Fogel, Effect of bias on bearing-only target location. IEEE Trans. on Aerospace Electron. Syst 26(1), 22–26 (1990). Publisher Full Text OpenURL

  15. M Gavish, AJ Weiss, Performance analysis of bearing-only target location algorithm. IEEE Trans. on Aerospace Electron. Syst 28(3), 817–828 (1992). Publisher Full Text OpenURL

  16. Y Oshman, P Davidson, Optimization of observer trajectories for bearings-only target localization. IEEE Trans. Aerospace Electron. Syst 35(3), 892–902 (1999). Publisher Full Text OpenURL

  17. N Patwari, JN Ash, S Kyperountas, AO Hero III, RL Moses, NS Correal, Locating the nodes: cooperative localization in wireless sensor networks. IEEE Signal Process. Mag 22(4), 54–69 (2005)

  18. H Krim, M Viberg, Two decades of array signal processing research: the parametric approach. IEEE Signal Process. Mag 13(4), 67–94 (1996). Publisher Full Text OpenURL

  19. R Schmidt, Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propagation 34(3), 276–280 (1986). Publisher Full Text OpenURL

  20. RG Stansfield, Statistical theory of DF fixing. J. IEE - Part IIIA: Radiocommunication 94(15), 762–770 (1947)

  21. P Stoica, A Nehorai, Music likelihood, maximum Acoustics, Cramer-Rao bound. IEEE Trans. Speech, Signal Process 37(5), 720–741 (1989). Publisher Full Text OpenURL

  22. H Gauvrit, JP Le Cadre, C Jauffret, A formulation of multitarget tracking as an incomplete data problem. IEEE Trans. Aerospace Electron. Syst 33(4), 1242–1257 (1997)

  23. D Marquardt, An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math 11(2), 431–441 (1963). Publisher Full Text OpenURL

  24. TF Coleman, Y Li, An interior, trust region approach for nonlinear minimization subject to bounds. SIAM J Optimization 6(2), 418–445 (1996). Publisher Full Text OpenURL

  25. T Ohira, K Gyoda, Electronically steerable passive array radiator antennas for low-cost analog adaptive beamforming. IEEE International Conference on Phased Array Systems and Technology, Dana Point May 2000 (New York: IEEE, 2000), pp. 101–104

  26. IF Akyildiz, WY Lee, MC Vuran, S Mohanty, Next generation / dynamic spectrum access / cognitive radio wireless networks: a survey. ELSEVIER Int. J. Comput. Telecommunications Netw. (ComNet) 50(13), 2127–2159 (2006). Publisher Full Text OpenURL

  27. N Nie, C Comaniciu, Adaptive channel allocation spectrum etiquette for cognitive radio networks. IEEE Symposium on Dynamic Spectrum Access Networks (DySPAN), Baltimore, November 2005 (New York: IEEE, 2005), pp. 269–278

  28. R Roy, T Kailath, ESPRIT - estimation of signal parameters via rotational invariance techniques. IEEE Trans. Signal Process 37(7), 984–995 (1989)

  29. GH Golub, Loan Van C F., Matrix Computations (Baltimore: John Hopkins University Press, 1996)

  30. M Kim, K Ichige, h Arai, Implementation of FPGA based fast DOA estimator using unitary MUSIC algorithm. IEEE Vehicular Technology Conference (VTC 2003-Fall), Orlando October 2003 (New York: IEEE, 2003), pp. 213–217

  31. J Lee, T Song, H Kwon, S Lee, Low-complexity estimation of 2D DOA for coherently distributed sources. ELSEVIER Signal Process 83(8), 1789–1802 (2003). Publisher Full Text OpenURL

  32. CS Marino, PM Chau, Reduced complexity covariance matrix estimate for subspace-based array processing. IEEE Asilomar Conference on Signals, Systems and Computers, Pacific Grove, November 2003 (New York: IEEE, 2003), pp. 790–794

  33. Y Fenggang, J Ming, X Qiao, Low-complexity DOA estimation based on compressed MUSIC and its performance analysis. IEEE Trans. Signal Process 61(8), 1915–1930 (2013)

  34. F Penna, D Cabric, Bounds and tradeoffs for cooperative DoA-only localization of primary users. IEEE Global Communications Conference (GLOBECOM), Houston, December 2011 (New York: IEEE, 2011), pp. 1–5