Cramér-Rao Bounds in Classical Direction of Arrival Estimation Problems

Date: 3-22 2017

Tags: direction-of-arrival, performance bound, statistics

The problem of direction-of-arrival (DOA) estimation is to find the directions of impinging far-field narrow-band signals using a sensor array (either acoustic or electromagnetic), as illustrated in the figure below. It is a classical problem in the field of array signal processing. For a comprehensive understanding of this field, it is recommended to check Dr. Van Trees's book Optimum Array Processing^[1].

doa

For a linear array, the time-domain snapshot model can be written as

\mathbf{y}(t) = \mathbf{A}(\mathbf{\theta})\mathbf{s}(t) + \mathbf{n}(t), t = 1,2,\ldots,T,

(1)

where $\mathbf{\theta} = [\theta_1, \theta_2, \cdots, \theta_K]^T$ denotes the DOAs, $\mathbf{y}(t) \in \mathbb{C}^M$ denotes the snapshot vector, $\mathbf{A}(\mathbf{\theta}) \in \mathbb{C}^{M\times K}$ denotes the array steering vector, $\mathbf{s}(t) \in \mathbb{C}^K$ denotes the source signal, and $\mathbf{n}(t) \in \mathbb{C}^M$ denotes the additive measurement noise. For a linear array whose sensors are located at $\mathcal{D}=\{d_1, d_2, \ldots, d_M\}$ , the steering matrix can be written as

\mathbf{A}(\mathbf{\theta}) = [\mathbf{a}(\theta_1), \mathbf{a}(\theta_2), \cdots, \mathbf{a}(\theta_K)],

where

\mathbf{a}(\theta_k) = \Big[\exp\Big(-j\frac{2\pi}{\lambda}d_1 \sin\theta_k\Big), \exp\Big(-j\frac{2\pi}{\lambda}d_2 \sin\theta_k\Big), \cdots, \exp\Big(-j\frac{2\pi}{\lambda}d_M \sin\theta_k\Big)\Big]^T,

and $\lambda$ denote the wavelength of the narrow-band signals.

Example: For a uniform linear array (ULA), the sensor locations are given by $\mathcal{D}=\{0, d_0, \ldots, (M-1)d_0\}$ , where $d_0$ is usually chosen to be $\lambda/2$ to avoid grating lobes.

Given the measurement model (1) and further assumptions on the statistics of $\mathbf{s}(t)$ and $\mathbf{n}(t)$ , the problem of finding the DOAs becomes a parameter estimation problem. One straight-forward approach is to use the maximum likelihood estimator (MLE). There also exists many other DOA estimation algorithms. Classical ones include beamforming, MUltiple SIgnal Classification (MUSIC)^[2], Estimation of Signal Parameters via Rotational Invariance (ESPRIT)^[3]. Recent developed ones include sparsity based algorithms and super resolution based algorithms^[4].

In parameter estimation problems, one of the important concepts is the Cramér-Rao Bound (CRB), which gives the lower bound of the variance of unbiased estimators. The CRB can be very useful in performance analysis, as well as in assessing the statistical efficiency of unbiased estimators. Compact expressions of the CRBs for the classical DOA estimation problem based on (1) were first derived and analyzed by P. Stoica and A. Nehorai^[5]^[6]^[7]. We will review these classical results as follows.

Concerning (1), there are two commonly used data model^[7]:

Conditional/deterministic model: the source signals $\mathbf{s}(t)$ are assumed to be nonrandom, and deterministic unknown.
Unconditional/stochastic model: the source signals $\mathbf{s}(t)$ are assumed to be random and unknown.

We further make the following assumptions:

The $K$ DOAs are distinct.
The additive noise is white circularly symmetric complex Gaussian noise, and temporarily uncorrelated.
The source signal sequence $\mathbf{s}(t)$ are temporarily uncorrelated.
The source signals and the additive noise are spatially and temporarily uncorrelated.

For the unconditional model, we further assume that $\mathbf{s}(t)$ follows the circularly symmetric complex Gaussian distribution with zero mean and non-singular covariance matrix $\mathbf{P}$ .

For the conditional model,

\mathbf{y}(t) \sim \mathcal{CN}(\mathbf{A}\mathbf{x}(t), \sigma^2\mathbf{I}).

(2)

The unknown parameters are

\mathbf{\eta}_\mathrm{C} = [\theta_1, \cdots, \theta_K, \Re(x(1)), \Im(x(1)), \cdots, \Re(x(T)), \Im(x(T)), \sigma^2] \in \mathbb{R}^{K+2T+1}.

(3)

Theorem 1. Under the aforementioned assumptions, the CRB of the DOAs for the conditional model is given by^[5]:

\mathrm{CRB}_\mathbf{\theta}^\mathrm{C} = \frac{\sigma^2}{2T} \big\{\Re[\mathbf{H} \circ \hat{\mathbf{P}}^T]\big\}^{-1},

(4)

where $\circ$ denotes the Hadamard product, and

\begin{aligned} \mathbf{H} &= \dot{\mathbf{A}}^H [\mathbf{I} - \mathbf{A}(\mathbf{A}^H\mathbf{A})^{-1}\mathbf{A}^H]\dot{\mathbf{A}}, \\ \dot{\mathbf{A}} &= \begin{bmatrix} \frac{\partial\mathbf{a(\theta_1)}}{\partial\theta_1} & \frac{\partial\mathbf{a(\theta_2)}}{\partial\theta_2} & \cdots & \frac{\partial\mathbf{a(\theta_K)}}{\partial\theta_K} \end{bmatrix}, \\ \hat{\mathbf{P}} &= \frac{1}{T}\sum_{t=1}^T \mathbf{x}(t)\mathbf{x}^H(t). \end{aligned}

For the unconditional model,

\mathbf{y}(t) \sim \mathcal{CN}(\mathbf{0}, \mathbf{A}\mathbf{P}\mathbf{A}^H + \sigma^2\mathbf{I})

(5)

The unknown parameters are

\mathbf{\eta}_\mathrm{U} = [\theta_1, \cdots, \theta_K, \cdots, \Re(P_{ij}), \Im(P_{ij}), \cdots, \sigma^2] \in \mathbb{R}^{K+K^2+1}.

(6)

Theorem 2. Under the aforementioned assumptions, the CRB of the DOAs for the unconditional model is given by^[7]:

\mathrm{CRB}_\mathbf{\theta}^\mathrm{U} = \frac{\sigma^2}{2T} \big\{\Re[\mathbf{H} \circ (\mathbf{P}\mathbf{A}^H\mathbf{R}^{-1}\mathbf{A}\mathbf{P})^T]\big\}^{-1},

(7)

where $\mathbf{R}=\mathbf{A}\mathbf{P}\mathbf{A}^H + \sigma^2\mathbf{I}$ , and $\mathbf{H}$ follows the same definition as in Theorem 1.

Note: In the original paper, (7) is obtained from the asymptotic ( $T\gg 0$ ) covariance matrix of the unconditional MLE. P. Stoica et. al. also gave a textbook derivation later on^[8].

Recently, there is a renewed interest in sparse linear arrays (thinned ULAs, such as minimum redundancy linear arrays^[9]^[10]) with the development of nested arrays^[11] and co-prime arrays^[12]. When the sources are uncorrelated, these arrays have the capability of identifying more sources than the the number of physical sensors ( $K > M$ ) under the unconditional model. However, careful inspection of (7) reveals that $\mathrm{CRB}^\mathrm{U}_\mathbf{\theta}$ does not exist when $K > M$ because $\mathbf{A}^H\mathbf{A}$ becomes singular. Hence, we cannot use (7) to characterize the performance of these sparse linear arrays when $K > M$ . Note that the derivation of (7) assumes a general source covariance matrix $\mathbf{P}$ and does not make use of the prior knowledge that the sources are uncorrelated. If we further assume the sources are uncorrelated in the unconditional model, the unknown parameters become:

\mathbf{\eta}_\mathrm{UU} = [\theta_1, \cdots, \theta_K, p_1, \cdots, p_K, \sigma^2] \in \mathbb{R}^{2K+1},

(8)

where $p_1, \ldots, p_K$ are the diagonal elements of $\mathbf{P}$ . The corresponding CRB is summarized in Theorem 3.

Theorem 3. Under the aforementioned assumptions, plus the prior knowledge of uncorrelated sources, the CRB of the DOAs for the unconditional model is given by^[1]^[13]^[14]^[15]:

\mathrm{CRB}_{\mathbf{\theta}}^\mathrm{UU} = \frac{1}{T} (\mathbf{M}_{\mathbf{\theta}}^H \Pi^\perp_{\mathbf{M}_{\mathbf{p}}} \mathbf{M}_{\mathbf{\theta}})^{-1},

(9)

where

\begin{aligned} \Pi^\perp_{\mathbf{M}_\mathbf{p}} &= \mathbf{I} - \mathbf{M}_\mathbf{p}(\mathbf{M}_\mathbf{p}^H \mathbf{M}_\mathbf{p}) \mathbf{M}_\mathbf{p}^H, \\ \mathbf{M}_\mathbf{\theta} &= (\mathbf{R}^T \otimes \mathbf{R})^{-1/2} \dot{\mathbf{A}}_\mathrm{d} \mathbf{P}, \\ \mathbf{M}_\mathbf{p} &= (\mathbf{R}^T \otimes \mathbf{R})^{-1/2} [\mathbf{A}_\mathrm{d}\ \mathrm{vec}(\mathbf{I}_M)], \\ \mathbf{A}_\mathrm{d} &= \mathbf{A}^* \odot \mathbf{A}, \\ \dot{\mathbf{A}}_\mathrm{d} &= \dot{\mathbf{A}}^* \odot \mathbf{A} + \mathbf{A}^* \odot \dot{\mathbf{A}}. \end{aligned}

Here $\mathrm{vec}(\mathbf{X})$ denotes the vectorization of $\mathbf{X}$ (i.e., stacking the columns of $\mathbf{X}$ into a long vector), $\dot{\mathbf{A}}$ follows the same definition as in Theorem 1, $\otimes$ denotes the Kronecker product, and $\otimes$ denotes the Khatri-Rao product (i.e., the column-wise Kronecker product).

Remark 1. The corresponding Fisher information matrix of (9) is given by^[1]

\mathbf{J} = T\begin{bmatrix} \mathbf{J}_{\mathbf{\theta}\mathbf{\theta}} & \mathbf{J}_{\mathbf{\theta}\mathbf{p}} & \mathbf{J}_{\mathbf{\theta}\mathbf{\sigma^2}} \\ \mathbf{J}_{\mathbf{p}\mathbf{\theta}} & \mathbf{J}_{\mathbf{p}\mathbf{p}} & \mathbf{J}_{\mathbf{p}\mathbf{\sigma^2}} \\ \mathbf{J}_{\mathbf{\sigma^2}\mathbf{\theta}} & \mathbf{J}_{\mathbf{\sigma^2}\mathbf{p}} & \mathbf{J}_{\mathbf{\sigma^2}\mathbf{\sigma^2}} \end{bmatrix},

(10)

where

\begin{aligned} \mathbf{J}_{\mathbf{\theta}\mathbf{\theta}} &= 2\Re[(\dot{\mathbf{A}}^H \mathbf{R}^{-1} \dot{\mathbf{A}})^* \circ (\mathbf{P} \mathbf{A}^H \mathbf{R}^{-1} \mathbf{A} \mathbf{P}) + (\dot{\mathbf{A}}^H \mathbf{R}^{-1} \mathbf{A})^* \circ (\mathbf{P} \mathbf{A}^H \mathbf{R}^{-1} \dot{\mathbf{A}} \mathbf{P})],\\ \mathbf{J}_{\mathbf{p}\mathbf{p}} &= (\mathbf{A}^H \mathbf{R}^{-1} \mathbf{A})^* \circ (\mathbf{A}^H \mathbf{R}^{-1} \mathbf{A}),\\ \mathbf{J}_{\mathbf{\sigma^2}\mathbf{\sigma^2}} &= \mathrm{tr}(\mathbf{R}^{-2}),\\ \mathbf{J}_{\mathbf{\theta}\mathbf{p}} &= 2\Re[(\dot{\mathbf{A}}^H \mathbf{R}^{-1} \mathbf{A})^* \circ (\mathbf{P} \mathbf{A}^H \mathbf{R}^{-1} \mathbf{A})],\\ \mathbf{J}_{\mathbf{\theta}\mathbf{\sigma^2}} &= 2\Re[\mathrm{diag}(\mathbf{P} \dot{\mathbf{A}}^H \mathbf{R}^{-2} \mathbf{A})],\\ \mathbf{J}_{\mathbf{p}\mathbf{\sigma^2}} &= \mathrm{diag}(\mathbf{A}^H \mathbf{R}^{-2} \mathbf{A}), \end{aligned}

(11)

and $\mathbf{J}_{\mathbf{p}\mathbf{\theta}} = \mathbf{J}_{\mathbf{\theta}\mathbf{p}}^H$ , $\mathbf{J}_{\mathbf{\sigma^2}\mathbf{\theta}} = \mathbf{J}_{\mathbf{\theta}\mathbf{\sigma^2}}^H$ , $\mathbf{J}_{\mathbf{\sigma^2}\mathbf{p}} = \mathbf{J}_{\mathbf{p}\mathbf{\sigma^2}}^H$ .

$\mathrm{CRB}_{\mathbf{\theta}}^\mathrm{UU}$ can still be valid even if $K > M$ , and can be used to evaluate the performance of sparse linear arrays in cases when $K > M$ . This is because the invertibility of the Fisher information matrix actually depends of the "difference coarray" steering matrix $\mathbf{A}_\mathrm{d}$ , which may still be full-rank when $K > M$ ^[13]. For certain sparse linear arrays, $\mathrm{CRB}_{\mathbf{\theta}}^\mathrm{UU}$ remains valid for up to $\mathcal{O}(M^2)$ sources^[14]^[15].

Theorem 1-3 summarizes three CRBs in the classical DOA estimation problem. These three CRBs are commonly used in the literature for performance comparison. It should be noted that in the field of array signal processing, there exists many other signal models equipped with difference statistical assumptions. Hence there exists many other CRBs for DOA estimation, and this article only covers a small corner.

Code: the MATLAB code for evaluating the above three CRBs are available on my GitHub. There is also an example that compares the three CRBs numerically.

H. L. Van Trees, Optimum array processing. New York: Wiley, 2002. ↩ ↩ ↩
R. Schmidt, "Multiple emitter location and signal parameter estimation," IEEE Transactions on Antennas and Propagation, vol. 34, no. 3, pp. 276–280, Mar. 1986. ↩
R. Roy and T. Kailath, "ESPRIT-estimation of signal parameters via rotational invariance techniques," IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 37, no. 7, pp. 984–995, Jul. 1989. ↩
Z. Tan, Y. C. Eldar, and A. Nehorai, "Direction of arrival estimation using co-prime arrays: a super resolution viewpoint," IEEE Transactions on Signal Processing, vol. 62, no. 21, pp. 5565–5576, Nov. 2014. ↩
P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and Cramer-Rao bound," IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 37, no. 5, pp. 720–741, May 1989. ↩ ↩
P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and Cramer-Rao bound: further results and comparisons," IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 38, no. 12, pp. 2140–2150, Dec. 1990. ↩
P. Stoica and A. Nehorai, "Performance study of conditional and unconditional direction-of-arrival estimation," IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 38, no. 10, pp. 1783–1795, Oct. 1990. ↩ ↩ ↩
P. Stoica, E. G. Larsson, and A. B. Gershman, "The stochastic CRB for array processing: a textbook derivation," IEEE Signal Processing Letters, vol. 8, no. 5, pp. 148–150, May 2001. ↩
M. Ishiguro, "Minimum redundancy linear arrays for a large number of antennas," Radio Sci., vol. 15, no. 6, pp. 1163–1170, Nov. 1980. ↩
A. Moffet, "Minimum-redundancy linear arrays," IEEE Transactions on Antennas and Propagation, vol. 16, no. 2, pp. 172–175, Mar. 1968. ↩
P. Pal and P. P. Vaidyanathan, "Nested arrays: a novel approach to array processing with enhanced degrees of freedom," IEEE Transactions on Signal Processing, vol. 58, no. 8, pp. 4167–4181, Aug. 2010. ↩
P. Pal and P. P. Vaidyanathan, "Coprime sampling and the MUSIC algorithm," in 2011 IEEE Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop (DSP/SPE), 2011, pp. 289–294. ↩
M. Wang and A. Nehorai, "Coarrays, MUSIC, and the Cramér-Rao bound," IEEE Transactions on Signal Processing, vol. 65, no. 4, pp. 933–946, Feb. 2017. ↩ ↩
C.-L. Liu and P. P. Vaidyanathan, "Cramér–Rao bounds for coprime and other sparse arrays, which find more sources than sensors," Digital Signal Processing, vol. 61, pp. 43–61, Feb. 2017. ↩ ↩
A. Koochakzadeh and P. Pal, "Cramér-Rao bounds for underdetermined source localization," IEEE Signal Processing Letters, vol. 23, no. 7, pp. 919–923, Jul. 2016. ↩ ↩

Mianzhi Wang

Cramér-Rao Bounds in Classical Direction of Arrival Estimation Problems