Parametric Spectrum Estimation: Yule-Walker, Burg, MUSIC, ESPRIT, and MVDR

Why Parametric?

Nonparametric methods are limited by the Rayleigh resolution limit: two tones closer than $F_s/N$ cannot be separated, no matter how much averaging you do. A parametric method assumes the signal is generated by a model with a few parameters, estimates those parameters, and computes the spectrum from the model. The result can resolve far below the Rayleigh limit.

The trade-off: if the model is wrong, the estimate is wrong — and confidently so.

---

AR Spectral Estimation

Assume the signal is an AR process of order $p$:

\[x[n]=-\sum_{k=1}^{p}a_kx[n-k]+w[n]\]

The spectrum is

\[\hat P_{AR}(e^{j\omega})=\frac{\hat\sigma_w^2}{\left\vert1+\sum_{k=1}^{p}\hat a_k e^{-j\omega k}\right\vert^2}\]

The denominator creates sharp peaks when the poles approach the unit circle. This is why AR gives high resolution from short records.

Choosing the Order

Criterion	Formula	Behaviour
AIC	$N\ln\hat\sigma_p^2+2p$	Tends to overestimate
MDL (BIC)	$N\ln\hat\sigma_p^2+p\ln N$	Consistent, penalizes more
CAT	$\frac{1}{N}\sum\dots$	Similar to AIC

For a single sinusoid in noise, order 2 per tone plus a few for the noise works well. For speech LPC, $p\approx F_s/1000+2$.

---

The Yule–Walker Method

Estimate the autocorrelation from the data, then solve the normal equations:

\[\hat r[l]=\frac{1}{N}\sum_{n=0}^{N-1-l}x[n]x[n+l],\qquad l=0,\dots,p\] \[\begin{bmatrix}\hat r[0]&\hat r[1]&\dots&\hat r[p-1]\\\hat r[1]&\hat r[0]&\dots&\hat r[p-2]\\\vdots&&\ddots&\vdots\\\hat r[p-1]&\hat r[p-2]&\dots&\hat r[0]\end{bmatrix}\begin{bmatrix}\hat a_1\\\vdots\\\hat a_p\end{bmatrix}=-\begin{bmatrix}\hat r[1]\\\vdots\\\hat r[p]\end{bmatrix}\]

Solve with Levinson-Durbin.

Properties

Advantages:

• Guaranteed stable all-pole model ($\vert k_m\vert<1$ for a valid autocorrelation).
• Computationally efficient ($O(p^2)$).
• Simple to implement.

Disadvantages:

• The biased autocorrelation estimate (dividing by $N$ regardless of lag) tapers the data, which biases and broadens the spectral peaks.
• The unbiased estimate can produce an invalid autocorrelation matrix that is not positive definite, and then Levinson-Durbin fails.
• Poor resolution for short records with high SNR, compared with Burg.

---

The Burg Method

Burg estimates the reflection coefficients directly by minimizing the forward and backward prediction errors at each stage, constrained to $\vert k_m\vert<1$.

At stage $m$:

\[\hat k_m=\frac{-2\sum_n e_f^{(m-1)}[n]\,e_b^{(m-1)}[n-1]}{\sum_n\left(\left(e_f^{(m-1)}[n]\right)^2+\left(e_b^{(m-1)}[n-1]\right)^2\right)}\]

Then update the errors with the lattice recursion.

Properties

Advantages:

• Higher resolution than Yule–Walker for short records.
• Stable by construction: $\vert k_m\vert<1$ is enforced.
• No autocorrelation estimate, so no windowing bias.
• Uses both forward and backward errors — effectively twice the data.

Disadvantages:

• Spectral peak splitting at high SNR with long records (the estimated pole splits into two).
• Spectral line shifting — the peak position drifts slightly with SNR and record length.
• Slightly more computation than Yule–Walker.

Yule–Walker vs Burg

Aspect	Yule–Walker	Burg
Autocorrelation needed	Yes	No
Stability	Guaranteed for valid $r$	Guaranteed always
Resolution	Lower	Higher
Peak splitting	No	Yes, at high SNR
Line shifting	Small	Present
Computation	$O(p^2)$	$O(pN)$ per stage

---

MA and ARMA Estimation

MA(q)

\[x[n]=\sum_{k=0}^{q}b_kw[n-k]\] \[P_{MA}(e^{j\omega})=\sigma_w^2\left\vert\sum_{k=0}^{q}b_ke^{-j\omega k}\right\vert^2\]

The autocorrelation is nonzero only for lags $0\le l\le q$, so the parameters come from the modified Yule–Walker equations. MA needs a high order to represent sharp peaks.

ARMA(p,q)

\[P_{ARMA}(e^{j\omega})=\sigma_w^2\frac{\left\vert\sum_{k=0}^{q}b_ke^{-j\omega k}\right\vert^2}{\left\vert1+\sum_{k=1}^{p}a_ke^{-j\omega k}\right\vert^2}\]

Estimation is nonlinear. Standard approaches:

1. Modified Yule–Walker: estimate the AR part from the autocorrelation at lags $q+1,\dots,q+p$, then estimate the MA part from the residual.
2. Durbin’s method: fit a high-order AR, then factor into AR and MA.
3. Nonlinear least squares / iterative: best accuracy, most computation.

ARMA fits peaks and notches with fewer parameters than pure AR or MA, but the estimation is harder.

---

Pisarenko Harmonic Decomposition

The first subspace method. Assume $M$ exponentials in white noise:

\[x[n]=\sum_{i=1}^{M}A_ie^{j(\omega_in+\phi_i)}+w[n]\]

The autocorrelation matrix $\mathbf{R}_{p}$ of size $(p+1)\times(p+1)$ with $p\ge M$ has eigenvalues:

• $M$ large eigenvalues (signal subspace),
• $p+1-M$ small eigenvalues equal to $\sigma_w^2$ (noise subspace).

The smallest eigenvector $\mathbf{v}$ is orthogonal to the signal subspace, and its polynomial

\[V(z)=\sum_{k=0}^{p}v_kz^{-k}\]

has roots on the unit circle at the signal frequencies.

Properties

• Resolution far beyond the Rayleigh limit.
• Needs only one eigenvector (the smallest).
• Uses only $M$ of the $p$ noise eigenvalues → statistically inefficient.
• Sensitive to noise and to an incorrect $M$.

---

MUSIC

MUltiple SIgnal Classification uses the whole noise subspace.

The Step

1. Estimate the $p\times p$ autocorrelation matrix $\hat{\mathbf{R}}$ from the data.
2. Eigendecompose: $\hat{\mathbf{R}}=\mathbf{U}_s\boldsymbol{\Lambda}_s\mathbf{U}_s^H+\mathbf{U}_n\boldsymbol{\Lambda}_n\mathbf{U}_n^H$.
3. The noise subspace $\mathbf{U}_n$ spans the $(p-M)$ smallest eigenvectors.
4. The steering vector at frequency $\omega$ is

\[\mathbf{a}(\omega)=\left[1,\ e^{j\omega},\ e^{2j\omega},\dots,e^{j(p-1)\omega}\right]^T\]

1. The pseudospectrum is

\[\hat P_{MUSIC}(e^{j\omega})=\frac{1}{\mathbf{a}^H(\omega)\mathbf{U}_n\mathbf{U}_n^H\mathbf{a}(\omega)}\]

At a true signal frequency, $\mathbf{a}(\omega)$ is orthogonal to the noise subspace, so the denominator is zero and the pseudospectrum peaks sharply.

Properties

Aspect	Value
Resolution	Far below Rayleigh, approaches the Cramér–Rao bound
Consistency	Yes, as $N\to\infty$
Requirement	Correct model order $M$
Cost	Eigendecomposition $O(p^3)$
Weakness	Coherent sources (multipath) decorrelate the matrix and break it

Coherent Sources

If two sources are perfectly correlated (multipath), the signal covariance matrix loses rank and MUSIC fails. The fix is spatial smoothing: average subarray covariance matrices to restore rank.

---

ESPRIT

Estimation of Signal Parameters via Rotational Invariance Techniques avoids the search over $\omega$.

The Rotational Invariance

Split the steering matrix into two overlapping subarrays:

\[\mathbf{U}_s=\begin{bmatrix}\mathbf{U}_1\\ \text{last row}\end{bmatrix}=\begin{bmatrix}\text{first row}\\ \mathbf{U}_2\end{bmatrix}\]

The two subarrays differ by a rotation $\boldsymbol{\Phi}=\operatorname{diag}(e^{j\omega_1},\dots,e^{j\omega_M})$:

\[\mathbf{U}_2=\mathbf{U}_1\boldsymbol{\Phi}\]

Solve in the least-squares sense:

\[\boldsymbol{\Psi}=\left(\mathbf{U}_1^H\mathbf{U}_1\right)^{-1}\mathbf{U}_1^H\mathbf{U}_2\]

The eigenvalues of $\boldsymbol{\Psi}$ are $e^{j\omega_i}$, so the frequencies are read directly from the eigenvalue phases — no spectral search.

MUSIC vs ESPRIT

Aspect	MUSIC	ESPRIT
Search	1-D spectral search over $\omega$	None — closed form
Cost	Eigendecomposition + search	Eigendecomposition + small eigenproblem
Accuracy	Slightly better	Slightly worse but close
Array requirement	Any	Requires shift-invariant structure
Coherent sources	Fails without smoothing	Also fails without smoothing

In practice ESPRIT is preferred when the array is shift-invariant (ULA), because it is much faster and the accuracy difference is small.

---

MVDR (Capon) Beamformer

The minimum variance distortionless response estimator minimizes output power subject to unity gain at the look frequency:

\[\min_{\mathbf{w}}\mathbf{w}^H\hat{\mathbf{R}}\mathbf{w}\quad\text{subject to}\quad\mathbf{w}^H\mathbf{a}(\omega)=1\]

The solution is

\[\mathbf{w}=\frac{\hat{\mathbf{R}}^{-1}\mathbf{a}(\omega)}{\mathbf{a}^H(\omega)\hat{\mathbf{R}}^{-1}\mathbf{a}(\omega)}\]

and the spectrum is

\[\hat P_{MVDR}(e^{j\omega})=\frac{1}{\mathbf{a}^H(\omega)\hat{\mathbf{R}}^{-1}\mathbf{a}(\omega)}\]

Properties

• Better resolution than AR, worse than MUSIC.
• Does not need the model order — no eigendecomposition.
• Robust to coherent sources, unlike MUSIC.
• Requires inverting $\hat{\mathbf{R}}$, which needs $N>p$ and can be ill-conditioned.

Method	Resolution	Needs order?	Coherent sources	Cost
Periodogram	Rayleigh	No	OK	Low
AR (Yule–Walker)	Good	Yes	OK	$O(p^2)$
AR (Burg)	Better	Yes	OK	$O(pN)$
Pisarenko	High	Yes	Fails	$O(p^3)$
MUSIC	Highest	Yes	Fails	$O(p^3)$ + search
ESPRIT	Highest	Yes	Fails	$O(p^3)$
MVDR	High	No	OK	$O(p^3)$

---

Model Order Selection

All subspace methods need $M$, the number of sources. Getting it wrong is the main failure mode.

Information Criteria

\(\text{AIC}(M)=-2\ln L(\hat\theta_M)+2k_M\) \(\text{MDL}(M)=-\ln L(\hat\theta_M)+\tfrac12k_M\ln N\)

For the eigenvalue problem, the practical form uses the noise eigenvalues:

\[\text{MDL}(M)=-\ln\left(\frac{\left(\prod_{i=M+1}^{p}\lambda_i\right)^{1/(p-M)}}{\frac{1}{p-M}\sum_{i=M+1}^{p}\lambda_i}\right)^{(p-M)N}+\tfrac12M(2p-M)\ln N\]

Choose the $M$ that minimizes the criterion.

Practical Rules

• Eigenvalue gap: plot the eigenvalues and look for the knee where they drop to the noise floor.
• Overestimate slightly rather than underestimate — an extra source adds a spurious peak but does not lose a real one.
• Coherent sources: the gap disappears; use spatial smoothing or a different method.
• Low SNR / short record: increase $p$ (the matrix size), not just $M$.

---

Practical Notes

Problem	Cause	Fix
Spurious peaks	Order too high	Use MDL, not AIC
Missed weak source	Order too low	Increase $M$
MUSIC fails on multipath	Coherent sources	Spatial smoothing
Burg peak splitting	High SNR, long record	Yule–Walker or reduce order
Unstable AR filter	Invalid autocorrelation	Biased estimator or Burg
MVDR matrix singular	$N\le p$	Increase $N$ or diagonal loading

Diagonal loading: add $\epsilon\mathbf{I}$ to $\hat{\mathbf{R}}$ before inverting. It regularizes the solution and makes MVDR robust to steering-vector errors.

---

Summary

1. Parametric methods beat the Rayleigh limit by fitting a model rather than smoothing data.
2. Yule–Walker is simple and stable but biased; Burg has higher resolution and is always stable but can split peaks.
3. Pisarenko uses the smallest noise eigenvector; MUSIC uses the whole noise subspace and gives a sharp pseudospectrum.
4. ESPRIT reads frequencies from the eigenvalues of a rotation operator — no search needed.
5. MVDR (Capon) needs no model order and handles coherent sources.
6. Model order selection with MDL is the make-or-break step for all subspace methods.