Wireless Sensor Array Network DoA Estimation from Compressed Array Data via Joint Sparse Representation

Abstract

A compressive sensing joint sparse representation direction of arrival estimation (CSJSR-DoA) approach is proposed for wireless sensor array networks (WSAN). By exploiting the joint spatial and spectral correlations of acoustic sensor array data, the CSJSR-DoA approach provides reliable DoA estimation using randomly-sampled acoustic sensor data. Since random sampling is performed at remote sensor arrays, less data need to be transmitted over lossy wireless channels to the fusion center (FC), and the expensive source coding operation at sensor nodes can be avoided. To investigate the spatial sparsity, an upper bound of the coherence of incoming sensor signals is derived assuming a linear sensor array configuration. This bound provides a theoretical constraint on the angular separation of acoustic sources to ensure the spatial sparsity of the received acoustic sensor array signals. The Cram$\acute{e}$r–Rao bound of the CSJSR-DoA estimator that quantifies the theoretical DoA estimation performance is also derived. The potential performance of the CSJSR-DoA approach is validated using both simulations and field experiments on a prototype WSAN platform. Compared to existing compressive sensing-based DoA estimation methods, the CSJSR-DoA approach shows significant performance improvement.

1. Introduction

Direction of arrival (DoA) estimation using acoustic sensor arrays has attracted significant interest due to its wide applications [1,2,3,4]. Traditionally, DoA is realized using sensor arrays, such as passive towed-array sonar systems, where all sensors are wire-connected to a fusion center (FC) [5]. For wireless sensor array networks (WSAN) [6,7], arrays (locally wire/wireless-connected clusters) of dispersed sensors are deployed over large sensor fields and communicate wirelessly to the FC. However, remote sensor arrays often suffer from restricted resources, such as local power supply, local computational capacity and wireless transmission bandwidth. Therefore, data samples at remote sensor arrays need to be compressed before transmitting to the FC. On the other hand, executing sophisticated data compression algorithms on remote sensors is also restricted by local power consumption and computational ability.

To meet those resource requirements of the WSAN, the newly-emerged compressive sensing (CS) technology [8,9,10] has great potential because of its ability to reconstruct the raw signal from only a small number of random measurements. In this work, we propose a novel DoA estimation approach that accomplishes low power, robust DoA estimation over a WSAN platform, which reduces the volume of transmitted data without complicated local data compression operations. Incorporating a CS-based formulation, this approach is able to directly estimate the incident angles of acoustic sources from randomly-sampled acoustic signals. This is made possible by fully exploiting the joint spatial and spectral sparse structure of the acoustic signals acquired by the sensor array. Hence, we refer to this approach as compressive sensing joint sparse representation direction of arrival (CSJSR-DoA) estimation method.

Existing CS-based source DoA estimation methods can be categorized into two different approaches: compressive bearing estimation (COBE) [11,12,13] and compressive sensing array DoA estimation (CSA-DoA) [14]. With COBE, the incident angles (DoAs) are modeled as a sparse angle vector, with non-zero entries indicating the presence of (a few well separated) acoustic sources. The acoustic signals received at each sensor can then be represented as the product of a redundant steering matrix and a sparse angle vector. The steering matrix, consisting of steering vectors of all possible incident angles, must be estimated using the received noisy acoustic signal at a reference node. Thus, the raw data at the reference node and randomly-projected measurements of the non-reference nodes need to be transmitted to the remote FC. This imposes a heavy communication cost and excessive energy consumption at remote sensor arrays. With CSJSR-DoA, individual sensor nodes only perform data acquisition and random subsampling, while the DoA estimation is performed at the FC. Therefore, no local computation is required at the remote sensor nodes, and the CSJSR-DoA approach is more suitable for a resource-restricted WSAN system.

With the CSA-DoA [15,16,17] method, analog sensor array data are randomly projected onto a lower dimensional subspace in the analog domain before being converted into digital samples using an analog to digital converter (ADC). As such, fewer ADCs are required. This approach is similar to the CS-camera projection and was originally proposed in [18]. The steering vectors of different incident angles at a certain frequency are used to form the reconstruction dictionary. Further work extends this type of array processing method to broadband scenarios [19,20]. Since the data volume reduction is realized in the analog domain, special analog electronic circuits are required to implement the CSA-DoA approach. The proposed CSJSR-DoA approach, on the other hand, requires no special-purpose analog hardware. Some similar work [21] uses the angle domain sparsity of sources and formulates the narrow band signal of an antenna array within the BCS [22,23] framework. However, it is utilized as an alternative DoA estimation approach, and no data compression is considered.

Another distinction between the proposed CSJSR-DoA approach and these existing CS-based DoA estimation approaches is that the CSJSR-DoA approach exploits both the spatial and spectral domain structure of the acoustic sensor array signals. We argue that many practical broadband acoustic signals can be characterized by a few dominant frequency entries. With the purpose of DoA estimation, it would be sufficient to focus on these few dominant frequency entries and to exploit their frequency sparsity to enhance the performance. For example, Figure 1 shows the spectrograms of two types of broadband sources: an engine sound of a Porsche vehicle and a bird chirping [24]. It is clear that these broadband acoustic signals are dominated by multiple dominant frequency entries. The frequency domain sparse structure [25] illustrated in these figures may be utilized to realize efficient compressive sensing.

The main contributions are:

• A joint sparse representation-based DoA estimation approach is proposed, that exposes the joint spatial and spectral domain sparse structure of array signals and incorporates the multiple measurement vector [26] approach to solve the DoA estimation problem;
• A theoretical mutual coherence bound of a uniform linear sensor array is provided, which defines the minimum angular separation of the sources that are required for the CSJSR-DoA approach to yield a reliable solution with a high probability;
• The Cram$\acute{e}$r–Rao bound of the CSJSR-DoA estimator is derived to quantify the theoretical DoA estimation performance;
• A prototype acoustic WSAN platform is developed to validate the effectiveness of the proposed CSJSR-DoA approach.

The remainder of this paper is organized as follows. In Section 2, the array signal model, the sparsity model of narrowband array processing and the CS theory are briefly reviewed. The CSJSR-DoA approach is derived in Section 3. In Section 4, the theoretical analysis of the DoA estimation performance, as well as the Cram$\acute{e}$r–Rao bound are derived. In Section 5, the performance of the CSJSR-DoA approach is evaluated by simulations and field experiments using a prototype wireless sensor array platform. At last, the conclusion is summarized in Section 6.

2. Background

In this paper, the transposition and complex-conjugate transposition operations are denoted by superscripts ${(·)}^T$ and ${(·)}^H$, respectively. Lowercase and uppercase bold symbols denote vectors and matrices, respectively. The major symbols used in this paper are summarized in

Table 1. Symbols and notations. CSJSR, compressive sensing joint sparse representation.

Symbol	Explanation
H	number of sensors within an array
Q	number of sources
$R_q$	sparsity of the q-th source signal
${\mathbf{s}}_q$	the q-th source signal in the time domain
${\mathbf{r}}_q$	dominant sparse vector in the frequency domain with $R_q$-sparsity
${\mathbf{\nu }}_q$	less prominent components of ${\mathbf{s}}_q$ in the frequency domain
${\tau }_{h,q}$	time delay of the q-th source between the h-th sensor and the array centroid
$x_h\left(t\right)$	received time domain signal in the h-th array
Ψ	$N\times N$ inverse DFT matrix
${\mathbf{x}}_h$	received time domain signal at the h-th sensor, ${\mathbf{x}}_h={[x_h\left(1\right),\cdots ,x_h\left(N\right)]}^T$
${\tilde{\mathbf{d}}}_h$	signal spectrum of ${\mathbf{x}}_h$, ${\tilde{\mathbf{d}}}_h={[d_h\left(0\right),d_h\left(1\right),\cdots ,d_h(N-1)]}^T$
$\mathbf{d}\left(k\right)$	array data spectrum at the k-th frequency, $\mathbf{d}\left(k\right)={[d_1\left(k\right),d_2\left(k\right),\cdots ,d_H\left(k\right)]}^T$
${\mathbf{a}}_{\ell }\left(k\right)$	steering vector of ${\theta }_{\ell }$ at the k-th frequency,
	${\mathbf{a}}_{\ell }\left(k\right)={[exp(-j{\omega }_k{\tau }_{1,q}),\cdots ,exp(-j{\omega }_k{\tau }_{H,q})]}^T$
$\mathbf{r}\left(k\right)$	source signal vector of L directions, $\mathbf{r}\left(k\right)={[r_1\left(k\right),r_2\left(k\right),\cdots ,r_L\left(k\right)]}^T$
${\mathbf{A}}_L\left(k\right)$	steering matrix of L directions at the k-th frequency
	${\mathbf{A}}_L\left(k\right)=[{\mathbf{a}}_1\left(k\right),{\mathbf{a}}_2\left(k\right),\cdots ,{\mathbf{a}}_L\left(k\right)]$
$\mathbf{D}$	array data spectrum matrix
$\mathbf{M}$	permutation matrix
$\overline{\mathbf{d}}$	wideband array data spectrum, $\overline{\mathbf{d}}={[{\mathbf{d}}^T\left(0\right),{\mathbf{d}}^T\left(1\right),\cdots ,{\mathbf{d}}^T(N-1)]}^T$
$\tilde{\mathbf{d}}$	array spectrum of H sensors, $\tilde{\mathbf{d}}={[{\tilde{\mathbf{d}}}_1^T,{\tilde{\mathbf{d}}}_2^T,\cdots ,{\tilde{\mathbf{d}}}_H^T]}^T$
$\mathbf{u}$	sample intervals, $\mathbf{u}={[u\left(1\right),u\left(2\right),\cdots ,u\left(M\right)]}^T$
$\mathrm{r}(·)$	rounding operation
Φ	random sub-sampling matrix
${\mathsf{\Phi }}_{loss}^h$	channel loss matrix
${\mathbf{y}}_h$	received measurement of the h-th sensor in the fusion center
$\mathbf{y}$	joint measurement vector of H sensors, $\mathbf{y}={[{\mathbf{y}}_1^T,\cdots ,{\mathbf{y}}_H^T]}^T$
$\mathbf{w}$	joint noise vector of N frequencies, $\mathbf{w}={[\mathbf{w}\left(0\right),\mathbf{w}\left(1\right),\cdots ,\mathbf{w}(N-1)]}^T$
Θ	joint measurement matrix of H sensors, $\mathsf{\Theta }=\mathrm{diag}({\mathsf{\Phi }}_e^1\mathsf{\Psi },{\mathsf{\Phi }}_e^2\mathsf{\Psi },\cdots ,{\mathsf{\Phi }}_e^H\mathsf{\Psi })$
$\mathrm{diag}(.)$	block diagonal matrix operation
Υ	joint sparse matrix
$\mathbf{s}$	direction indicative vector, $\mathbf{s}\left(\ell \right)={\sum }_{n=0}^{N-1}{\mathbf{r}}_{\ell }{\left(n\right)}^2$
$\mathbf{supp}(.)$	nonzero index of a vector, $\mathrm{supp}\left(\mathbf{x}\right)=\left\{n\left|\mathbf{x}\right(n)>0\right\}$
$\widetilde{\mathsf{{\rm Y}}}$	pruned joint sparse matrix, $\widetilde{\mathsf{{\rm Y}}}=\left[\mathsf{{\rm Y}}\left[n_1\right],\mathsf{{\rm Y}}\left[n_2\right],\cdots ,\mathsf{{\rm Y}}\left[n_{\left|\mathbf{T}\right|}\right]\right],r\in \mathbf{T}$
$\mathbf{F}\left(\mathsf{\Lambda }\right)$	Fisher information matrix of parameter Λ
$\mathbf{CRB}$	Cram$\acute{e}$r–Rao bound of the CSJSR algorithm

.

2.1. Signal Model and Joint Spatial-Spectral Sparse Structure

A WSAN consists of one or more sensor arrays and one FC. Sensors within the same sensor array are arbitrarily deployed and connected to the FC via wireless channels. Sensor arrays are battery powered and have limited processing capabilities. The FC is connected to the infrastructure and has stronger processing capabilities. To conserve energy at remote sensor arrays, it is desirable to reduce the volume of data to be transmitted via wireless channels and to move as much of the processing tasks to the FC as possible.

For simplicity of notation, in this work, we assume a single sensor array consisting of H acoustic sensor nodes, and it communicates wirelessly via radio channels to a single FC. Each sensor node is equipped with an omni-directional microphone. The received acoustic signals will be sampled at a sampling frequency of $f_s$ Hz. N samples will be grouped into a frame (a snapshot) from which a DoA estimate will be made. We assume that state-of-the-art wireless synchronization protocols, such as Reference Broadcast Synchronization (RBS) [7], are applied so that clocks at sensor nodes can be synchronized at a precision in the order of micro-seconds, which is sufficiently accurate for acoustic DoA estimation.

We assume there are Q targets. The q-th source emits an acoustic signal and contains no more than $R_q$ ($N>R_q$) dominant frequency entries. Denote ${\mathbf{s}}_q$ to be a $N\times 1$ frame vector of the q-th source received at the reference node of the sensor array. Its N-point DFT (discrete Fourier transform) may be decomposed into two components:

$${\mathbf{s}}_q=\mathsf{\Psi }({\mathbf{r}}_q+{\mathbf{\nu }}_q)$$

(1)

where Ψ is an $N\times N$ inverse DFT matrix, ${\mathbf{r}}_q$ is an $N\times 1$ and $R_q$-sparse vector, containing no more than $R_q$ DFT coefficients with larger magnitude, and ${\mathbf{\nu }}_q$ contains the remaining smaller DFT coefficients.

The h-th sensor node of the same sensor array will receive the same ${\mathbf{s}}_q$ with a relative delay ${\tau }_{h,q}$. Thus, the data received at the h-th sensor are the time-delayed summation of Q sources. That is:

$$x_h\left(t\right)=\sum _{q=1}^Q{s}_q(t-{\tau }_{h,q})+{\omega }_h\left(t\right)$$

(2)

where ${\omega }_h\left(t\right)$ is the zero mean white Gaussian noise at the h-th sensor.

Figure 2 is the time delay model of a sensor array. Define ${\theta }_q$ as the incident angle of the q-th source, c as the speed of acoustic signal and $[u_h,v_h]$ as the position of the h-th sensor with respect to the array centroid (${\sum }_{h=1}^H[u_h,v_h]=[0,0]$). The time delay ${\tau }_{h,q}$ is given by:

$${\tau }_{h,q}=d_{h}cos({\theta }_q-{\rho }_h)/c=d_h\left[cos\left({\theta }_q\right)cos\left({\rho }_h\right)+sin\left({\theta }_q\right)sin\left({\rho }_h\right)\right]/c=\frac{1}{c}\left(cos\left({\theta }_q\right)u_h+sin\left({\theta }_q\right)v_h\right)$$

(3)

where ${\rho }_h$ is the direction of the h-th sensor with respect to the array centroid. After exploiting N-point DFT to ${\mathbf{x}}_h={[x_h\left(1\right),\cdots ,x_h\left(N\right)]}^T$, its frequency domain expression is given by:

$${\tilde{\mathbf{d}}}_h={\mathbf{\Psi }}^{-1}{\mathbf{x}}_h$$

(4)

where ${\tilde{\mathbf{d}}}_h={\left[d_h\left(0\right),d_h\left(1\right),\cdots ,d_h(N-1)\right]}^T$, and each component $d_h\left(k\right)$ is equivalent to multiplying the k-th entry of the acoustic spectrum ${\mathbf{r}}_q+{\mathbf{\nu }}_q$ by a phase shift $exp(-j{\omega }_k{\tau }_{h,q})$. That is:

$$d_h\left(k\right)=\sum _{q=1}^{Q}exp(-j{\omega }_k{\tau }_{h,q})(r_q\left(k\right)+{\nu }_q\left(k\right))+w_h\left(k\right)$$

(5)

where $w_h\left(k\right)$ is the frequency domain expression of ${\omega }_h\left(t\right)$. Since ${\mathbf{x}}_h$ consists of all acoustic signals, the $N\times 1$ vector ${\tilde{\mathbf{d}}}_h$ will be R-sparse ($max\left(R_q\right)\le R\le {\sum }_{q=1}^Q{R}_q$). This is because the dominant components of Q sources may fall within overlapped or different spectrum bands.

Consider the array data spectrum $\mathbf{d}\left(k\right)={[d_1\left(k\right),d_2\left(k\right),\cdots ,d_H\left(k\right)]}^T$ of H sensors at the k-th frequency; one may write:

$$\mathbf{d}\left(k\right)=\sum _{q=1}^Q{\mathbf{a}}_q\left(k\right)r_q\left(k\right)+\mathbf{w}\left(k\right)$$

(6)

where ${\mathbf{a}}_q\left(k\right)={[exp(-j{\omega }_k{\tau }_{1,q}),\cdots ,exp(-j{\omega }_k{\tau }_{H,q})]}^T$ is the steering vector of the q-th source, $\mathbf{w}\left(k\right)=[w_1^{\prime }\left(k\right),w_2^{\prime }\left(k\right),\cdots$$,w_H^{\prime }{\left(k\right)]}^T$, $w_h^{\prime }\left(k\right)={\sum }_{q=1}^{Q}exp(-j{\omega }_k{\tau }_{h,q}){\nu }_q\left(k\right)+w_h\left(k\right)$.

The relative delay ${\tau }_{h,q}$ is the time for the acoustic wave of the q-th source traveling from the array centroid to the h-th node in the array along the incident angle ${\theta }_q$. Suppose that the entire range of incident angles is divided into $L(\gg Q)$ divisions with each division corresponding to a quantized incident angle ${\theta }_{\ell }=\ell /L$, $0\le \ell \le L-1$. We can construct a redundant $H\times L$ steering matrix ${\mathbf{A}}_L\left(k\right)=[{\mathbf{a}}_1\left(k\right),{\mathbf{a}}_2\left(k\right),\cdots ,{\mathbf{a}}_L\left(k\right)]$ that includes L steering vectors. Here, ${\mathbf{a}}_{\ell }\left(k\right)={[exp(-j{\omega }_k{\tau }_{1,\ell }),exp(-j{\omega }_k{\tau }_{2,\ell }),\cdots ,exp(-j{\omega }_k{\tau }_{H,\ell })]}^T$ is defined as the steering vector corresponding to an incident angle ${\theta }_{\ell }$. Furthermore, define a $L\times 1$, Q-sparse source energy vector $\mathbf{r}\left(k\right)$ such that:

$$r_{\ell }\left(k\right)=\left\{\begin{array}{cc}{r}_q\left(k\right)\hfill & {\theta }_{\ell }\le {\theta }_q\le {\theta }_{\ell +1}\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(7)

If L is sufficiently large, it is safe to assume that different ${\theta }_q$ will fall into different angle bins. With this change of notation, Equation (6) can be rewritten as:

$$\mathbf{d}\left(k\right)={\mathbf{A}}_L\left(k\right)\mathbf{r}\left(k\right)+\mathbf{w}\left(k\right)$$

(8)

where $\mathbf{r}\left(k\right)={[r_1\left(k\right),r_2\left(k\right),\cdots ,r_L\left(k\right)]}^T$. With Equations (4) and (6), one may construct an $H\times N$ matrix:

$$\mathbf{D}=\left[\mathbf{d}\left(0\right),\mathbf{d}\left(1\right),\cdots ,\mathbf{d}(N-1)\right]={\left[{\tilde{\mathbf{d}}}_1,{\tilde{\mathbf{d}}}_2,\cdots ,{\tilde{\mathbf{d}}}_H\right]}^T$$

(9)

By performing row-vectorization of matrix $\mathbf{D}$, we have an $H·N\times 1$ vector:

$$\tilde{\mathbf{d}}={\left[{\tilde{\mathbf{d}}}_1^T,{\tilde{\mathbf{d}}}_2^T,\cdots ,{\tilde{\mathbf{d}}}_H^T\right]}^T$$

(10)

Similarly, by performing column-vectorization of matrix $\mathbf{D}$, we have another $H·N\times 1$ vector:

$$\overline{\mathbf{d}}={\left[{\mathbf{d}}^T\left(0\right),{\mathbf{d}}^T\left(1\right),\cdots ,{\mathbf{d}}^T(N-1)\right]}^T$$

(11)

Since the entries of both $\tilde{\mathbf{d}}$ and $\overline{\mathbf{d}}$ are associated with the same $\mathbf{D}$ matrix, there exists an $H·N\times H·N$ permutation matrix $\mathbf{M}$ such that:

$$\tilde{\mathbf{d}}=\mathbf{M}\overline{\mathbf{d}}$$

(12)

Equation (12) provides a link between $\tilde{\mathbf{d}}$, which exhibits the R-sparsity in the frequency domain, and $\overline{\mathbf{d}}$, which exploits the Q-sparsity in the spatial (incident angle) domain. In other words, the array samples, as summarized in matrix $\mathbf{D}$, exhibit a joint frequency and spatial (angle) domain sparsity that will be exploited later in this paper.

2.2. Compressive Sensing and Random Sub-Sampled Measurements

The presence of sparsity, as described in the models above, promises great potential for introducing CS to reduce the amount of transmitted sensor data.

The compressive sensing theory [27] states that a signal $\mathbf{x}$ may be reconstructed perfectly from a dimension reduced measurement $\mathbf{y}\in R^{M\times 1}$ through a linear measurement system $\mathbf{y}=\mathsf{\Phi }\mathbf{x}+\mathbf{n}$ when $\mathbf{\alpha }={\mathsf{\Psi }}^{-1}\mathbf{x}$ is a K-sparse vector. This is often accomplished by solving a constrained optimization problem:

$${min\left|\right|\mathbf{\alpha }\left|\right|}_{1}s.t.\left|\right|\mathbf{y}-\mathsf{\Phi }\mathsf{\Psi }\mathbf{\alpha }\left|\right|<\sigma$$

(13)

where ${\left|\right|\mathbf{\alpha }\left|\right|}_1$ is the ${\ell }_1$ norm of vector $\mathbf{\alpha }$, Φ is the measurement matrix, Ψ is the sparse matrix and σ is a noise threshold. To confirm stable reconstruction, the measurement matrix Φ is usually chosen to be a random Gaussian matrix or a random binary matrix.

The measurement vector in traditional compressive sensing problem settings is often computed as a weighted linear combination of the observed high dimensional signals. For example, in COBE [13], the $N\times 1$ sensor data vector in a frame will be multiplied to an $M\times N$ measurement matrix ($M<N$) digitally to yield a measurement vector y and then transmitted via a wireless channel to the FC. This obviously requires much computation and is likely to further deplete the energy reserve on sensor nodes. In CSA-DoA [15], an analog filter needs to be inserted before the ADC to realize the measurement operation. This approach requires special-purpose hardware and can be quite expensive.

In this work, time domain sensor data will be purposely discarded and will not be transmitted to the FC. By doing so, data reduction may be achieved without incurring expensive compression computations on the sensor nodes. Specifically, we will use a random sub-sampling matrix [28] defined on a set of non-uniform, yet known random sampling intervals $\{{\tau }_m;1\le m\le M\}$. Let a sequence be $\mathbf{u}=\left\{u\right(1),u(2),\cdots ,u(M\left)\right\}$ such that $u\left(1\right)=1$, and:

$$u\left(m\right)=u(m-1)+\mathrm{r}\left({\tau }_m\right),2\le m\le M$$

(14)

where $\mathrm{r}(·)$ is the rounding operation. Then, the $(m,n)$-th element of this proposed random subsampling matrix is given by:

$$\begin{array}{c}\mathsf{\Phi }(m,n)=\left\{\begin{array}{cc}1\hfill & 1\le n=u\left(m\right)\le N\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(15)

The selection of the proposed random subsampling matrix is because the product of $\mathsf{\Phi }\mathsf{\Psi }$ is a partial Fourier matrix, which has been proven to satisfy the restricted isometry property (RIP) property with a high probability [29].

Furthermore, in this work, the lossy nature of wireless channels will be exploited by modeling the packet loss during wireless transmission over noisy channels as a form of (involuntary) random data sub-sampling. It is well known that due to different link types (amplify forward and direct link) [30], power allocation [31] and time-varying channel conditions, wireless transmissions are likely to suffer from packet loss, and the identities of lost data packets are known at the FC. Hence, this type of packet loss can be modeled as another form of random sub-sampling. In our recent work [32], the data packet loss is modeled with a random selection matrix.

In this work, both the random sub-sampling matrix and the random selection matrix will be incorporated as a combined measurement matrix that requires no computation on the sensor nodes while achieving the desired data reduction. Denote the random sampling matrix at the h-th sensor node as ${\mathsf{\Phi }}_h$ and the random selection matrix over the wireless channels between the h-th sensor node and the FC as ${\mathsf{\Phi }}_{loss}^h$. The combined equivalent measurement matrix may then be obtained as:

$${\mathsf{\Phi }}_e^h={\mathsf{\Phi }}_{loss}^h·{\mathsf{\Phi }}_h.$$

(16)

3. CSJSR-DoA Formulation

The proposed CSJSR-DoA algorithm consists of three major steps: (1) at the h-th wireless sensor node, the received and digitized acoustic signal ${\mathbf{x}}_h$ will be sub-sampled using random compressive sampling as described in Equation (15) to yield a node measurement vector ${\mathsf{\Phi }}_h{\mathbf{x}}_h$; (2) these node measurement vectors will then be transmitted through lossy wireless channels to the FC; (3) The FC will receive an overall measurement vector $\mathbf{y}={\left[{\mathbf{y}}_1^T,{\mathbf{y}}_2^T,\cdots ,{\mathbf{y}}_H^T\right]}^T$, which will then be processed by the CSJSR-DoA approach to directly estimate a sparse DoA indicator vector. The schematic diagram of the proposed CSJSR-DOA approach is summarized in Figure 3.

Note that at the FC, the received compressive measurement of the h-th sensor node may be expressed as:

$${\mathbf{y}}_h={\mathsf{\Phi }}_{loss}·{\mathsf{\Phi }}_h{\mathbf{x}}_h={\mathsf{\Phi }}_e^h\mathsf{\Psi }{\tilde{\mathbf{d}}}_h$$

(17)

Hence, the overall measurement vector $\mathbf{y}$ may be expressed as:

$$\left[\begin{array}{c}{\mathbf{y}}_1\\ {\mathbf{y}}_2\\ ⋮\\ {\mathbf{y}}_H\end{array}\right]=\left[\begin{array}{cccc}{\mathsf{\Phi }}_e^1\mathsf{\Psi }& & & \\ & {\mathsf{\Phi }}_e^2\mathsf{\Psi }& & \\ & & \ddots & \\ & & & {\mathsf{\Phi }}_e^H\mathsf{\Psi }\end{array}\right]\left[\begin{array}{c}{\tilde{\mathbf{d}}}_1\\ {\tilde{\mathbf{d}}}_2\\ ⋮\\ {\tilde{\mathbf{d}}}_H\end{array}\right]$$

(18)

or in matrix form,

$$\mathbf{y}=\mathsf{\Theta }\tilde{\mathbf{d}}$$

(19)

where $\mathsf{\Theta }=\mathrm{diag}({\mathsf{\Phi }}_e^1\mathsf{\Psi },{\mathsf{\Phi }}_e^2\mathsf{\Psi },\cdots ,{\mathsf{\Phi }}_e^H\mathsf{\Psi })$ is the joint sparse representation matrix. As discussed earlier, $\{{\tilde{\mathbf{d}}}_h,1\le h\le H\}$ shares a common R-sparse structure, since the signal measured at one sensor node would be a time-delayed version of those at other sensors.

Next, from Equations (8) and (11), one may write:

$$\overline{\mathbf{d}}={\mathbf{A}}_L{\mathbf{r}}_L+\mathbf{w},$$

(20)

where ${\mathbf{r}}_L={[{\mathbf{r}}_L{\left(0\right)}^T,{\mathbf{r}}_L{\left(1\right)}^T,\cdots ,{\mathbf{r}}_L{(N-1)}^T]}^T$, $\mathbf{w}={[\mathbf{w}{\left(0\right)}^T,\mathbf{w}{\left(1\right)}^T,\cdots ,\mathbf{w}{(N-1)}^T]}^T$, ${\mathbf{A}}_L$ is the $HN\times NL$ joint steering matrix:

$${\mathbf{A}}_L=\left[\begin{array}{cccc}{\mathbf{A}}_L\left(0\right)& & & \\ & {\mathbf{A}}_L\left(1\right)& & \\ & & \ddots & \\ & & & {\mathbf{A}}_L(N-1)\end{array}\right]$$

Substituting Equations (12) and (20) into Equation (19) leads to:

$$\begin{array}{cc}\hfill \mathbf{y}& =\mathsf{\Theta }\tilde{\mathbf{d}}=\mathsf{\Theta }\mathbf{M}\overline{\mathbf{d}}\hfill \\ \hfill & =\mathsf{\Theta }\mathbf{M}({\mathbf{A}}_L{\mathbf{r}}_L+\mathbf{w})=\mathsf{{\rm Y}}{\mathbf{r}}_L+{\mathbf{w}}_c\hfill \\ \hfill \end{array}$$

(21)

where ${\mathbf{w}}_c=\mathsf{\Theta }\mathbf{H}\mathbf{w}$ is the subsampled Gaussian white noise received in the FC. Here, $\mathsf{{\rm Y}}=\mathsf{\Theta }\mathbf{M}{\mathbf{A}}_L$ is the joint sparse representation matrix that combines the joint sparse matrices Θ and ${\mathbf{A}}_L$. The details of this joint sparse representation matrix are:

$$\mathsf{{\rm Y}}=\underset{\mathsf{{\rm Y}}\left[1\right]}{\underbrace{\left[\begin{array}{c}{\mathbf{\varphi }}_{e,1}^1{\mathbf{a}}_{L,1}^T\left(0\right)\\ {\mathbf{\varphi }}_{e,1}^2{\mathbf{a}}_{L,2}^T\left(0\right)\\ ⋮\\ {\mathbf{\varphi }}_{e,1}^H{\mathbf{a}}_{L,H}^T\left(0\right)\end{array}\right.}}\underset{\mathsf{{\rm Y}}\left[2\right]}{\underbrace{\begin{array}{c}{\mathbf{\varphi }}_{e,2}^1{\mathbf{a}}_{L,1}^T\left(1\right)\\ {\mathbf{\varphi }}_{e,2}^2{\mathbf{a}}_{L,2}^T\left(1\right)\\ ⋮\\ {\mathbf{\varphi }}_{e,2}^H{\mathbf{a}}_{L,H}^T\left(1\right)\end{array}}}\begin{array}{c}\\ \\ \ddots \\ \end{array}\underset{\mathsf{{\rm Y}}\left[N\right]}{\underbrace{\left.\begin{array}{c}{\mathbf{\varphi }}_{e,N}^1{\mathbf{a}}_{L,1}^T\left(N-1\right)\\ {\mathbf{\varphi }}_{e,N}^2{\mathbf{a}}_{L,2}^T\left(N-1\right)\\ ⋮\\ {\mathbf{\varphi }}_{e,N}^H{\mathbf{a}}_{L,H}^T\left(N-1\right)\end{array}\right]}}$$

(22)

where ${\mathbf{\varphi }}_{e,n}^h$ is the n-th column of ${\mathsf{\Phi }}_e^h\mathsf{\Psi }$ and ${\mathbf{a}}_{L,h}\left(n\right)$ is the h-th row of ${\mathbf{A}}_L\left(n\right)$. The significance of Equation (22) is that one may reconstruct the joint sparse indicating vector ${\mathbf{r}}_L$ from the joint measurement vector $\mathbf{y}$ directly without reconstructing the raw data ${\mathbf{x}}_h,1\le h\le H$ at individual sensor nodes. ${\mathbf{r}}_L$ is a joint sparse vector with no more than R non-zero blocks, and there are at most Q non-zero entries in each of the R non-zero blocks.

With these relations, we may now formulate the joint sparse representation-based DoA estimation problem (group sparse reconstruction [33]) in a multiple measurements vector (MMV) formulation [26,34,35]:

$$\underset{{\mathbf{r}}_L}{min}\sum _{\ell =1}^L\left|\mathbf{s}\left(\ell \right)\right|\text{subject}\text{to}\left|\right|\mathbf{y}-\mathbf{{\rm Y}}{\mathbf{r}}_L{\left|\right|}_2\le \gamma$$

(23)

where $\mathbf{s}\left(\ell \right)={\sum }_{n=0}^{N-1}{\mathbf{r}}_{\ell }{\left(n\right)}^2$ and γ is a noise threshold. To solve this optimization problem, a second-order cone programming (SOCP) [36] approach has been proposed. The joint sparse representation of the wideband array signal in Equation (23) helps to improve the DoA estimation performance of wideband signals, as long as the same DoA of a signal is shared by different frequencies [37]. Note that the group sparsity-based DoA estimation for wideband signals is advantageous as compared to the traditional wideband DoA estimation methods, because it achieves near-coherent DoA estimation performance without the requirement of coherent signal projection through, e.g., the well-known focusing technique [38].

However, the computation cost is rather expensive. In [39], an MMV-Proxmethod has been proposed. We adopt this method for the underlying problem. In particular, we solve Equation (23) by a two-step process:

• Estimate the frequency domain support $\mathbf{T}$ of array signals by reconstructing the spectral sparse indicative vector $\tilde{\mathbf{d}}$ from the joint measurements $\mathbf{y}$ and prune the joint reconstruction matrix Υ by selecting $\mathsf{{\rm Y}}\left[r\right]$ with non-zero frequency bins;
• reconstruct the DoA indicative vector $\mathbf{s}$ directly from the joint measurements $\mathbf{y}$ using the pruned joint reconstruction matrix.

In summary, the first step is to estimate the non-zero block ${\mathbf{r}}_L\left(k\right)$ and to simplify the joint sparse matrix by selecting $\mathsf{{\rm Y}}\left[r\right]$ with non-zero ${\mathbf{r}}_L\left(r\right)$. The second step is to reconstruct the DoA using the pruned joint sparse representation matrix $\widetilde{\mathsf{{\rm Y}}}$. These steps are summarized in the listing of Algorithm 1.

The computational complexity of Equation (23) in the SOCP framework using an interior point implementation is $\mathcal{O}\left(L^3{N}^3\right)$ [39]. If we decouple the original problem into two sub-problems (Steps 2, 3 and 4 in the above algorithm), its computational complexity can be reduced to $\mathcal{O}(L^3{R}^3+N^3{H}^3)$.

Algorithm 1 Decoupled joint sparse reconstruction.

Input: Received joint random measurement vector $\mathbf{y}$;

Equivalent random measurement matrix ${\mathsf{\Phi }}_e^h$;

Output: DoA indicating vector $\mathbf{s}$;

1. Estimate the noise level by random sampling in source free scenario, $\gamma =\left|\right|\mathbf{y}|{|}_2^2$;

2. Estimate the common support $\mathbf{T}=\mathbf{supp}\left(\tilde{\mathbf{d}}\right)$ by solving:

$\begin{array}{c}\underset{\tilde{\mathbf{d}}}{min}{\displaystyle \sum _{n=0}^{N-1}}{\displaystyle \sum _{h=1}^H}{d}_h{\left(n\right)}^2\hfill \\ \text{subject to }\left|\right|\mathbf{y}-\mathsf{\Theta }\tilde{\mathbf{d}}{\left|\right|}_2^2\le \gamma ;\hfill \end{array}$

3. Construct the pruned joint reconstruction matrix $\widetilde{\mathsf{{\rm Y}}}$;

$\widetilde{\mathsf{{\rm Y}}}=\left[\mathsf{{\rm Y}}\left[n_1\right],\mathsf{{\rm Y}}\left[n_2\right],\cdots ,\mathsf{{\rm Y}}\left[n_{\left|\mathbf{T}\right|}\right]\right],r\in \mathbf{T}$;

4. Solve the pruned reconstruction problem:

$\begin{array}{c}\underset{{\overline{\mathbf{r}}}_L}{min}{\parallel \mathbf{s}\parallel }_1\text{subject to }\left|\right|\mathbf{y}-\widetilde{\mathsf{{\rm Y}}}\overline{\mathbf{r}}{\left|\right|}_2^2\le \gamma ,\hfill \\ \mathbf{s}\left(\ell \right)={\sum }_{n\in \mathbf{T}}{\mathbf{r}}_{\ell }{\left(n\right)}^2,\overline{\mathbf{r}}={[{\mathbf{r}}^T\left(n_1\right),\cdots ,{\mathbf{r}}^T\left(n_{\left|T\right|}\right)]}^T;\hfill \end{array}$

4. Performance Analysis

In this section, the performance of the proposed CSJSR-DoA approach is investigated. First, the relevant coherence, its upper bound and the angular separation problem are studied. Second, the Cram$\acute{e}$r–Rao bound of the proposed CSJSR-DoA approach will be discussed.

4.1. Sparse Reconstruction Analysis and Angle Separation

The existence of solution to Equation (23) depends on the property of the reconstruction dictionary Υ. Usually, the mutual coherence and the restricted isometry property (RIP) of a sparse representation matrix are used as the criteria that indicate the reconstruction performance. Consider that the verification of the RIP property requires combinational computational complexity [40]; it is preferable to use the property of coherence that is easily computable to provide more concrete recovery guarantees. In this paper, we try to analyze the coherence of the CSJSR-DoA approach.

In the CS theory, the coherence of a matrix is the worst-case linear dependence of any pair of its column vectors. Similarly, in the block sparse signal case, in which a signal is composed by several blocks and only a small number of these blocks is non-zero, the block-coherence [41] has been proposed.

Recall from Equation (21):

$$\mathsf{{\rm Y}}=\mathsf{\Theta }\mathbf{M}{\mathbf{A}}_L$$

(24)

where Θ is a block diagonal joint sparse sampling matrix, $\mathbf{M}$ is a full rank, unitary permutation matrix and ${\mathbf{A}}_L$ is a block diagonal steering matrix. To study the coherence of Υ, we need to study the block-coherence of Υ.

Theorem 1. 

Assume that H sensor nodes form a sensor array, and the received data in the fusion center are ${\mathbf{y}}_h={\mathsf{\Phi }}_e^h{\mathbf{x}}_h({\mathsf{\Phi }}_e^h\in R^{M\times N}),h=1,2,\cdots ,H$. Then, the block-coherence of the joint sparse representation matrix Υ satisfies:

$${\mu }_B\left(\mathsf{{\rm Y}}\right)\le \sqrt{\frac{N-M}{(N-1)M}}$$

(25)

Proof. 

See Appendix A. ☐

The block-coherence ${\mu }_B$ guarantees stable reconstruction of the block sparse signal ${\mathbf{r}}_L$ with a high probability. However, each block ${\mathbf{r}}_L\left(n\right)$ indicates the DoA of sources. Thus, the coherence of each sub-matrix $\mathsf{{\rm Y}}\left[n\right]$ should be considered. In this paper, we derive the coherence of $\mathsf{{\rm Y}}\left[n\right]$. This result will provide a lower bound of the source separation (in turns of DoA angle) issue to ensure a reliable DoA estimation. However, the coherence of $\mathsf{{\rm Y}}\left[n\right]$ depends on array geometry, and it is difficult to obtain an analytical bound.

In this paper, we consider the typical uniform linear array geometry and find an interesting result of the coherence bound. The derived result is similar to the array pattern [42] analysis and is validated for both narrowband and broadband acoustic sources. For broadband sources, the highest dominant frequency may be used.

Theorem 2. 

Assume that H sensor nodes form a uniformly-spaced linear array with $d(\le \lambda /2)$ being the spacing between adjacent sensors where λ is the wavelength of the acoustic signal. If the difference of the incident angles of any pair of sources satisfies:

$$\mathsf{\Delta }\theta \ge \frac{c}{Hdf}$$

(26)

then the coherence of $\mathsf{{\rm Y}}\left[n\right]$ is bounded by:

$$\mu \left(\mathsf{{\rm Y}}\left[n\right]\right)\le \frac{1}{Hsin(\pi /H)}\sqrt{\frac{N-M}{(N-1)M}}$$

(27)

Proof. 

See Appendix B. ☐

4.2. Cram$\acute{e}$r–Rao Bound of CSJSR-DoA

In this subsection, we derive the Cram$\acute{e}$r–Rao bound (CRB) of the proposed CSJSR-DoA method and study its impact on the DoA estimation accuracy due to factors such as the number of measurements, the number of sensors and the noise level ${\sigma }^2$, which is defined in Equation (13).

Recall that the CRB is defined as the inverse of the Fisher information matrix (FIM) [43],

$$\mathbf{CRB}\left(\mathsf{\Lambda }\right)={\mathbf{F}}^{-1}\left(\mathsf{\Lambda }\right)$$

(28)

where Λ is the unknown parameter set. In this paper, we only derive the CRB of a single source case. Thus, the unknown parameter is $\mathsf{\Lambda }={[\theta ,{\mathbf{r}}^T,{\mathbf{f}}^T]}^T$, where $\mathbf{f}={[f_1,f_2,\cdots ,f_R]}^T$ and $\mathbf{r}={[r_0\left(k_1\right),r_0\left(k_2\right),\cdots ,r_0\left(k_R\right)]}^T$ are the sparse frequency bins and their corresponding source spectrum, respectively.

In the case of Gaussian white noise, the FIM is:

$$\mathbf{F}\left(\mathsf{\Lambda }\right)=\frac{{\sigma }^2}{2}\mathrm{Re}\left[{\mathbf{G}}^H\mathbf{G}\right]$$

(29)

where $\mathbf{G}=\partial \mathbf{y}/\partial {\mathsf{\Lambda }}^T$. For simplification, we assume that all H sensors have the same ${\mathsf{\Phi }}_e({\mathsf{\Phi }}_e^1=,\cdots ,={\mathsf{\Phi }}_e^H={\mathsf{\Phi }}_e)$. Recall Equations (4), (6) and (17); the measurement vector of the CSJSR-DoA approach in a single source case is:

$$\mathbf{y}\left(\mathsf{\Lambda }\right)=\sum _{i=1}^R{r}_0\left(k_i\right){\mathbf{a}}_0\left(k_i\right)\otimes {\mathbf{\varphi }}_{e,k_i}+{\mathbf{w}}_c$$

(30)

where ${\mathbf{\varphi }}_{e,k_i}$ is the $k_i$-th column of ${\mathsf{\Phi }}_e\mathsf{\Psi }$. After some calculations, the CRB of the CSJSR-DoA approach is given by:

$$CRB\left(\theta \right)=\frac{{\sigma }^2}{{2M\alpha \parallel \mathbf{e}\parallel }_2^2}$$

(31)

where $\mathbf{e}={[f_1{r}_0\left(k_1\right),f_2{r}_0\left(k_2\right),\cdots ,f_R{r}_0\left(k_R\right)]}^T$, $\alpha =\frac{4{\pi }^2}{c^2}{\sum }_{h=0}^{H-1}{\left(u_{h}cos\left(\theta \right)-v_{h}sin\left(\theta \right)\right)}^2$.

Proof. 

See Appendix C. ☐

5. Performance Evaluation

In this section, extensive simulations are carried out to: (1) study the performance of the CSJSR-DoA approach; and (2) compare the performance of the CSJSR-DoA approach against COBE and CSA-DoA, using the L1-SVD [26] algorithm as the baseline. In addition, we also developed a hardware prototype platform and collected data from outdoor experiments to validate the practical applicability of the CSJSR-DoA approach.

5.1. Simulation Settings

We synthesized acoustic signals based on Equation (1). The synthesized signals contain four dominant frequency entries at 300 Hz, 500 Hz, 600 Hz and 800 Hz. These correspond to wavelengths of $1.14$, $0.7$, $0.57$ and $0.43$ m, respectively, assuming the sound speed at 343 m/s. The received acoustic signal also contains additive Gaussian white noise with zero mean. The variances are set so that the resulting SNR ranges from $-10$ dB to 10 dB in 5-dB increments. Such an acoustic signal resembles the dominant entries’ distribution of some practical acoustic sources, such as truck engine sounds. The sampling rate is 2 ksamples/s. At this rate, a frame of 125 ms is selected. The sampling rate is chosen so that it satisfies the Nyquist sampling rate and avoids the frequency aliasing issue.

We assume a single uniform linear array consisting of six ($H=6$) acoustic sensor nodes deployed in a sensor field. The spacing between adjacent sensors is 0.2 m, which is smaller than half of the wavelength ($0.43$ m) of the 800-Hz entry. All acoustic sources are located at the broad side of this linear array, so that the incident angles (directions of arrivals) are constrained within a $-{90}^{\circ }$ to ${90}^{\circ }$ range. We divide this angle range into ($L=360$) partitions, so that each partition equals $0.5^{\circ }$.

Consider that data packets transmitted from the sensor array to the FC may suffer from data packet loss. We will simulate data packet loss rates, denoted by $r_{\mathrm{loss}}$, at 0, $5\%$, $10\%$, $20\%$ and $30\%$, respectively. To mitigate burst data transmission loss, the acoustic data stream will be interleaved while being assembled into the data packet for transmission. However, other than packet headers, there will be no additional channel coding bits appended. The data loss rate is computed as the percentage of packets that are lost during transmission through wireless channels versus the total number of packets that are sent from the sensor array.

With the CSJSR-DoA approach, both random subsampling and wireless channel packet loss will reduce the amount of received acoustic samples compared to what was originally sampled at the sensor array. In particular, at each sensor node, the acoustic data samples will be randomly subsampled according to Equations (14) and (15). The ratio $M/N$ in these equations can be regarded, in the context of CS, as the ratio of the number of measurements versus the total number of data samples. With channel data packet loss taken into account, we instead define $r_{\mathrm{dc}}(\le 1)$ to be the ratio of the number of acoustic samples successfully received at the fusion center versus the total number of samples acquired by the H acoustic sensors. Ignoring the overhead of data packet header, we have:

$$r_{\mathrm{dc}}\approx (1-r_{\mathrm{loss}})\times M/N$$

(32)

In this section, the DoA performance will be reported against different values of $r_{\mathrm{dc}}$. In practice, if $r_{\mathrm{loss}}$ can be estimated in real time, we may adjust M to achieve the desired DoA accuracy.

Four algorithms are implemented: L1-SVD, CSJSR-DoA, CSA-DoA and COBE. The L1-SVD algorithm [26] is applied to the entire set of acoustic samples without random subsampling and is used as a baseline for both simulation and experiment. The COBE algorithm [13] requires a specific reference node to be sampled at the Nyquist rate without any subsampling in order to provide a reference source signal. To have fair comparisons, we performed trial runs to empirically obtain the best parameters for the CSA-DoA and COBE algorithms. All four algorithms are implemented using MATLAB Version 7.9.0. The optimization is performed using the Sedumi 1.2.1 toolbox [44].

5.2. CSJSR-DoA Performance Evaluation

In this simulation, two ($Q=2$) stationary acoustic sources are placed in the far field from the sensor array at $[-{60}^{\circ },{30}^{\circ }]$. The noise levels of the acoustic signal correspond to SNR = $-10$ dB, $-5$ dB, 0 dB, 5 dB and 10 dB, and it is assumed that $r_{\mathrm{loss}}=0\%$. Five hundred independent trials are performed. Consider that sparse reconstruction is a probability event; a DoA estimate is considered a success detection if a sparse indicative vector can be reconstructed successfully. In other words, some trials fail to give a solution to the sparse DoA indicative vector. Based on the results of detected trials, the number of impinging signals, Q, and the corresponding DoAs can be estimated. Similar to some spatial spectrum search-based approaches, the spatial spectrum will be calculated using ${\mathbf{r}}_{\mathrm{ss}}=20{log}_{10}\left(\mathbf{s}/max\left(\mathbf{s}\right)\right)$, and some conventional peak-finding approaches [45,46] are used to find $\tilde{Q}$ peaks of the obtained spatial spectrum. However, the estimated $\tilde{Q}$ directions can be divided into three categories: $\tilde{Q}<Q$, $\tilde{Q}=Q$ and $\tilde{Q}>Q$. In this case, a proper performance criterion is hard to choose. Following a similar work [22], which takes into account both the errors in estimating the signal number $\tilde{Q}$ and the corresponding DoAs, the root mean square errors (RMSE) of the n-th detected trial will also be reported using the formula:

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{n=1}^{N_e}}{\left({\mathrm{RMSE}}^{\left(n\right)}\right)}^2}{N_e}}$$

(33)

$$\mathrm{RMS}{\mathrm{E}}^{\left(n\right)}=\left\{\begin{array}{c}\sqrt{\frac{{\sum }_{q=1}^{\tilde{Q}\left(n\right)}{\left|{\theta }_q-{\tilde{\theta }}_q^{\left(n\right)}\right|}^2+\left|Q-{\tilde{Q}}^{\left(n\right)}\right|{\left(\mathsf{\Delta }{\theta }_{max}\right)}^2}{Q}},{\tilde{Q}}^{\left(n\right)}\le Q\\ \sqrt{\frac{{\sum }_{q=1}^Q{\left|{\theta }_q-{\tilde{\theta }}_q^{\left(n\right)}\right|}^2+{\sum }_{j=Q+1}^{{\tilde{Q}}^{\left(n\right)}}{\left|{\tilde{\theta }}_j^{\left(n\right)}-{\overline{\theta }}_j^{\left(n\right)}\right|}^2}{Q}},{\tilde{Q}}^{\left(n\right)}>Q\end{array}\right.$$

(34)

where $N_e$ is the number of detected trials, ${\overline{\theta }}_j^{\left(n\right)}=arg\left\{\underset{{\theta }_q,q\in [1,Q]}{min}\left|{\theta }_q-{\tilde{\theta }}_j^{\left(n\right)}\right|\right\}$. Here, $\mathsf{\Delta }{\theta }_{max}$ is a penalty term of the maximum admissible DoA error (i.e., ${180}^{\circ }$ for a liner array), and ${\overline{\theta }}_j^{\left(n\right)}$ is the least error of additional false DoA.

The averaged detection rate and the RMSE are summarized in Figure 4. It may be noted that when the SNR is greater or equal to 0 dB, $r_{\mathrm{dc}}\ge 0.25$ yields very satisfactory DoA performance. To put it into context, $r_{\mathrm{dc}}=0.25$ roughly equals transmitting 500 samples per second without data loss. The corresponding data volume of Nyquist rate sampling is 1600.

It is also interesting to examine how the RMSE obtained using simulation is compared against the theoretical Cram$\acute{e}$r–Rao bound. In Figure 5, it can be seen that for $\mathrm{SNR}>0$ dB, the RMSE is very close to the theoretically-computed CRB.

Next, we want to investigate the impact of varying the data packet loss rate for different values of random subsampling ratio $M/N$. The noise level of the acoustic signal is set at SNR = 0 dB. Two sets of plots are provided in Figure 6. On the top, the averaged probability of successful detection and the RMSE are plotted against $M/N$. The performance degradation due to data packet loss can be seen clearly. On the bottom, the averaged probability of successful detection and the RMSE are plotted against $r_{\mathrm{dc}}$, which, according to Equation (32), already takes into account the impact of $r_{\mathrm{loss}}$. Hence, the different curves on the bottom panel are merged into an identical one. This is because the random data loss is equivalent to the random subsampling of the raw data.

Our next goal is to investigate the DoA resolution of the proposed CSJSR-DoA approach when two acoustic sources become closer in terms of DoAs. Under the same simulation configuration, we fix one source at $0^{\circ }$ and vary the incident angle of the second source from $2^{\circ }$ to ${40}^{\circ }$ in $2^{\circ }$ increments. We perform 300 independent trials for each setting. The $\mathrm{SNR}=5\mathrm{dB}$, and the peaks of the $360\times 1$ sparse vector indicate the DoAs of two sources.

We report the averaged angle estimation error versus target DoA separation angles in Figure 7. Based on the discussion on Theorem 2, the steering vectors of two nearby sources are strongly correlated or have larger coherence. Therefore, the reconstruction of a correct sparse DoA indicative vector does not perform well. In simulations, the probability (false DoA probability) that only one source is found when two DoAs exist and the RMSE of each DoA estimation versus target DoA separation angles are shown in Figure 8. Note that when the separation angles are small, the averaged estimation errors, the RMSEs and the probability of false DoA estimation are rather large. However, when the targets are well separated ($\ge {20}^{\circ }$), the DoAs of both targets are accurately estimated (the averaged estimation error and the false DoA probability are close to zero, and the RMSEs decrease to a constant). This validates the theoretical angle separation result of $\frac{c}{Hdf}=\frac{343}{6\times 0.2\times 800}=20.5^{\circ }$.

Finally, we use simulation to compare the proposed CSJSR-DoA approach against two other CS-based DoA methods: COBE and CSA-DoA. The simulation conditions are identical to those used in Figure 4. While these three methods differ in how the measurements are made, the definition of $r_{\mathrm{dc}}$ would still be applicable to all of them; namely, the number of samples received at the FC versus the total amount of raw data samples at the sensor array. Specifically, the analog projection of the CSA-DoA approach is emulated by projecting the digitized signal. The results are summarized in Figure 9. It shows that the CSJSR-DoA method has favorable performance over others in a number of situations. Additional comparison of these methods using data collected from a prototype WSAN platform will be discussed next.

5.3. Prototype WSAN Platform and Field Experiment

To further validate the capability of the CSJSR-DoA approach, a WSAN [47] was developed. Each sensor node (shown in Figure 10a) is equipped with an omni-directional microphone, a 16-bit ADC and a wireless transceiver operating at a rate of 31.25 KBps using the ZigBee protocol [48]. Random sub-sampled acoustic data samples are transmitted to the FC through the wireless transceiver to the FC (laptop). To guarantee time synchronization within the sensor array, the RBS time synchronization scheme is used. Figure 10b shows the experiment configuration, in which a laptop equipped with a wireless receiver was used as the FC. The sensor array is shown in the background close to the upper right corner.

In the first experiment, the task is to track the DoAs of a single moving acoustic source (digital acoustic signal played through a speaker). The distance between the source and the array varies from 20 m to 40 m, and the measured SNR in the array was about 6 dB. Each sensor node took samples during a period of 125 ms at a 2 Ksampling rate and transmitted them to the FC using the IEEE 802.15.4 protocol. In the FC, the DoA estimate was updated every second. The measured data loss ratio $r_{\mathrm{loss}}$ was $0.27$, and the corresponding $r_{\mathrm{dc}}$ of the received data in the FC was $0.23$. The result is shown in Figure 11. The red line is the ground truth, and the green circles are CSJSR-DoA estimates. The averaged probability of success detection is $90\%$, and the corresponding RMSE is $2.{76}^{\circ }$. This experiment provides an example of the practical use of the CSJSR-DoA approach.

In the second experiment, two stationary acoustic sources are used. They are placed at distances of $12.94$ m and $14.66$ m with (ground truth) DoA angles at $-6.5^{\circ }$ and $31.0^{\circ }$ respectively. Six sensors were incorporated into the array with an inter-sensor spacing equal to 0.2 m. Using data collected in this experiment, the performance of the CSJSR-DoA approach is compared against COBE, CSA-DoA and the baseline algorithm L1-SVD. The results are shown in Figure 12 and

Table 2. Result of the prototype system experiment. CSA, compressive sensing array; COBE, compressive bearing estimation.

	L1-SVD	CSJSR-DoA	CSA-DoA	COBE
DoA (${}^{\circ }$)	[$-7$,29]	[$-7.5$,27.5]	[$-2.5$,30]	[$-9.5$,25]

. The data compression ratio of all of these algorithms is fixed at $r_{\mathrm{dc}}=0.3$. To facilitate the comparison, the corresponding $L\times 1$ sparse DoA angle vector $\mathbf{s}$ is normalized and expressed in units of dB. From Figure 12, it is shown that the CSJSR-DoA approach yields sharper DoA estimates compared to the other three methods.

6. Conclusions

In this paper, a compressive sensing joint sparse representation approach (CSJSR-DoA) is presented for DoA estimation on a WSAN platform. By exploiting both frequency domain and spatial domain sparsity, the CSJSR-DoA approach provides a direct DoA angle estimation at the FC, while requiring almost no computation at power-constrained remote sensor nodes. We provided performance analysis in terms of DoA angle resolution and the Cram$\acute{e}$r–Rao bound of the estimates. We further conducted extensive simulation and built a prototype experimental WSAN platform to investigate the impacts of various parameters on the DoA performance. We also compared the performance of the proposed CSJSR-DoA approach against two other compressive sensing DoA estimation methods and showed that the CSJSR-DoA approach provides superior performance in both simulation runs and real-world experiments.