Design and performance analysis of a fully distributed source detection algorithm for WSNs

TL;DR

This paper presents a fully distributed source detection algorithm for wireless sensor networks that efficiently cooperates among nodes to detect signals without central coordination, with performance close to the optimal global test.

Contribution

It introduces a novel distributed GLR-based detection scheme that requires minimal communication and local network knowledge, suitable for resource-constrained WSNs.

Findings

Distributed GLR test matches global test performance.

Algorithm requires few communication exchanges.

Effective with limited measurements per sensor.

Abstract

In this article, we consider the detection of a localized source emitting a signal using a wireless sensor network (WSN). We consider that geographically distributed sensor nodes obtain energy measurements and compute cooperatively and in a distributed manner a statistic to decide if the source is present or absent without the need of a central node or fusion center (FC). We first start from the continuous-time signal sensed by the nodes and obtain an equivalent discrete-time hypothesis testing problem. Secondly, we propose a fully distributed scheme, based on the well-known generalized likelihood ratio (GLR) test, which is suitable for a WSN, where resources such as energy and communication bandwidth are typically scarce. In third place, we consider the asymptotic performance of the proposed GLR test. The derived results provide an excellent matching with the scenario in which only a finite amount of measurements are available at each sensor node. We finally show that the proposed distributed algorithm performs as well as the global GLR test in the considered scenarios, requiring only a small number of communication exchanges between nodes and a limited knowledge about the network structure and its connectivity.

Equations

$$\left\{\begin{array}[]{ll}\!\!\!\mathcal{H}_{0}:y_{k}(t)=v_{k}(t),&\!\!\!t\in[0,\tau],k\in[1:N]\\ \!\!\!\mathcal{H}_{1}:y_{k}(t)=h_{k}s(t)+v_{k}(t),&\!\!\!t\in[0,\tau],k\in[1:N]\end{array}\right.$$

$${\tilde{z}}_k(l)=\frac{1}{T}{\int }_{(l-1)T}^{lT}|y_k(t){|}^{2}dt,l\in [1:L],k\in [1:N].$$

$${\int }_0^T\frac{\sin (2\pi W(s-T/2))}{\pi (s-T/2)}{\varphi }_p(s)ds={\lambda }_p{\varphi }_p(s),$$

$$\left\{\begin{array}[]{llll}\mathcal{H}_{0}:&\tilde{z}_{k}(l)\approx\frac{1}{M}\sum^{M}_{p=1}|a_{l,p,k}|^{2}\\ \mathcal{H}_{1}:&\tilde{z}_{k}(l)\approx\frac{1}{M}\sum^{M}_{p=1}|h_{k}s_{l,p}+a_{l,p,k}|^{2}.\end{array}\right.$$

$$c_k\equiv {\sigma }_s^2|h_k{|}^2/{\sigma }_v^2>0.$$

$$\left\{\begin{array}[]{llll}\mathcal{H}_{0}:&\boldsymbol{\{}\tilde{\boldsymbol{z}}_{l}\}_{l=1}^{L}\overset{iid}{\sim}\mathcal{N}(\sigma_{v}^{2}\boldsymbol{1}_{N},\frac{\sigma_{v}^{4}}{M}\mathbf{I}_{N})\\ \mathcal{H}_{1}:&\{\tilde{\boldsymbol{z}}_{l}\}_{l=1}^{L}\overset{iid}{\sim}\mathcal{N}(\tilde{\boldsymbol{\mu}}(\mathbf{c}),\tilde{\mathbf{\Sigma}}(\mathbf{c})),\end{array}\right.$$

$$\left\{\begin{array}[]{llll}\mathcal{H}_{0}:&\boldsymbol{\{}\boldsymbol{z}_{l}\}_{l=1}^{L}\overset{iid}{\sim}\mathcal{N}(\boldsymbol{0},\mathbf{I}_{N})\equiv p(\boldsymbol{z}_{l};\boldsymbol{\theta}=\boldsymbol{0}_{N})\\ \mathcal{H}_{1}:&\{\boldsymbol{z}_{l}\}_{l=1}^{L}\overset{iid}{\sim}\mathcal{N}(\boldsymbol{\mu}(\boldsymbol{c}),\mathbf{\Sigma}(\boldsymbol{c}))\equiv p(\boldsymbol{z}_{l};\boldsymbol{\theta}=\boldsymbol{c}),\end{array}\right.$$

$$\begin{array}[]{ll}\boldsymbol{\mu}(\boldsymbol{c})&=\sqrt{M}\boldsymbol{c}\\ \mathbf{\Sigma}(\boldsymbol{c})&=\boldsymbol{c}\boldsymbol{c}^{T}+2\text{diag}(\boldsymbol{c})+\mathbf{I}_{N}.\end{array}$$

$$\begin{gathered} {\hat{c}}_{\mathrm{L}-\mathrm{MLE},k}=\frac{1}{2}(\sqrt{(M+2+\sqrt{M}{m}_k{)}^2+4(p_z+\sqrt{M}{m}_k-1)} \\ -(M+2+\sqrt{M}{m}_k))k\in [1:N], \end{gathered}$$

$$2\log T_{G}(\boldsymbol{z})\overset{a}{\sim}\left\{\begin{array}[]{ll}\chi^{2}_{N}&\text{under }\mathcal{H}_{0}\\ \chi^{\prime 2}_{N}(\lambda_{g})&\text{under }\mathcal{H}_{1},\end{array}\right.$$

$$T_L(\mathbf{z})\equiv \frac{p(\mathbf{z};\mathbf{\theta }={\hat{\mathbf{c}}}_{\mathrm{L}-\mathrm{MLE}})}{p(\mathbf{z};\mathbf{\theta }=0)}.$$

$$\int \cdots \int p({\mathbf{z}}_l;\mathbf{\theta })d{z}_1(l)\dots d{z}_{k-1}(l)d{z}_{k+1}(l)\dots d{z}_N(l)$$

$$=p_k(z_k(l);{\theta }_k).$$

$${\hat{\theta }}_{loc,k}^L=\arg \underset{{\theta }_k\in {\mathcal{A}}^k}{\max }\frac{1}{L}\sum _{l=1}^L\log p_k(z_k(l);{\theta }_k),$$

$$2\log T_L(\mathbf{z})$$

$$\log T_L(\mathbf{z})$$

$$=-\frac{L}{2}\log \det (\hat{\mathbf{\Sigma }})+\frac{1}{2}\sum _{l=1}^L\left\{\Vert {\mathbf{z}}_l{\Vert }^2-\Vert {\hat{\mathbf{\Sigma }}}^{-\frac{1}{2}}({\mathbf{z}}_l-\hat{\mathbf{\mu }}){\Vert }^2\right\}$$

$$\stackrel{(a)}{=}-\frac{L}{2}(\log (1+{\bar{c}}_1)+{\bar{c}}_2)$$

$$+\frac{1}{2}\sum _{l=1}^L\sum _{k=1}^N\left\{z_k^2(l)-\frac{(z_k(l)-\sqrt{M}{\hat{c}}_k{)}^2}{1+2{\hat{c}}_k}\right\}$$

$$+\frac{1}{2(1+{\bar{c}}_1)}\sum _{l=1}^L{\left\{\sum _{k=1}^N(z_k(l)-\sqrt{M}{\hat{c}}_k)\frac{{\hat{c}}_k}{1+2{\hat{c}}_k}\right\}}^2$$

$$u_k\equiv \sum _{l=1}^L\left\{z_k^2(l)-\frac{(z_k(l)-\sqrt{M}{\hat{c}}_k{)}^2}{1+2{\hat{c}}_k}\right\}.$$

$${\left\{\sum _{k=1}^N\frac{1}{L}\sum _{l=1}^L(z_{kl}-\sqrt{M}{\hat{c}}_k)\frac{{\hat{c}}_k}{1+2{\hat{c}}_k}\right\}}^2=$$

$${\left\{\sum _{k=1}^N(m_k-\sqrt{M}{\hat{c}}_k)\frac{{\hat{c}}_k}{1+2{\hat{c}}_k}\right\}}^2,$$

$$T_{L-FD}(\mathbf{z})$$

$$a_k(t)=W_{kk}{a}_k(t-1)+\sum _{j\in {\mathcal{N}}_k}{W}_{kj}{a}_j(t-1),k\in \mathcal{N},t\in \mathbb{N}$$

$$\mathbf{a}(t)=\mathbf{W}\mathbf{a}(t-1)={\mathbf{W}}^t\mathbf{a}(0),t=1,2,\dots$$

$$A_{kj}=(\mathbf{A})_{kj}=\left\{\begin{array}[]{rl}1,&\text{if edge $j$ starts from node $k$,}\\ -1,&\text{if edge $j$ ends at node $k$,}\\ 0,&\text{otherwise.}\end{array}\right.$$

$$\mathbf{W}={\mathbf{I}}_N-\mathbf{A}\text{diag}(\mathbf{w}){\mathbf{A}}^T.$$

$$\log p({\mathbf{z}}^k;c_k)\propto -L\log (1+c_k)-\frac{1}{2}{\Vert \frac{{\mathbf{z}}^k-\sqrt{M}{c}_k{1}_L}{1+c_k}\Vert }^2.$$

$$c_k^2+\left(M+2+\sqrt{M}{m}_k\right)c_k-\left(p_k+\sqrt{M}{m}_k-1\right)=0.$$

$$\mathbf{i}(0)={\mathbb{E}}_{\mathbf{\theta }=\mathbf{c}}{\left[\frac{\partial \log p({\mathbf{z}}_l;\mathbf{\theta })}{\partial \mathbf{\theta }}\frac{\partial \log p({\mathbf{z}}_l;\mathbf{\theta })}{\partial {\mathbf{\theta }}^T}\right]|}_{\mathbf{\theta }=0}$$

$$[\mathbf{i}(\mathbf{\theta }){]}_{kj}$$

$$+\frac{1}{2}\text{tr}\left[\mathbf{\Sigma }(\mathbf{\theta }{)}^{-1}\frac{\partial \mathbf{\Sigma }(\mathbf{\theta })}{\partial {\theta }_k}\mathbf{\Sigma }(\mathbf{\theta }{)}^{-1}\frac{\partial \mathbf{\Sigma }(\mathbf{\theta })}{\partial {\theta }_j}\right].$$

$${\mathbf{a}}^T(\mathbf{f}(\mathbf{y})-\mathbf{f}(\mathbf{x}))={\mathbf{a}}^{T}J(\mathbf{w})(\mathbf{y}-\mathbf{x}),$$

Full text

I Introduction

In the near past, Wireless Sensor Networks (WSN) have received considerable attention from the research and industrial community because of their remote monitoring and control capabilities [1, 2, 3]. More recently, they have become an essential part of the emerging technology of Internet of Things (IoT) [4, 5]. Among the different tasks to be done by WSNs, distributed detection is an actively researched topic [6, 7, 8].

Distributed detection architectures can be broadly classified in two classes. In the first class all sensors transmit their local measurements to a fusion center (FC), where some processing tasks are done and the final decision about the underlying phenomenon is made [9, 10, 11, 12]. In many applications it is unfeasible or expensive to develop an infrastructure with a FC. This centralized architecture also presents some weaknesses as, for example, its lack of robustness against the malfunctioning of a single device, given that a failure in the FC may severely degrade the performance of the system. Additionally, it requires that sensor nodes, typically battery-powered devices, communicate through orthogonal channels with the FC, consuming excessive energy and bandwidth. A way to circumvent this issue is to quantize the measurements to few bits (binary quantization is a popular choice) to save bandwidth. However, this strategy involves the design of the quantizers, which can be a hard task when the observations are correlated typically resulting in complicated decision rules [13, 14] even under the assumption of Gaussian data and networks with only a few nodes [15].

An alternative to the above described architecture is to consider distributed strategies for which there is not a central processing unit or FC. In this type of detection architectures, sensors distributed geographically, collect measurements from the phenomenon of interest, make some processing, exchange information with their neighbors and, finally, execute some consensus or diffusion algorithm to achieve their respective decisions. This option is robust against node failures, and the communications between nodes are done locally, over typically short distances, saving energy and also bandwidth, by employing spatial reuse of the frequency bands. Thus, the quantization of the measurements can be done with more levels and it becomes a less relevant problem.

Many works have considered the second option also known as a fully distributed detection architecture [16, 8, 17, 18]. Nevertheless, most of the work found in the literature assumes that the spatial measurements are independent or, they ignore the statistical dependence of the data when designing the distributed detection algorithms [19, 20, 21, 22, 23, 24]. For example, Cattivelli et al. proposed in [25] a distributed detection algorithm to detect a known deterministic signal under Gaussian noise, where the noise is assumed to be independent across the sensors, and thus the observations at each node are independent under each hypothesis. However, in many applications of interest, the measurements taken by spatially distributed nodes are statistically dependent, and disregarding this effect markedly degrades the detection performance of the network [26].

Other works have considered dependent observations using Gaussian Markov Random Fields [27, 28] to design a Neyman-Pearson detector in a centralized scenario. However, the design of distributed detection algorithms with dependent measurements in a decentralized scenario deserves more investigation [29].

In this work we deal with spatially correlated observations and propose a fully distributed algorithm to detect the presence or absence of a localized source emitting a stochastic signal. This problem is important in its own with multiple applications in fields as cognitive radio [30], massive MIMO wireless networks [31] and acoustic source detection, separation and localization [32], among others.

I-A Main contributions

The main contributions of the work can be summarized as follows. First, we develop a model for the problem of source detection. We assume that under the null hypothesis ($\mathcal{H}_{0}$) the signal is absent and under the alternative hypothesis ($\mathcal{H}_{1}$) it is present. The location of the source is unknown along with other parameters of the stochastic signal $s(t)$ that models the signal emitted by the source. Also, our network model does not assume the presence of a FC. The desired goal is the distributed detection of the source signal if present (that is, all the sensor nodes have to reach the same decision). Assuming that the nodes sense the energy of a signal, we are able to model the statistical dependence between samples in different nodes under both hypotheses. To make the problem tractable, we use the Central Limit Theorem to approximate the statistics of the observations by a multivariate Gaussian distribution under each hypothesis, where the covariance matrix under $\mathcal{H}_{1}$ has a particular structure that can be exploited to simplify the detection algorithm.

Secondly, we compute a modified version of the generalized likelihood ratio (GLR) detector, which estimates the unknown parameters locally at each node, instead of doing that globally, which would consume more network resources and would require a distributed solution of a complex optimization problem. We also derive its asymptotic performance and prove that under mild conditions, it coincides with the asymptotic performance of the statistic that uses the global estimation.

In third place, we provide a fully distributed detector that can be efficiently computed using a spatial averaging algorithm where the communication between sensors is done locally and where the required prior knowledge at each sensor about the network connectivity is minimal. Its performance is evaluated using numerical simulations showing excellent results for a wide range of signal-to-noise ratios values.

I-B Organization

The paper is organized as follows. We present the detection problem and compute the statistics of the measurements taken by the nodes in Section II. In Section III, we first propose to estimate the unknown network parameters at each node locally, and then, we compute the asymptotic distribution of this statistic under each hypothesis, which allows to characterize its asymptotic performance. In Section IV, we simplify this detector to one that can be efficiently implemented in WSNs via a consensus algorithm. In Section V, we evaluate the performance of the algorithm numerically and finally, in Section VI, we draw the main conclusions of this work. The proofs of some of the presented mathematical results are relegated to the appendices.

I-C Notation

We will denote with $\boldsymbol{1}_{N}$ the $N$-dimensional vector with all its entries equal to one, with $\mathbf{0}_{N}$ the $N$-dimensional null vector and with $\mathbf{I}_{N}$ the $N$-dimensional identity matrix. Given a vector $\mathbf{a}$ we denote with $\text{diag}(\mathbf{a})$ a diagonal matrix with diagonal entries given by the components of $\mathbf{a}$. Similarly, given a square matrix $\mathbf{A}$, we denote with $\text{diag}(\mathbf{A})$ a diagonal matrix which preserves the diagonal of $\mathbf{A}$. With $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ we denote a multivariate normal distribution with mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. Given two $N$-dimensional vectors $\mathbf{a}$ and $\mathbf{b}$ we write $\mathbf{a}\succeq\mathbf{b}$ ($\mathbf{a}\succ\mathbf{b}$) if $a_{i}\geq b_{i}$ ($a_{i}>b_{i}$) for all $i\in[1:N]$.

II Model

We consider a WSN with $N$ nodes with sensing capabilities distributed in a bounded geographical area. Each sensor position is denoted by $\mathbf{x}_{k}\in\mathbb{R}^{2}$ with $k\in[1:N]$. We will assume that at an unknown position $\mathbf{x}_{0}$ there is a possibility of having a source emitting a signal $s(t)$. Each sensor has a observation window of duration $\tau$ in which observes a signal $y_{k}(t)$, $t\in[0,\tau]$ $k\in[1:N]$. Through the processing of their observations the network looks for the correct decision regarding the presence or absence of the source in a fully distributed manner. This means that each sensor node has to reach the same decision about the presence or absence of the source without the help of a FC. This leads us to the following binary hypothesis testing problem

[TABLE]

where $v_{k}(t)$ is a zero-mean base-band Gaussian complex circular111We will work with low-pass complex equivalent signals. sensing noise with flat spectra with value $N_{0}$ and independent through the sensors. The source signal $s(t)$ is also assumed to be a zero-mean base-band Gaussian complex circular stationary stochastic process, independent from the sensing noise signals $v_{k}(t),\ k\in[1:N]$. It also assumed that the spectrum of $s(t)$ is again flat with value $N_{s}$ (this model can be generalized in several ways, see Remark 1 at the end of the section). We consider that the sensing system at each node has a limited two-sided bandwidth of $2W$, which leads us to stochastic signals with limited bandwidth under both hypothesis. The value of $h_{k}$ is assumed to be constant during the whole observation window and takes into account the characteristic of the wireless path between the source position $\mathbf{x}_{0}$ and the $k$-th sensor one, $\mathbf{x}_{k}$. It is in general a complex value which models attenuation and delay between the source position and node node $k$. For example, if we assume that a power-law attenuation is valid, then $|h_{k}|=\frac{1}{\epsilon+\|\mathbf{x}_{k}-\mathbf{x}_{0}\|^{\frac{\alpha}{2}}}$ where $\alpha$ is the path-loss exponent and $\epsilon$ is small constant. See in Fig. 1 a sketch of the network model and the sensing node scheme.

We will consider that each sensor has an energy detector being able to compute the energy of the received signals over some time window. The use of the energy detector is justified by practical considerations in a distributed setting, where a coherent detector at each sensor will require a very precise network-wide synchronization in order to take advantage of the phase or delay information of the measurements. Although some basic synchronization is always needed, energy detectors do not need such a fine synchronization as in the coherent case which could be expensive in terms of resources and difficult to achieve, specially in settings where the network have a large number of nodes. In more precise terms, we assume that the full observation window of length $\tau$ can be divided in $L=\lfloor\frac{\tau}{T}\rfloor$ time slots with duration $T$. In each slot, each sensor uses its energy detector for computing:

[TABLE]

As the signals under both hypotheses are bandwidth-limited and with flat spectrum we can use the Karhunen-Loève (KL) expansion for band-limited processes [33]. More precisely, $y_{k}(t)=\sum_{p=1}^{\infty}a_{l,p,k}\phi_{p}(t-(l-1)T)$ under $\mathcal{H}_{0}$, and $y_{k}(t)=\sum_{p=1}^{\infty}(a_{l,p,k}+s_{l,p})\phi_{p}(t-(l-1)T)$ under $\mathcal{H}_{1}$, for every $k\in[1:N]$, and for $t\in[(l-1)T,lT]$, where the eigenfunctions $\left\{\phi_{p}(t)\right\}_{p=1}^{\infty}$ are orthonormal and complete in $L_{2}\left([0,T]\right)$, the space of square-integrable functions, and satisfy the integral equation:

[TABLE]

where for each $p\in\mathbb{N}$, the eigenvalues satisfies $\lambda_{p}>0$. The quantities $a_{l,p,k}$ and $s_{l,p}$ are complex and circular Gaussian independent random variables with zero mean and second order moment equal to $2\lambda_{p}N_{0}W$ and $2\lambda_{p}N_{s}W$ respectively. If the time-bandwidth product verifies $2TW\equiv M\gg 1$ we can mimic the ideas in [34], and use the orthonormality of $\left\{\phi_{p}(t)\right\}_{p=1}^{\infty}$ to show that, for $k\in[1:N]$ and $\ l\in[1:L]$:

[TABLE]

It is well-known [35, 36] that when $M\gg 1$, $\lambda_{p}\approx 1$ when $p\in[1:M]$ and $\lambda_{p}\approx 0$ when $p>M$ which implies that $\mathbb{E}\left[|a_{l,p,k}|^{2}\right]\approx 2N_{0}W\equiv\sigma_{v}^{2}$ and $\mathbb{E}\left[|s_{l,p}|^{2}\right]\approx 2N_{s}W\equiv\sigma_{s}^{2}$ for all $k\in[1:N],\ l\in[1:L],\ p\in[1:M]$.

Using the above detailed properties about the random coefficients, it is straightforward to show that under $\mathcal{H}_{0}$, $\mathbb{E}\left[\tilde{z}_{k}(l)\right]=\sigma_{v}^{2}$ and $\text{Cov}\left[\tilde{z}_{k}(l),\tilde{z}_{n}(l)\right]=\frac{\sigma_{v}^{4}}{M}\delta_{kn}$, for all $l\in[1:L]$, where $\delta_{kn}$ is the Kronecker delta.

Similarly, under $\mathcal{H}_{1}$, we have $\mathbb{E}\left[\tilde{z}_{k}(l)\right]=\sigma_{v}^{2}(1+c_{k})$, and $\text{Cov}\left[\tilde{z}_{k}(l),\tilde{z}_{n}(l)\right]=\frac{\sigma_{v}^{4}}{M}(c_{k}c_{n}+(2c_{k}+1)\delta_{kn})$, for all $l\in[1:L]$ and $k,n\in[1:N]$, where

[TABLE]

We can group the measurements of all sensors in the time window $l$ in the vector $\tilde{\boldsymbol{z}}_{l}\equiv[\tilde{z}_{1}(l),\dots,\tilde{z}_{N}(l)]^{T}\in\mathbb{R}^{N}$ and define $\boldsymbol{c}=[c_{1},\dots,c_{N}]^{T}\succ\boldsymbol{0}_{N}$. At this point, and assuming again that $M$ is sufficiently large222We also assume that the temporal correlation between $s_{l,p}$ and $a_{l,p,k}$ in different time slots decreases sufficiently fast (mixing property [37])., we use the multidimensional Central Limit Theorem (CLT) [37] to show that

[TABLE]

where for easy reference we defined $\tilde{\boldsymbol{\mu}}(\mathbf{c})\equiv\sigma_{v}^{2}(\boldsymbol{c}+\boldsymbol{1}_{N})$ and $\tilde{\mathbf{\Sigma}}(\boldsymbol{c})\equiv\frac{\sigma^{4}_{v}}{M}\left(\boldsymbol{c}\boldsymbol{c}^{T}+2\text{diag}(\boldsymbol{c})+\mathbf{I}_{N}\right)$. For further notational simplicity we will work on the following equivalent test obtained from (4) using the change of variables: $\boldsymbol{z}_{l}=(\tilde{\boldsymbol{z}}_{l}-\sigma_{v}^{2}\boldsymbol{1}_{N})/(\sigma^{2}_{v}/\sqrt{M})$.

[TABLE]

where

[TABLE]

We will consider that the noise variance $\sigma_{v}^{2}$ is known (or estimated) at each node, and that $\boldsymbol{c}$ is unknown. This is because of the lack of knowledge of the true position of the source, but also due to the fact that the exact nature of the wireless path between each sensor is not known exactly (e.g. the value of the path-loss parameter $\alpha$) and it is also influenced by several complex phenomena (e.g. shadowing, fading, etc) which are difficult to know and model in advance. Therefore, we need a statistic that avoids the use of the unknown parameters or estimates them in some way. We will attack this issue in the next section.

It is important to observe that the vectors $\boldsymbol{z}_{l}$ contains the measurements taken in each sensor at the time slot $l\in[1:L]$. As we do not assume the presence of a FC and the sensors are geographically separated, we need to allow cooperation and communication between them in order to obtain a common and distributed decision about the presence or absence of the source. More precisely, we assume that sensors can communicate through error-free channels with other sensors in their neighborhood. The concept of neighborhood is naturally introduced modeling the WSN as a graph (where the edges weights can be defined taking into account the distance between the two connected nodes and the resources that each node can put on the communication) which is assumed to be connected. We will return to this in Section IV-A.

Remark 1**.**

It should be noted that the assumed hypotheses about the temporal correlation structure of the noise and the source signal can be relaxed in several aspects. In first place, there is no need to assume that the spectrum of $s(t)$ is flat. Using again the formalism of Karhunen-Loève expansion the non-flat spectrum case can also be treated. Unfortunately, the final model is slightly more complex requiring a few more unknown parameters besides $\boldsymbol{c}$. Moreover, we can abandon the Gaussianity hypothesis of $s(t)$ if we include another unknown parameter related with the fourth-order moment of the expansion coefficients of $s(t)$. We have chosen to present the more restricted Gaussian model of $s(t)$ with flat spectrum over the system bandwidth to simplify presentation. As our main goal is to exploit the spatial correlation in the measurements in a distributed setting, we have assumed the simpler temporal correlation model for the noise and source signals presented above.

III GLR with local estimation

III-A Local estimate $\boldsymbol{c}_{\rm L-MLE}$

The test in (5) is basically a composite hypothesis testing problem. In particular, it is a parameter test [38] (over $\boldsymbol{c}\succ\mathbf{0}_{N}$) because under both hypotheses the distribution is the same but with a different vector parameter: $\boldsymbol{\theta}=\boldsymbol{0}_{N}$ under $\mathcal{H}_{0}$, and $\boldsymbol{\theta}=\boldsymbol{c}$ under $\mathcal{H}_{1}$. In order to perform the test, we have to build a statistic without using the unknown parameter $\boldsymbol{c}$, like the Wald or Rao test [38, 39], or to estimate it as in the generalized likelihood ratio (GLR) test. We follow the later approach, which has some asymptotic guarantees [38, 39].

The classical GLR statistic to test the hypotheses in (5) is $T_{G}(\boldsymbol{z})\equiv\frac{p(\boldsymbol{z};\boldsymbol{\theta}=\hat{\boldsymbol{c}}_{\text{G-MLE}})}{p(\boldsymbol{z};\boldsymbol{\theta}=\boldsymbol{0}_{N})}$, where $\hat{\boldsymbol{c}}_{\text{G-MLE}}$ is the (global) maximum likelihood estimator (MLE)333We call this estimator global MLE to differentiate it from the local one $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ to be defined next. of $\boldsymbol{c}$ under $\mathcal{H}_{1}$ and $\boldsymbol{z}=\{\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{L}\}$. We see from (3) that each entry of the vector $\boldsymbol{c}$ is physically related with: i) the wireless link characteristics between the source and the corresponding sensor node, and ii) the second order statistical moment of the source. Also, from (6), there is correlation between the energy measurements taken at different nodes, given by the term $\boldsymbol{c}\boldsymbol{c}^{T}$.

It is also observed that the MLE for this model is difficult to compute even in a centralized scenario (where all measurements can be conveyed to a FC) given that the estimated value of each entry of $\boldsymbol{c}$ is a function of the energy measurements in all sensor nodes, that is, $\hat{c}_{{\rm G-MLE},k}\equiv f_{k}(\boldsymbol{\{}\boldsymbol{z}_{l}\}_{l=1}^{L})$ for each $k\in[1,N]$. In a distributed setting this would imply that each individual energy measurement in each sensor $k$ and time slot $l$ should be made available to all nodes in the network which is clearly not a practical solution. Even in the hypothetical case that all sensor measurements could be conveyed to each node across the network which would allow the computation of the MLE at each node, we would have the additional difficulty that not closed form mathematical solution for the MLE is available. Although each node could perform a numerical procedure to find the MLE, this would have a large computational load (that scales with network size $N$) which could impose a serious practical constraint, specially for nodes with limited computational capacity.

Looking for a simpler approach to the computation of the MLE and taking into account that the local energy samples at each node should be sufficiently informative about the corresponding true value of $c_{k}$ for $k\in[1:N]$, we consider a local MLE in sensor node $k$ which only use the locally sensed values $\left\{z_{k}(l)\right\}_{l=1}^{L}$. Therefore, we obtain the estimate $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ of $\boldsymbol{c}$ using the local estimates $\hat{c}_{{\rm L-MLE},k}\equiv g_{k}(\left\{z_{k}(l)\right\}_{l=1}^{L})$, $k\in[1:N]$. In more precise terms, $\hat{c}_{{\rm L-MLE},k}$, is estimated using the model (5) under $\mathcal{H}_{1}$ where sensor $k$ has access to its $L$ measurements whose distribution under $\mathcal{H}_{1}$ is $\{z_{k}(l)\}_{l=1}^{L}\overset{iid}{\sim}\mathcal{N}(\sqrt{M}c_{k},(c_{k}+1)^{2})$, that is the marginal distribution from $\mathcal{N}(\boldsymbol{\mu}(\mathbf{c}),\mathbf{\Sigma}(\mathbf{c}))$ for measurements at sensor node $k$. It is shown in the Appendix A that $\hat{c}_{{\rm L-MLE},k}$ is given by:

[TABLE]

where $m_{k}=\textstyle\frac{1}{L}\sum_{l=1}^{L}z_{k}(l)$ and $p_{k}=\textstyle\frac{1}{L}\sum_{l=1}^{L}z_{k}^{2}(l)$. It is not difficult to see that the above computed local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ is the MLE for a signal model given by $\{\boldsymbol{z}_{l}\}_{l=1}^{L}{\sim}\prod_{l=1}^{L}\mathcal{N}(\boldsymbol{\mu}(\mathbf{c}),\text{diag}(\Sigma(\boldsymbol{c})))$. That is, a model in which the correlation between the measurements at different sensor nodes is neglected. However, it is important to notice that, although we are neglecting this correlation, $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ is still an asymptotically consistent estimator of the true parameter $\mathbf{c}$. More importantly, it does not introduce a penalty in the asymptotic performance of the corresponding GLR statistic. We will show this in the next section.

III-B Asymptotic performance

In this section, we consider the asymptotic distribution, when $L\rightarrow\infty$, of the so called local GLR statistic in which we use the local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$. We also provide a comparison with the full GLR statistic $T_{G}$ which uses the true global MLE $\hat{\boldsymbol{c}}_{\text{G-MLE}}$. We prove that the local GLR statistic has exactly the same asymptotic distribution. This strongly motivates the use of the local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ which is clearly easier to compute in a distributed setting.

III-B1 Full GLR

It is well known that the distribution of the global GLR statistic $T_{G}$ is given by: [38]:

[TABLE]

where the symbol $\overset{a}{\sim}$ means “asymptotically distributed as when $L$ tends to infinity”, $\chi^{2}_{N}$ is the chi-square distribution with $N$ degrees of freedom and $\chi^{\prime 2}_{N}(\lambda_{g})$ is the non-central chi-square distribution with $N$ degrees of freedom and non-centrality parameter $\lambda_{g}=L\boldsymbol{c}^{T}\boldsymbol{i}(\boldsymbol{0})\boldsymbol{c}$, where $\boldsymbol{i}(\boldsymbol{0})$ is the Fisher information matrix [40] evaluated at $\boldsymbol{\theta}=\boldsymbol{0}$. As shown in Appendix B, $\lambda_{g}=L(M+2)\|\boldsymbol{c}\|^{2}$.

III-B2 GLR with local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$

We will obtain the asymptotic distribution of the local MLE and the local GLR statistic $T_{L}(\boldsymbol{z})$ defined by

[TABLE]

We first note the particular structure of the joint pdf $p(\boldsymbol{z}_{l};\boldsymbol{\theta})$ with $l\in[1:L]$ under $\mathcal{H}_{1}$ as shown in (5) and (6). From these equations and the discussion in the preceding section, it is clear that the marginalization of $p(\boldsymbol{z}_{l};\boldsymbol{\theta})$ over all components but the $k$-th, only depends on the $k$-th component of $\boldsymbol{\theta}$, $\theta_{k}$, and not on the whole vector $\boldsymbol{\theta}$, that is:

[TABLE]

However, the components of $\boldsymbol{z}_{l}$ are not independent, i.e., $p(\boldsymbol{z}_{l};\boldsymbol{\theta})\neq\prod_{k=1}^{N}p_{k}(z_{k}(l);\theta_{k})$. Let $\mathcal{A}$ be a set of feasible points of $\boldsymbol{\theta}$ that includes the true parameter $\boldsymbol{c}$. Define $\hat{\boldsymbol{\theta}}_{loc}^{L}\equiv\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ and consider the local MLE estimate $\hat{\boldsymbol{\theta}}_{loc}^{L}=[\hat{\theta}_{loc,1}^{L},\dots,\hat{\theta}_{loc,N}^{L}]$. It is easy to see that its components are computed solving the following problem assuming that each sensor only has access to its own measurements:

[TABLE]

where $\mathcal{A}^{k}$ is the the projection of the set $\mathcal{A}$ on the $k$-th component. The following lemma states the asymptotic distribution of $\hat{\boldsymbol{\theta}}_{loc}^{L}$ and the local GLR statistic $T_{L}$. The proof is presented in Appendix C.

Lemma 1**.**

Consider the following assumptions:

1. A1.

The first and second-order derivatives of the log-likelihood function are well defined and continuous functions. 2. A2.

$\mathbb{E}[\partial\log p_{k}(z_{k}(l);\theta_{k})/\partial\theta_{k}]=0$, $k\in[1:N]$, $\forall l$. 3. A3.

The signal is weak, i.e., $\|\boldsymbol{c}\|\!\leq\!c_{0}/\sqrt{L}$ for a constant $c_{0}$.

We have that the asymptotic distribution of $T_{L}$ is:

[TABLE]

where $\lambda_{loc}=L(M+2)\|\boldsymbol{c}\|^{2}$ is the non-centrality parameter of the non-central chi-square distribution.

Notice that this result shows us that the distribution of the global (8) and the local (11) GLR statistics are asymptotically equal and, therefore, the asymptotic performance of both statistics is the same. This implies that our distributed hypothesis testing procedure has not penalties with respect to the centralized procedure, at least asymptotically.

Remark 2**.**

Notice that the analysis of the asymptotic behavior of $\hat{\boldsymbol{\theta}}_{loc}^{L}\equiv\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ is an instance of the mismatched ML asymptotic performance problem [41], [42]. In our case, $\hat{\boldsymbol{\theta}}_{loc}^{L}$ is the so-called local MLE estimate based on the model $\prod_{k=1}^{N}p_{k}(z_{k}(l);\theta_{k})$, which is clearly different from the true data model $p(\boldsymbol{z}_{l};\boldsymbol{\theta})$. Our main interest is not in the asymptotic behavior of $\hat{\boldsymbol{\theta}}_{loc}^{L}$ but in the asymptotic behavior of the GLR statistic using this estimate. It is in this respect that we point out the relevance of Lemma 1.

IV Distributed computation of the testing statistic

In this section, we will consider the problem of adapting the GLR statistic that uses the local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ (denoted also as $\hat{\boldsymbol{c}}$ to keep notation uncluttered) in order to be efficiently computed and distributed across the network. From (5), the GLR statistic can be written as:

[TABLE]

where both $\hat{\boldsymbol{\mu}}$ and $\hat{\boldsymbol{\Sigma}}$ are built using $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ given in (7). In the $(a)$-step, we used the Woodbury matrix inversion formula and the matrix determinant lemma to compute the closed forms of $\hat{\boldsymbol{\Sigma}}^{-1}$ and $\det(\hat{\boldsymbol{\Sigma}})$, respectively. We also defined $\bar{c}_{1}\equiv\sum_{k=1}^{N}\frac{\hat{c}_{k}^{2}}{1+2\hat{c}_{k}}$, $\bar{c}_{2}\equiv\sum_{k=1}^{N}\log(1+2\hat{c}_{k})$.

In first place, notice that the computation of $\bar{c}_{1}$ and $\bar{c}_{2}$ requires the spatial sum (through the index $k$) over the sensors of the quantities $\frac{\hat{c}_{k}^{2}}{1+2\hat{c}_{k}}$ and $\log(1+2\hat{c}_{k})$ (which can be computed at each sensor using the local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$). This spatial sum (proportional to the spatial averaging of the same quantities) can be computed with algorithms already developed in the literature [43, 17, 44] or some of their variants, where the spatial average of a quantity is obtained by propagating local neighborhood averages computed at each node. We formalize this in the next section.

The second term in (12) can also be written as a spatial sum of terms that can be computed locally at each sensor node. To show this, define for each $k\in[1:N]$:

[TABLE]

As these terms can be computed locally at each sensor, the second term in (12) can be written as $\frac{1}{2}\sum_{k=1}^{N}u_{k}$, and again we can resort to a spatial averaging algorithm for its computation.

The third term in (12), however is more complicated. As we can see, the summation in time and space do not commute in general. This implies that if we want to exactly compute this term we need to implement $L$ runs of the spatial averaging algorithm. When $L$ is large this could be inefficient in terms of energy, delay and bandwidth. Therefore, in order to simplify the computation of the distributed algorithm we do commute the summations and replace $\frac{1}{L}\sum_{l=1}^{L}\left\{\sum_{k=1}^{N}(z_{kl}-\sqrt{M}\hat{c}_{k})\frac{\hat{c}_{k}}{1+2\hat{c}_{k}}\right\}^{2}$ by

[TABLE]

where $m_{k}$ was defined in Section III-A and where the term inside the square can be computed as an spatial sum of terms that can be locally computed at each sensor. We will see in Section V that this replacement does not introduce a severe penalty in the algorithm performance in a wide range of signal-to-noise ratios. Defining for each $k\in[1:N]$: $w_{k}\equiv(m_{k}-\sqrt{M}\hat{c}_{k})\frac{\hat{c}_{k}}{1+2\hat{c}_{k}},$ $\bar{u}\equiv\sum_{k=1}^{N}u_{k}$ and $\bar{w}\equiv\sum_{k=1}^{N}w_{k}$, we can write the new statistic (which we call the fully distributed statistic) as:

[TABLE]

IV-A Spatial averaging algorithm

The above distributed statistic requires the computation of quantities $\bar{c}_{1}$, $\bar{c}_{2}$, $\bar{u}$ and $\bar{w}$ which as explained above are spatial sums over the different sensors in the network. Next we will generically refer to the sum $\bar{a}\equiv\sum_{k=1}^{N}a_{k}$, which will represent $\bar{c}_{1},\bar{c}_{2},\bar{u}$ or $\bar{w}$ accordingly. Each sensor node accesses to only a scalar value $a_{k}\in\mathbb{R}$, $k\in\mathcal{N}\equiv[1:N]$ and it is desired to compute the average $\tilde{a}=\frac{1}{N}\sum_{k=1}^{N}a_{k}$ (or the sum $\bar{a}=N\tilde{a}$) at each node in a distributed manner and with minimal resources allocated to the exchanges between the nodes.

Assuming that the nodes only communicates with their neighbors through error-free channels, the spatial averages can be computed via a distributed diffusive procedure such as in [43, 17, 44]. Between all the existing possibilities, we will consider a simple but effective algorithm usually called local-degree weights distributed averaging algorithm [43]. In more precise terms, consider a network (modeled as a connected graph) $\mathcal{G}=(\mathcal{N},\mathcal{E})$ consisting of a set of nodes $\mathcal{N}$ and a set of edges $\mathcal{E}$, where each edge $\{i,j\}\in\mathcal{E}$ is an unordered pair of distinct nodes. The set of neighbors of node $i$ is denoted by $\mathcal{N}_{i}=\{j\in\mathcal{N}|\{i,j\}\in\mathcal{E}\}$. Notice that the sensor $i\notin\mathcal{N}_{i}$. See Fig. 1 (a) for a graphical description.

The average value $\tilde{a}$ can be computed iteratively as:

[TABLE]

where $a_{k}(t)$ is the average after $t$ iterations (or message exchanges between the nodes), $a_{k}(0)=a_{k}$ is the initial value and $W_{kj}$ is the weight on $a_{j}(t-1)$ at the node $k$. These set of equations can be succinctly written using using a matrix formulation. To this purpose, let $\boldsymbol{a}(t)\equiv[a_{1}(t),\dots,a_{N}(t)]^{T}$. Then, the iterative equation in its matrix form is:

[TABLE]

where $\mathbf{W}$ is the matrix of weights with elements $W_{kj}=(\mathbf{W})_{kj}$. Considering local communication only, i.e., each node broadcasts its local value at iteration $t$ only to the nodes in its neighborhood, we have that for each $k\in\mathcal{N}$, $W_{kj}=0$ for $j\notin\mathcal{N}_{k}$ and $j\neq k$. Thus, the feasible weight matrices must satisfy a sparsity pattern given by the network connectivity: $\mathbf{W}\in\mathcal{S}$, where $\mathcal{S}=\{\mathbf{W}\in\mathbb{R}^{N\times N}|W_{kj}=0\text{ if }\{k,j\}\notin\mathcal{E}\text{ and }k\neq j\}$.

The optimum weights matrix $\mathbf{W}$ in terms of the asymptotic convergence factor can be computed by solving semidefinite program (SDP), assuming that $\mathbf{W}$ is a symmetric matrix [43]. Although this procedure guaranties the fastest convergence of $\boldsymbol{a}(t)$ to the average vector $\tilde{a}\boldsymbol{1}_{N}$ when $t\rightarrow\infty$, each node should be aware of its corresponding weights to perform the average. This would require to solve the mentioned SDP program in a distributed manner or, optionally, in a centralized manner and then to communicate the corresponding weights to each node. In both cases, the complexity of this procedure is high.

A simpler way is to select the weights directly, without any optimization procedure. This, for example, can be achieved with local-degree weights [43] where the convergence to the required average is guaranteed given that graph is not bipartite. Although with this choice we are sacrificing speed of convergence to the desired average, we will use it because it requires a minimal knowledge, at each node, about the network topology. This is certainly a very much desired feature, specially for large and/or rapidly changing networks.

The weights are defined as follows. Assume arbitrarily a direction for each edge of the graph. Let $P\equiv|\mathcal{E}|$ be amount of edges of the graph, and define the incident matrix $\mathbf{A}\in\mathbb{R}^{N\times P}$ as444It can be proved that neither the election of the edges direction nor the sign of the ones in the incident matrix modify the weights defined in (17).

[TABLE]

Considering symmetric weights, each edge is associated with a unique weight $w_{l}=W_{kj}=W_{jk}=1/\max(d_{k},d_{j})$, where edge $l\in[1:P]$ connects nodes $k$ and $j$ and $d_{k}$ is the degree of node $k$, i.e., the number of neighbors of node $k$. Letting $\boldsymbol{w}\in\mathbb{R}^{P}$ with components $(\boldsymbol{w})_{l}=w_{l}$, the matrix of weights can be written as

[TABLE]

It should be easy to notice that, as explained above, this construction of the weights depends only of local information about network connectivity at each node (only the degree of each node is required). Clearly, in a WSN, this information is available in each node at the network layer of the communication stack.

Algorithm 1 summarizes the steps required to compute the fully distributed statistic (15). Several stopping criteria can be considered in the iterative computation of the spatial average (16). For example, we can consider stopping criteria as a fixed number of exchanges, or a fixed number of exchanges after no significant changes in each $a_{k}(t)$, $k\in[1:N]$, is observed. Although this a important aspect of the distributed calculation procedure for the spatial averaging, in this work, we will consider the former option to evaluate numerically the algorithm performance in Section V.

It is important to consider the total number of messages that the nodes need to exchange in order to compute $T_{L-FD,k}$, for each $k\in[1:N]$. Let $N_{it}$ be the predefined number of message exchanges or iterations established for computing the fully distributed statistic. As already mentioned, it needs to compute 4 spatial sums: $\bar{c}_{1,k}$, $\bar{c}_{2,k}$, $\bar{u}_{k}$ and $\bar{w}_{k}$, and for each of them, $N\times N_{it}$ message transmissions are needed, given that each node broadcasts its data to its neighborhood. Then, a total of $4NN_{it}$ transmissions are needed. On the other hand, the local statistic $T_{L}$ in (12) requires to compute $\bar{c}_{1,k}$, $\bar{c}_{2,k}$, $\bar{u}_{k}$ (as in $T_{L-FD}$) and $L\times N\times N_{it}$ broadcast transmissions to compute the last term in (12). It makes a total of $(3+L)NN_{it}$ transmissions, which is typically much greater that $4NN_{it}$ when the time slots for energy computation at each node satisfy $L\gg 1$. This analysis shows the advantage of $T_{L-FD}$ over $T_{L}$, in terms of communication resources, for the source detection problem in a distributed scenario.

It is important to remember how the use of the local MLE $\hat{\boldsymbol{c}}_{\rm{L-MLE}}$ allowed us to have a fully and efficient distributed statistic $T_{L-FD}(\boldsymbol{z})$ in terms of communication exchanges over the network. It should be also clear, that this would have been impossible with the global MLE solution $\hat{\boldsymbol{c}}_{\rm G-MLE}$. At this point we would ask ourselves if the use of the local MLE, although important from a practical point of view in the distributed setting, will bring some performance penalization of the hypothesis testing problem, in the non-asymptotic regime, with respect to the case in which the global MLE solution is employed. This will be analyzed in the following section via numerical simulations.

V Numerical results

In this section, we compare the performance of the proposed fully distributed statistic with a finite amount of samples per node against the asymptotic results. We also numerically evaluate the ability of the network to achieve consensus about the final decision (source present or absent) between its nodes.

We consider the network represented through its graph shown in the Fig. 2 with $N=10$ nodes and $|\mathcal{E}|=20$ edges, and the source located in the center $(0,0)$. This network was randomly generated following [43]. First we randomly generated 10 nodes, uniformly distributed on a square of $200\times 200\ \rm m^{2}$. We impose that two nodes are connected by an edge if their distance is less than a predefined threshold. Then we increase the threshold until the total number of edges is $20$ and check that the resulting graph is connected.

We set the following parameters which could be assumed for sensing, for example, a TV signal in the 400-800 MHz UHF band, where the bandwidth of each channel is 6 MHz [45]. Therefore, we take $W=6$ MHz555A Cognitive Radio receiver with energy detector could be implemented as in [45] using a analog-to-digital converter followed by a N-FFT operation, an averaging and a square device. In this context, the frequency bin spacing is $W=6$ MHz.. The observation window is set to $\tau=41.6\mu$s and is divided in $L=50$ time-slots of duration each $T=\frac{\tau}{L}=0.83\mu$s. Then, the time-bandwidth product $M=2WT=10$. Finally, we take the path loss to be $\alpha=4$. The non-centrality parameter $\lambda=\lambda_{loc}=\lambda_{g}$ depends on the quotient $\sigma^{2}_{s}/\sigma^{2}_{v}$ through $\boldsymbol{c}$. Thus, when the source is present, $\sigma^{2}_{s}/\sigma^{2}_{v}$ is adjusted to achieve the value $\lambda$ shown (in dB) in each figure.

Before showing the numerical results, we define the miss-detection and the false alarm probability of a statistic $T$ for a predefined threshold $\gamma$ as $P_{\text{md}}\equiv\Pr(T<\gamma|\mathcal{H}_{1})$ and $P_{\text{fa}}\equiv\Pr(T>\gamma|\mathcal{H}_{0})$. The detection probability is $P_{\text{d}}=1-P_{\text{md}}$.

In Fig. 3, we set $\lambda=12$ dB and plot several complementary receiver operating characteristics (CROC) for the presented statistics. First, we plot the (theoretical) asymptotic performance of the GLR test using global MLE estimation $\hat{\boldsymbol{c}}_{G-MLE}$, $T_{G}$, and local MLE estimation $\hat{\boldsymbol{c}}_{L-MLE}$, $T_{L}$ (c.f. (8) and (11), respectively). We also evaluate the performance of all statistics for a finite amount of measurements (i.e. $L=50$) generating $10^{4}$ Monte Carlo runs. We see that the performance of the statistic $T_{L}$ matches very well with the theoretical asymptotic performance. We also see that the fully distributed statistic $T_{L-FD}$ with $N_{it}=20$ iterations (see Algorithm 1) has a similar performance to $T_{L}$ (the curves are superposed). This shows that the replacement in (14) works well, allowing for a fully distributed computation without introducing any significant loss in the performance. Finally, the behavior of the optimal likelihood ratio (LR) test, computed through the method of Monte Carlo, is also plotted only to have a purely theoretical reference. It is not possible to implement this test in practice, due to the fact that it requires the exact knowledge of the parameters under $\mathcal{H}_{1}$.

In Fig. 4, we plot the detection probability for a wide range of $\lambda$ for a fixed false alarm probability $P_{\text{fa}}=10^{-2}$. As in the previous figure, the curves match very well in all the range validating again the performance of $T_{L-FD}$. Again, $N_{it}=20$ iterations through the nodes are used to compute the fully distributed algorithm. We emphasize that the fully distributed algorithm has almost the same performance than the global GLR statistic and has a loss of about 3 dB in the parameter $\lambda$ for $P_{\mathrm{d}}=0.9$ with respect to the unrealizable likelihood ratio test. Notice also that the parameter $\lambda$ determines the performance of the tests and can be written as $\lambda=L(M+2)N\rho_{\text{avg}}$, where $\rho_{\text{avg}}\equiv\|\boldsymbol{c}\|^{2}/N$ can be interpreted as a signal-to-noise ratio averaged through the sensor nodes. In general, $\rho_{\text{avg}}$ depends on the coverage range of the network, the signal propagation model, the source power and the noise power, and cannot be chosen freely. However, the designer of the sensor network has freedom to select the remaining parameters: $L$ (number of sensing time slots), $M$ (time-bandwidth product) and $N$ (number of sensors). Thus, the analytical characterization of the network performance obtained in this work gives a good starting point for designing wireless sensor networks for the application of source detection without the need of a FC and using simple and cheap energy detectors at each sensor node.

Finally, we evaluate the ability of the network to achieve consensus between its nodes using the fully distributed statistic. As each sensor computes, in principle and for a finite number of iterations (or communication exchanges with its neighbors), a different statistic, the decision at each node could be not the same. In order to quantify this, we define the probability of consensus of the network as the probability that all sensors make the same decision about the presence or absence of the source. Obviously, this probability depends strongly on the graph of the network. In Fig. 5, we plot this probability computed through the method of Monte Carlo against the number of iterations used to build the fully distributed statistic. It can be seen that for the network of Fig. 2, only 4 iterations (or communication exchanges at each node) are needed to achieve a 90% probability of consensus and that for 10 iterations the consensus is practically a certain event. This allows us to conclude that the spatial average procedure proposed is economical in terms of resources consumption.

VI Conclusions

In this paper we derived a signal model for nodes in a WSN that implement energy detectors to perform a test and decide if a source emitting a stochastic signal is present or not. We then built a cooperatively fully distributed algorithm suitable for a WSN and analyzed its performance. We showed that even though the measurements from different sensors are correlated, the parameters of the test can be estimated locally without asymptotically performance loss. This was analytically proved by computing the asymptotic distribution of the GLR test with the proposed local MLE estimator of the unknown parameters, showing that it is the same as the global GLR distribution which consider the global MLE estimation of the same quantities. This allows to quantify the performance of the proposed test and evaluate its dependence with the main parameters of the problem to design a WSN for this kind of applications. Finally, we showed that few iterations of the fully distributed algorithm are sufficient to achieve a high probability of consensus between the nodes, showing the usefulness of the proposed algorithm in a WSN with limited resources.

Appendix A Derivation of the local MLE

Consider sensor $k$ and its energy measurements denoted by $\boldsymbol{z}^{k}=[z_{k}(1),\dots,z_{k}(L)]^{T}$. The distribution of $\boldsymbol{z}^{k}$ under $\mathcal{H}_{1}$ is $\boldsymbol{z}^{k}\sim\mathcal{N}(\sqrt{M}c_{k}\boldsymbol{1}_{L},(c_{k}+1)^{2}\mathbf{I}_{L})$. The log-likelihood is given by (neglecting terms which do not depend on $c_{k}$):

[TABLE]

Deriving and setting to zero it is easy to show that the optimal solution for the maximum of the log-likelihood has to satisfy the following quadratic equation:

[TABLE]

Keeping the positive root, we obtain the MLE estimator (7).

Appendix B Computation of $\lambda_{g}$

To evaluate $\lambda_{g}$ we need to compute the following Fisher information matrix evaluated at $\boldsymbol{\theta}=\boldsymbol{0}$

[TABLE]

where the expectation is taken with respect to $p(\boldsymbol{z}_{l};\boldsymbol{\theta}=\boldsymbol{c})$. When $p(\boldsymbol{z}_{l};\boldsymbol{\theta})$ with $l\in[1:L]$ is Gaussian, as in (5), the $(i,j)$-th component of the Fisher information matrix can be computed as [40]

[TABLE]

From (6), we have $\frac{\partial\boldsymbol{\mu}(\boldsymbol{\theta})}{\partial\theta_{k}}=\sqrt{M}\boldsymbol{e}_{k}$, and $\frac{\partial\mathbf{\Sigma}(\boldsymbol{\theta})}{\partial\theta_{k}}=(\boldsymbol{\theta}\boldsymbol{e}_{k}^{T}+\boldsymbol{e}_{k}\boldsymbol{\theta}^{T})+2\boldsymbol{e}_{k}\boldsymbol{e}_{k}^{T}$, where $\boldsymbol{e}_{k}\in\mathbb{R}^{N}$ is the canonical vector with value $1$ in the $k$-th coordinate and [math] otherwise. Then, $\boldsymbol{i}(\boldsymbol{0})=(M+2)I_{N}$ and $\lambda_{g}=L(M+2)\|\boldsymbol{c}\|^{2}$.

Appendix C Proof of Lemma 1

The following proof follows the spirit of that one found in [38, 40], although it has important modifications. In the following we will make use of the following result:

Theorem 1** (Mean value theorem).**

[46, Th. 12.9] Let $B$ be an open subset of $\mathbb{R}^{N}$ and assume that $\boldsymbol{f}:B\rightarrow\mathbb{R}^{N}$ is differentiable at each point of $B$. Let $\boldsymbol{x}$ and $\boldsymbol{y}$ be two points in $B$ such that the segment $S(\boldsymbol{x},\boldsymbol{y})=\{t\boldsymbol{x}+(1-t)\boldsymbol{y}:t\in[0,1]\}\in B$. Then for every vector $\boldsymbol{a}$ in $\mathbb{R}^{N}$ there is a point $\boldsymbol{z}\in S(\boldsymbol{x},\boldsymbol{y})$ such that

[TABLE]

where $J(\boldsymbol{w})$ is the Jacobian matrix of $\boldsymbol{f}$ evaluated in $\boldsymbol{w}$, i.e., $[J]_{i,j}=\frac{\partial f_{i}}{\partial x_{j}}$, where $\{f_{i}\}$ are the components of $\boldsymbol{f}$. It is important to remark that $\boldsymbol{w}$ depends on $\boldsymbol{a}$.

C-A Asymptotic distribution of the local MLE

Firstly, we prove the consistency of the local MLE. Consider $k\in[1:N]$. Given that $\{z_{k}(l)\}_{l=1}^{L}$ are i.i.d. under both $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$, we can apply the weak law of large numbers (LLN):

[TABLE]

where the sum converges in probability and $\theta_{k}^{*}$ is the true parameter of the distribution ($\theta_{k}^{*}=0$ under $\mathcal{H}_{0}$ and $\theta_{k}^{*}=c_{k}$ under $\mathcal{H}_{1}$). Consider now the following inequality derived from the non-negative property of the Kullback-Leibler divergence [47] between $p_{k}(u;\theta_{1})$ and $p_{k}(u;\theta_{2})$ for arbitrary $\theta_{1}$ and $\theta_{2}$:

[TABLE]

Then, it is clear that (23) is maximized for $\theta=\theta_{k}^{*}$, and by a suitable continuity argument [48], the LHS of (22) must also be maximized for $\theta=\theta_{k}^{*}$. Therefore, $\hat{\theta}_{loc,k}^{L}\overset{p}{\rightarrow}\theta_{k}^{*}$ as $L\rightarrow\infty$ $\forall k=[1:N]$. In vector form we have $\hat{\boldsymbol{\theta}}_{loc}^{L}\overset{p}{\rightarrow}\boldsymbol{\theta}^{*}$ as $L\rightarrow\infty$. Thus, the local MLE estimator is consistent.

Now, we derive the asymptotic distribution of $\hat{\boldsymbol{\theta}}_{loc}^{L}$. Let

[TABLE]

By definition, the MLE must satisfy

[TABLE]

Consider Theorem 21 with $\boldsymbol{f}(\boldsymbol{\theta})=\frac{1}{L}\sum_{l=1}^{L}\boldsymbol{\psi}(\boldsymbol{z}_{l};\boldsymbol{\theta})$, $\boldsymbol{x}=\boldsymbol{\theta}^{*}$ and $\boldsymbol{y}=\hat{\boldsymbol{\theta}}_{loc}^{L}$, then we have

[TABLE]

where $\boldsymbol{w}^{L}$ belongs to the segment $S(\hat{\boldsymbol{\theta}}^{L}_{loc},\boldsymbol{\theta}^{*})$ and depends on $\boldsymbol{a}$. Assuming that the Jacobian matrix of $\boldsymbol{f}(\boldsymbol{\theta})$, $J(\boldsymbol{w}^{L})$ is invertible and using (24) in the previous equation,

[TABLE]

By consistency of the estimator $\hat{\boldsymbol{\theta}}^{L}_{loc}$, the segment $S(\hat{\boldsymbol{\theta}}^{L}_{loc},\boldsymbol{\theta}^{*})$ becomes the point $\boldsymbol{\theta}^{*}$ and $\boldsymbol{w}^{L}\overset{p}{\rightarrow}\boldsymbol{\theta}^{*}$ as $L\rightarrow\infty$. Thus, the expression inside the parenthesis in (25) becomes independent of $\boldsymbol{a}$ as $L\rightarrow\infty$, and therefore, it must converge in probability to $\boldsymbol{0}$. Then, using the continuity of the second-order partial derivatives of the log-likelihood function to apply the Continuous Mapping Theorem (CMT)[48], we have $J(\boldsymbol{w}^{L})\overset{p}{\rightarrow}J(\boldsymbol{\theta}^{*})$ where

[TABLE]

Clearly, $J(\boldsymbol{\theta}^{*})$ is a diagonal matrix given the fact that $p_{i}(z_{i}(l);\theta_{i})$ is a function only on $\theta_{i}$. Additionally, by the central limit theorem (CLT)

[TABLE]

where the mean of the Gaussian distribution is $\boldsymbol{0}$ by assumption A2 and its covariance matrix is defined by the local Fisher information matrix, given its resemblance to the Fisher information matrix (19)666Notice that the difference between the Fisher information matrix and the local one is the probability distribution inside the expectation: the full joint distribution appears in the former, while the marginal distribution corresponding to each node are considered in the local case.

[TABLE]

Using assumption A2 it easy to prove that $[J(\boldsymbol{\theta}^{*})]_{ii}=\mathbb{E}_{\theta_{i}^{*}}\left[\frac{\partial^{2}\log p_{i}}{\partial\theta_{i}^{2}}\right]=-\mathbb{E}_{\theta_{i}^{*}}\left[\left(\frac{\partial\log p_{i}}{\partial\theta_{i}}\right)^{2}\right]=-[\tilde{\boldsymbol{i}}(\boldsymbol{\theta}^{*})]_{ii}$. Finally,

[TABLE]

Solving for $\hat{\boldsymbol{\theta}}^{L}_{loc}$ we obtain the following result

[TABLE]

We now evaluate (26) for the model (5). The marginal distribution under $\mathcal{H}_{1}$ is $p_{i}(z_{i}(l);c_{i})=\mathcal{N}(\sqrt{M}c_{i},(1+c_{i})^{2})$. Given that the argument of the expectation only depends on $z_{i}(l)$ and $z_{j}(l)$, the expectation can be computed with respect to the jointly distribution on this pair of random variables

[TABLE]

Finally, making straightforward calculations we arrive to

[TABLE]

where $\sigma_{k}=1+c_{k}$ and $\sigma_{kj}=c_{k}c_{j}+(1+2c_{k})\delta_{kj}$.

C-B Asymptotic distribution of the local GLR

Next we start with the proof of the asymptotic distribution of the local GLR. First, we recall that the global MLE attains asymptotically the Cramer-Rao bound, i.e., it is asymptotically efficient, and therefore it satisfies:

[TABLE]

where $\hat{\boldsymbol{\theta}}^{L}$ is the global MLE. We know that this estimator is consistent, i.e. $\hat{\boldsymbol{\theta}}^{L}\overset{p}{\rightarrow}\boldsymbol{\theta}^{*}$ as $L\rightarrow\infty$. As the local MLE is also consistent, we have that $\hat{\boldsymbol{\theta}}_{loc}^{L}\overset{p}{\rightarrow}\hat{\boldsymbol{\theta}}^{L}$ as $L\rightarrow\infty$. Thus, $(\ref{eq:CReq})$ is also satisfied with $\hat{\boldsymbol{\theta}}_{loc}^{L}$ instead of $\hat{\boldsymbol{\theta}}^{L}$ when $L\rightarrow\infty$. Then, using a first-order Taylor expansion of $\boldsymbol{i}(\boldsymbol{\theta}^{*})$ around $\boldsymbol{\theta}^{*}$ and discarding the second order terms as $L\rightarrow\infty$, we have

[TABLE]

Integrating this equation with respect to $\boldsymbol{\theta}$ and evaluating at $\boldsymbol{\theta}=\boldsymbol{\theta}^{*}$:

[TABLE]

where the integration constant must be $c(\hat{\boldsymbol{\theta}}_{loc}^{L})=\log p(\boldsymbol{z};\hat{\boldsymbol{\theta}}_{loc}^{L})$ given that (31) is satisfied asymptotically by the consistence of $\hat{\boldsymbol{\theta}}_{loc}^{L}$ when $L\rightarrow\infty$. Therefore,

[TABLE]

Using the previous equation in the expression of local GLR, we obtain

[TABLE]

or

[TABLE]

which is the Wald test but using the local estimator instead of the global one. Using again the CMT and the continuity of the second-order partial derivatives, the following is satisfied when $L\rightarrow\infty$:

[TABLE]

Finally,

[TABLE]

where, under $\mathcal{H}_{0}$ we use (29) and the fact that $\bar{\boldsymbol{i}}(\boldsymbol{0})=\boldsymbol{i}(\boldsymbol{0})=(M+2)I_{N}$. Under $\mathcal{H}_{1}$, we define the generalized chi-square distribution $\chi^{{}^{\prime}2}_{g,N}$ with $N$ degrees of freedom as the distribution of the square norm of the Gaussian vector $\sqrt{L}\boldsymbol{i}(\boldsymbol{c})^{\frac{1}{2}}\hat{\boldsymbol{\theta}}_{loc}$ with asymptotic distribution $\mathcal{N}\left(\sqrt{L}\boldsymbol{i}(\boldsymbol{c})^{\frac{1}{2}}\boldsymbol{c},\boldsymbol{i}(\boldsymbol{c})^{\frac{1}{2}}\text{diag}(\tilde{\boldsymbol{i}}(\boldsymbol{c}))^{-1}\tilde{\boldsymbol{i}}(\boldsymbol{c})\text{diag}(\tilde{\boldsymbol{i}}(\boldsymbol{c}))^{-1}\boldsymbol{i}(\boldsymbol{c})^{\frac{1}{2}}\right)$.

So far, we have assumed nothing about the asymptotic behavior of the true parameter $\boldsymbol{c}$. In the particular case that there exists a constant $c_{0}$ such that $\|\boldsymbol{c}\|\!\leq\frac{c_{0}}{\sqrt{L}}$ (assumption A3), when $L\rightarrow\infty$ the covariance matrix of $\sqrt{L}\boldsymbol{i}(\boldsymbol{c})^{\frac{1}{2}}\hat{\boldsymbol{\theta}}_{loc}$ becomes the identity matrix $\mathbf{I}_{N}$ and $\chi^{{}^{\prime}2}_{g,N}$ becomes the non-central chi-square with $N$ degrees of freedom and non-centrality parameter $\lambda_{loc}\!=L(M+2)\|\boldsymbol{c}\|^{2}$, as is stated in (11).