Mathematical Foundations of DoA Estimation
Signal and Array Model
The mathematical model of the input signal for DoA estimation consists of the array steering vector, the multi-source signal model, and the spatial covariance matrix.
For a uniform linear array with $M$ elements, inter-element spacing $d$, and wavelength $\lambda$, the steering vector is
$$ \mathbf{a}(\theta)= \begin{bmatrix} 1 \\ e^{-j 2\pi \frac{d}{\lambda}\sin\theta} \\ \vdots \\ e^{-j 2\pi (M-1)\frac{d}{\lambda}\sin\theta} \end{bmatrix} $$
For half-wavelength spacing, $d = \frac{\lambda}{2}$
$$ \mathbf{a}(\theta)= \begin{bmatrix} 1 \\ e^{-j\pi\sin\theta} \\ \vdots \\ e^{-j\pi(M-1)\sin\theta} \end{bmatrix} $$
$$ \begin{aligned} \mathbf{a}(\theta) &\in \mathbb{C}^{M} && \text{array steering vector} \\ M &&& \text{number of antenna elements} \\ d &&& \text{inter-element spacing} \\ \lambda &&& \text{signal wavelength} \\ \theta &&& \text{direction of arrival measured from broadside} \\ j^2 &= -1 && \text{imaginary unit} \end{aligned} $$
For $K$ narrowband sources, the received array snapshot is
$$ \mathbf{x}\left[l\right]= \mathbf{A}\mathbf{s}\left[l\right]+ \mathbf{n}\left[l\right] $$
The array manifold matrix is formed from the source steering vectors:
$$ \mathbf{A}= \begin{bmatrix} \mathbf{a}(\theta_1) & \mathbf{a}(\theta_2) & \cdots & \mathbf{a}(\theta_K) \end{bmatrix} $$
$$ \begin{aligned} \mathbf{x}\left[l\right] &\in \mathbb{C}^{M} && \text{received array snapshot} \\ \mathbf{A} &\in \mathbb{C}^{M \times K} && \text{array manifold matrix} \\ \mathbf{s}\left[l\right] &\in \mathbb{C}^{K} && \text{source signal vector} \\ \mathbf{n}\left[l\right] &\in \mathbb{C}^{M} && \text{additive noise vector} \\ K &&& \text{number of sources} \\ l &&& \text{snapshot index} \end{aligned} $$
The steering vector defines the phase response of the array for a given direction.
Covariance Matrix
The theoretical spatial covariance matrix is
$$ \mathbf{R}= \mathbb{E} \left[ \mathbf{x}\left[l\right]\mathbf{x}^{H}\left[l\right] \right] $$
For mutually uncorrelated sources and spatially white noise,
$$ \mathbf{R}= \mathbf{A}\mathbf{R}_s\mathbf{A}^{H}+ \sigma_n^2\mathbf{I}_M $$
For $L$ measured snapshots, the sample covariance matrix is
$$ \widehat{\mathbf{R}}= \frac{1}{L} \sum_{l=0}^{L-1} \mathbf{x}\left[l\right]\mathbf{x}^{H}\left[l\right] $$
$$ \begin{aligned} \mathbf{R} &\in \mathbb{C}^{M \times M} && \text{theoretical covariance matrix} \\ \widehat{\mathbf{R}} &\in \mathbb{C}^{M \times M} && \text{sample covariance matrix} \\ \mathbf{R}_s &\in \mathbb{C}^{K \times K} && \text{source covariance matrix} \\ \sigma_n^2 &&& \text{noise variance} \\ \mathbf{I}_M &\in \mathbb{R}^{M \times M} && \text{identity matrix} \\ L &&& \text{number of snapshots} \\ (\cdot)^H &&& \text{Hermitian transpose} \end{aligned} $$
The covariance matrix represents the spatial power and correlation structure of the received signals.
Spatial FFT Beamformer
The spatial FFT beamformer estimates the angular power distribution by transforming the signals across the antenna elements into the spatial-frequency domain.
For one array snapshot, the beamformer output is
$$ Y[k]= \sum_{m=0}^{M-1} x_m e^{j2\pi m \nu_k} $$
$$ \nu_k = \frac{k-\frac{N_{\mathrm{FFT}}}{2}} {N_{\mathrm{FFT}}} $$
The spatial power spectrum is
$$ P_{\mathrm{FFT}}\left[k \right]= \left|Y[k]\right|^2 $$
For $L$ snapshots, the averaged spectrum is
$$ P_{\mathrm{FFT}}\left[k \right] = \frac{1}{L} \sum_{l=0}^{L-1} \left|Y_l[k]\right|^2 $$
The spatial-frequency bins are mapped to angles as
$$ \theta_k = \arcsin \left( \frac{\lambda}{d}\nu_k \right) $$
For half-wavelength spacing,
$$ \theta_k = \arcsin \left( 2\nu_k \right) $$
$$ \begin{aligned} x_m &&& \text{signal received by antenna element } m \\ M &&& \text{number of antenna elements} \\ N_{\mathrm{FFT}} &&& \text{FFT length} \\ \nu_k &&& \text{normalized spatial frequency} \\ Y[k] &&& \text{spatial FFT output} \\ P_{\mathrm{FFT}}[k] &&& \text{spatial power spectrum} \\ \theta_k &&& \text{angle corresponding to FFT bin } k \end{aligned} $$
The spatial FFT beamformer is computationally efficient and straightforward to implement, but its angular resolution is limited by the array aperture.
Bartlett Beamformer
The Bartlett beamformer estimates the angular power spectrum by steering the array toward each scan angle and measuring the received power.
Bartlett spectrum
$$ P_{\mathrm{Bartlett}}(\theta)= \mathbf{a}^{H}(\theta) \widehat{\mathbf{R}} \mathbf{a}(\theta) $$
Using the snapshot matrix, the same spectrum can be written as
$$ P_{\mathrm{Bartlett}}(\theta) = \frac{1}{L} \sum_{l=0}^{L-1} \left| \mathbf{a}^{H}(\theta) \mathbf{x}\left[ l \right] \right|^2 $$
For a normalized steering vector,
$$ \widetilde{\mathbf{a}}(\theta)= \frac{\mathbf{a}(\theta)} {\sqrt{\mathbf{a}^{H}(\theta)\mathbf{a}(\theta)}} $$
$$ P_{\mathrm{Bartlett}}(\theta)= \widetilde{\mathbf{a}}^{H}(\theta) \widehat{\mathbf{R}} \widetilde{\mathbf{a}}(\theta) $$
$$ \begin{aligned} P_{\mathrm{Bartlett}}(\theta) &&& \text{Bartlett spatial power spectrum} \\ \mathbf{a}(\theta) &\in \mathbb{C}^{M} && \text{steering vector} \\ \widehat{\mathbf{R}} &\in \mathbb{C}^{M \times M} && \text{sample covariance matrix} \\ \mathbf{x}\left[l \right] &\in \mathbb{C}^{M} && \text{received array snapshot} \\ L &&& \text{number of snapshots} \\ (\cdot)^H &&& \text{Hermitian transpose} \end{aligned} $$
The Bartlett beamformer is simple, robust, and applicable to arbitrary array geometries, but its angular resolution remains limited by the physical array aperture.
Capon MVDR Beamformer
The Capon MVDR (Minimum Variance Distortionless Response) beamformer estimates the angular power spectrum by minimizing the array output power while preserving unity gain in the scan direction.
MVDR optimization
The resulting beamforming weights are
$$ \mathbf{w}_{\mathrm{MVDR}}(\theta)= \frac{ \widehat{\mathbf{R}}^{-1} \mathbf{a}(\theta) }{ \mathbf{a}^{H}(\theta) \widehat{\mathbf{R}}^{-1} \mathbf{a}(\theta) } $$
Capon spectrum
$$ P_{\mathrm{Capon}}(\theta)= \frac{ 1 }{ \mathbf{a}^{H}(\theta) \widehat{\mathbf{R}}^{-1} \mathbf{a}(\theta) } $$
For improved numerical stability, diagonal loading can be applied:
$$ \widehat{\mathbf{R}}_{\mathrm{DL}}= \widehat{\mathbf{R}} + \delta \mathbf{I}_M $$
$$ P_{\mathrm{Capon}}(\theta)= \frac{ 1 }{ \mathbf{a}^{H}(\theta) \widehat{\mathbf{R}}_{\mathrm{DL}}^{-1} \mathbf{a}(\theta) } $$
$$ \begin{aligned} P_{\mathrm{Capon}}(\theta) &&& \text{Capon spatial power spectrum} \\ \mathbf{w}_{\mathrm{MVDR}}(\theta) &\in \mathbb{C}^{M} && \text{adaptive beamforming weights} \\ \mathbf{a}(\theta) &\in \mathbb{C}^{M} && \text{steering vector} \\ \widehat{\mathbf{R}} &\in \mathbb{C}^{M \times M} && \text{sample covariance matrix} \\ \delta &&& \text{diagonal loading factor} \\ \mathbf{I}_M &\in \mathbb{R}^{M \times M} && \text{identity matrix} \\ \left(\cdot\right)^H &&& \text{Hermitian transpose} \end{aligned} $$
The Capon MVDR beamformer provides higher angular resolution and interference suppression than conventional beamforming, but it requires covariance-matrix inversion and is sensitive to estimation errors and array mismatch.
$$ \begin{aligned} \underset{\mathbf{w}}{\operatorname{minimize}} \quad & \mathbf{w}^{H} \widehat{\mathbf{R}} \mathbf{w} \\ \operatorname{subject\ to} \quad & \mathbf{w}^{H} \mathbf{a}(\theta)= 1 \end{aligned} $$
MUSIC
MUSIC (MUltiple SIgnal Classification) estimates the directions of arrival by exploiting the orthogonality between the steering vectors of the sources and the noise subspace of the covariance matrix.
Eigendecomposition
$$ \widehat{\mathbf{R}}= \mathbf{E} \boldsymbol{\Lambda} \mathbf{E}^{H} $$
$$ \mathbf{E}= \begin{bmatrix} \mathbf{e}_1 & \mathbf{e}_2 & \cdots & \mathbf{e}_M \end{bmatrix} $$
$$ \boldsymbol{\Lambda}= \operatorname{diag} \begin{bmatrix} \lambda_1 & \lambda_2 & \cdots & \lambda_M \end{bmatrix} $$
For $K$ sources, the eigenvectors associated with the $M−K$ smallest eigenvalues form the noise subspace:
$$ \mathbf{E}_n = \begin{bmatrix} \mathbf{e}(K+1) & \mathbf{e}(K+2) & \cdots & \mathbf{e}(M) \end{bmatrix} $$
Orthogonality condition
$$ \mathbf{E}_n^{H} \mathbf{a}(\theta_k)= \mathbf{0} $$
MUSIC pseudospectrum
$$ P_{\mathrm{MUSIC}}(\theta)= \frac{ 1 }{ \mathbf{a}^{H}(\theta) \mathbf{E}_n \mathbf{E}_n^{H} \mathbf{a}(\theta) } $$
$$ \begin{aligned} P_{\mathrm{MUSIC}}(\theta) & && \text{MUSIC pseudospectrum} \\ \mathbf{E} &\in \mathbb{C}^{M \times M} && \text{eigenvector matrix} \\ \boldsymbol{\Lambda} &\in \mathbb{R}^{M \times M} && \text{diagonal eigenvalue matrix} \\ \mathbf{E}_n &\in \mathbb{C}^{M \times (M-K)} && \text{noise-subspace matrix} \\ \lambda_m & && \text{eigenvalue} \\ K & && \text{number of sources} \end{aligned} $$
MUSIC provides very high angular resolution, but it requires the number of sources to be known and is sensitive to covariance errors, low SNR, coherent sources, and array mismatch.
ESPRIT
ESPRIT (Estimation of Signal Parameters via Rotational Invariance Techniques) estimates the directions of arrival from the rotational relationship between two overlapping, identically spaced subarrays.
Signal subspace
For $K$ sources, the signal subspace is formed by the eigenvectors associated with the $K$ largest eigenvalues:
$$ \mathbf{E}_s= \begin{bmatrix} \mathbf{e}_1 & \mathbf{e}_2 & \cdots & \mathbf{e}_K \end{bmatrix} $$
Shift-invariant subarrays
$$ \mathbf{E}_1 = \mathbf{J}_1 \mathbf{E}_s $$
$$ \mathbf{E}_2 = \mathbf{J}_2 \mathbf{E}_s $$
The two subspaces are related by
$$ \mathbf{E}_2 = \mathbf{E}_1 \boldsymbol{\Psi} $$
The rotation matrix is estimated as
$$ \boldsymbol{\Psi}= \mathbf{E}_1^{\dagger} \mathbf{E}_2 $$
Eigenvalues and DoA estimates
$$ \boldsymbol{\Psi}\mathbf{v}_k = \mu_k\mathbf{v}_k $$
For the steering-vector convention used above,
$$ \mu_k = e^{-j2\pi\frac{d}{\lambda}\sin\theta_k} $$
$$ \theta_k=\arcsin\left(-\frac{\lambda}{2\pi d}\arg(\mu_k)\right) $$
For half-wavelength spacing,
$$ \theta_k=\arcsin\left(-\frac{1}{\pi}\arg(\mu_k)\right) $$
$$ \begin{aligned} \mathbf{E}_s &\in \mathbb{C}^{M \times K} && \text{signal-subspace matrix} \\ \mathbf{J}_1,\mathbf{J}_2 &&& \text{subarray selection matrices} \\ \boldsymbol{\Psi} &\in \mathbb{C}^{K \times K} && \text{rotation matrix} \\ \mu_k &&& \text{eigenvalue of the rotation matrix} \\ \left(\cdot\right)^{\dagger} &&& \text{Moore--Penrose pseudoinverse} \end{aligned} $$
ESPRIT provides direct high-resolution angle estimates without an angular scan, but it requires the number of sources to be known and an array with shift-invariant geometry.
References
Van Trees, H. L., Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory. New York: John Wiley & Sons, 2002. DOI: 10.1002/0471221104.
J. Li and P. Stoica, eds., MIMO Radar Signal Processing. Hoboken, NJ: John Wiley & Sons, 2009. DOI: 10.1002/9780470391488.