Mianzhi Wang

Ph.D. in Electrical Engineering

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

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].


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

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

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

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


a(θk)=[exp(j2πλd1sinθk),exp(j2πλd2sinθk),,exp(j2πλdMsinθk)]T,\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 D={0,d0,,(M1)d0}\mathcal{D}=\{0, d_0, \ldots, (M-1)d_0\}, where d0d_0 is usually chosen to be λ/2\lambda/2 to avoid grating lobes.

Given the measurement model (1) and further assumptions on the statistics of s(t)\mathbf{s}(t) and n(t)\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]:

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

We further make the following assumptions:

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

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

For the conditional model,

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

The unknown parameters are

ηC=[θ1,,θK,(x(1)),(x(1)),,(x(T)),(x(T)),σ2]RK+2T+1.\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}.

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

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

where \circ denotes the Hadamard product, and

H=A˙H[IA(AHA)1AH]A˙,A˙=[a(θ1)θ1a(θ2)θ2a(θK)θK],P^=1Tt=1Tx(t)xH(t).\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,

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

The unknown parameters are

ηU=[θ1,,θK,,(Pij),(Pij),,σ2]RK+K2+1.\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}.

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

CRBθU=σ22T{[H(PAHR1AP)T]}1,\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},

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

Note: In the original paper, (7) is obtained from the asymptotic (T0T\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>MK > M) under the unconditional model. However, careful inspection of (7) reveals that CRBθU\mathrm{CRB}^\mathrm{U}_\mathbf{\theta} does not exist when K>MK > M because AHA\mathbf{A}^H\mathbf{A} becomes singular. Hence, we cannot use (7) to characterize the performance of these sparse linear arrays when K>MK > M. Note that the derivation of (7) assumes a general source covariance matrix P\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:

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

where p1,,pKp_1, \ldots, p_K are the diagonal elements of P\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]:

CRBθUU=1T(MθHΠMpMθ)1,\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},


ΠMp=IMp(MpHMp)MpH,Mθ=(RTR)1/2A˙dP,Mp=(RTR)1/2[Ad vec(IM)],Ad=AA,A˙d=A˙A+AA˙.\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 vec(X)\mathrm{vec}(\mathbf{X}) denotes the vectorization of X\mathbf{X} (i.e., stacking the columns of X\mathbf{X} into a long vector), A˙\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]

J=T[JθθJθpJθσ2JpθJppJpσ2Jσ2θJσ2pJσ2σ2],\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},


Jθθ=2[(A˙HR1A˙)(PAHR1AP)+(A˙HR1A)(PAHR1A˙P)],Jpp=(AHR1A)(AHR1A),Jσ2σ2=tr(R2),Jθp=2[(A˙HR1A)(PAHR1A)],Jθσ2=2[diag(PA˙HR2A)],Jpσ2=diag(AHR2A),\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}

and Jpθ=JθpH\mathbf{J}_{\mathbf{p}\mathbf{\theta}} = \mathbf{J}_{\mathbf{\theta}\mathbf{p}}^H, Jσ2θ=Jθσ2H\mathbf{J}_{\mathbf{\sigma^2}\mathbf{\theta}} = \mathbf{J}_{\mathbf{\theta}\mathbf{\sigma^2}}^H, Jσ2p=Jpσ2H\mathbf{J}_{\mathbf{\sigma^2}\mathbf{p}} = \mathbf{J}_{\mathbf{p}\mathbf{\sigma^2}}^H.

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

  1. H. L. Van Trees, Optimum array processing. New York: Wiley, 2002.

  2. R. Schmidt, "Multiple emitter location and signal parameter estimation," IEEE Transactions on Antennas and Propagation, vol. 34, no. 3, pp. 276–280, Mar. 1986.

  3. 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.

  4. 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.

  5. 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.

  6. 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.

  7. 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.

  8. 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.

  9. M. Ishiguro, "Minimum redundancy linear arrays for a large number of antennas," Radio Sci., vol. 15, no. 6, pp. 1163–1170, Nov. 1980.

  10. A. Moffet, "Minimum-redundancy linear arrays," IEEE Transactions on Antennas and Propagation, vol. 16, no. 2, pp. 172–175, Mar. 1968.

  11. 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.

  12. 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.

  13. 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.

  14. 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.

  15. 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.