Parametric Spectral Estimation | DSP-101

Beyond Non-Parametric Methods

The periodogram, Welch's method, and the spectrogram are all non-parametric — they make no assumptions about the underlying signal model. They estimate the PSD directly from the data using the DFT. Their great strength is generality: they work for any signal. Their fundamental limitation is resolution: the frequency resolution is bounded below by 1/(NT), where N is the data length and T is the sample interval. You cannot resolve two sinusoids closer than this limit, no matter how much you process the data.

Parametric methods take a different approach. Instead of computing the spectrum directly from data, they first fit a signal model to the data — typically an autoregressive (AR), moving-average (MA), or ARMA model — and then compute the PSD analytically from the fitted model parameters. Because the model imposes structure on the signal, parametric methods can achieve super-resolution: resolving spectral features that lie closer than 1/(NT).

When to Use Parametric Methods

Parametric methods shine when the signal is well-described by the assumed model. They excel for short data records, closely spaced sinusoids, and applications where the number and approximate location of spectral peaks are known in advance. They can fail dramatically when the assumed model does not match the true signal structure.

Parametric = super-resolution but model-dependent. Non-parametric = lower resolution but assumption-free.

The Autoregressive (AR) Model

The most widely used parametric model is the autoregressive (AR) model of order p, also called the all-pole model. In an AR(p) model, each sample of the signal is a linear combination of the p previous samples plus a white noise input:

AR(p) Signal Model

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

a₁, a₂, …, aₚ are the AR coefficients to be estimated. w[n] ~ WN(0, σ²) is white noise with variance σ². The AR model is equivalent to an all-pole IIR filter driven by white noise.

The AR PSD is obtained by computing the power spectrum of the all-pole filter response. If A(ω) = 1 + a₁e^{-jω} + a₂e^{-j2ω} + … + aₚe^{-jpω}, then:

AR Power Spectral Density

P_{AR}(\omega) = \frac{\sigma^2}{\left|A(\omega)\right|^2} = \frac{\sigma^2}{\left|1 + \sum_{k=1}^{p} a_k\, e^{-j\omega k}\right|^2}

The AR PSD is a smooth, rational function of frequency with peaks where |A(ω)|² is small — i.e., near the poles of the all-pole filter. The resolution is not limited by data length — it is limited only by model order and SNR.

Estimating AR Coefficients: The Yule-Walker Method

To use the AR PSD, we must first estimate the coefficients a₁, …, aₚ and the noise variance σ² from the observed data. The most classical approach is the Yule-Walker method, which exploits the relationship between AR coefficients and the autocorrelation sequence of the signal.

The Yule-Walker equations link the AR parameters directly to the autocorrelation values r[0], r[1], …, r[p]:

Yule-Walker Equations (matrix form)

\mathbf{R}\,\mathbf{a} = -\mathbf{r}, \quad \mathbf{R} = \begin{bmatrix} r[0] & r[1] & \cdots & r[p-1] \\ r[1] & r[0] & \cdots & r[p-2] \\ \vdots & & \ddots & \vdots \\ r[p-1] & r[p-2] & \cdots & r[0] \end{bmatrix}

𝐑 is the p×p Toeplitz autocorrelation matrix with R_{ij} = r[|i−j|]. 𝐚 = [a₁, …, aₚ]ᵀ. 𝐫 = [r[1], …, r[p]]ᵀ. Solve by the Levinson-Durbin algorithm in O(p²) operations. The noise variance is σ² = r[0] + 𝐚ᵀ𝐫.

In practice, the true autocorrelation is unknown and must be estimated from the data. Substituting the sample autocorrelation r̂[k] into the Yule-Walker equations gives the Yule-Walker AR estimate. An alternative, often preferred for short records, is the Burg method, which minimizes forward and backward prediction errors simultaneously and guarantees a stable AR filter.

Model Order Selection

The AR model order p is a critical parameter — it controls the trade-off between resolution and variance. Too small and the PSD is oversmoothed; too large and spurious peaks appear. Three principled selection criteria are widely used:

Criterion	Formula	Bias
AIC (Akaike Information Criterion)	AIC(p) = N ln σ̂²(p) + 2p	Tends to overestimate (asymptotically inconsistent)
MDL (Minimum Description Length)	MDL(p) = N ln σ̂²(p) + p ln N	Consistent — converges to true order as N→∞
FPE (Final Prediction Error)	FPE(p) = σ̂²(p) · (N+p+1)/(N−p−1)	Similar to AIC, minimized at optimal p

In practice, compute the criterion for p = 1, 2, …, p_max and choose the p that minimizes it. For noisy data, the MDL is preferred. When prior knowledge about the number of sinusoidal components is available, set p ≈ 2 × (number of sinusoids) as a starting point.

Subspace Methods: MUSIC and ESPRIT

The most powerful parametric techniques for estimating frequencies of sinusoids in noise are the subspace methods. These exploit the eigenstructure of the data covariance matrix, decomposing it into a signal subspace and a noise subspace.

Consider a signal consisting of M complex sinusoids in white noise: x[n] = ∑ Aₖe^{jωₖn} + w[n]. The covariance matrix 𝐑_x has M large eigenvalues (signal subspace) and N−M equal eigenvalues equal to the noise variance σ² (noise subspace). The key insight: the steering vectors at the true frequencies e(ωₖ) = [1, e^{jωₖ}, …, e^{j(N-1)ωₖ}]ᵀ are orthogonal to the noise subspace.

The MUSIC (MUltiple SIgnal Classification) algorithm exploits this orthogonality to form a pseudospectrum:

MUSIC Pseudospectrum

P_{MUSIC}(\omega) = \frac{1}{\mathbf{e}^H(\omega)\,\mathbf{E}_N\,\mathbf{E}_N^H\,\mathbf{e}(\omega)}

𝐄_N is the matrix of noise eigenvectors (columns). The pseudospectrum P_MUSIC(ω) has sharp peaks at the true frequencies because 𝐞(ω)ᴴ𝐄_N𝐄_Nᴴ𝐞(ω) → 0 when ω = ωₖ. The MUSIC peaks are not true power estimates — they locate frequencies only.

ESPRIT (Estimation of Signal Parameters via Rotational Invariance Techniques) is a related method that extracts frequency estimates directly from the signal subspace without computing a pseudospectrum. It exploits a shift invariance property of the array manifold and solves a small eigenvalue problem, making it computationally more efficient than MUSIC for high-resolution frequency estimation.

MUSIC vs. ESPRIT

MUSIC produces a pseudospectrum (evaluated on a frequency grid) and requires a peak search. ESPRIT gives direct frequency estimates by solving an eigenvalue problem — no spectral search needed. Both require knowledge of the number of signal components M; both achieve super-resolution at high SNR. ESPRIT is preferred for computational efficiency; MUSIC for visual interpretation of the pseudospectrum.

Both methods require M (number of sinusoids) as prior input — misspecifying M degrades performance severely.

Parametric vs. Non-Parametric: When to Use Each

Parametric Methods

Super-resolution — resolves closely spaced spectral lines
Performs well with short data records (N small)
Smooth, analytically derived PSD (AR)
Best for line spectra: sinusoids in noise
Fails if model order or signal structure is wrong
Sensitive to model mismatch and noise coloring

Non-Parametric Methods

No model assumptions — works for any signal type
Resolution bounded by 1/(NT) — requires long records
Welch's method: stable, low-variance estimates
Best for broadband, continuous, or unknown spectra
Immune to model mismatch
Standard first choice in exploratory analysis

Applications of Parametric Spectral Estimation

Domain	Application	Method Typically Used
Radar / Sonar	DOA (direction-of-arrival) estimation from antenna arrays	MUSIC, ESPRIT
Speech	Linear predictive coding (LPC) — speech compression and synthesis	AR (Burg or Levinson-Durbin)
Biomedical	Heart rate variability analysis, EEG spectral peaks	AR (Yule-Walker)
Seismology	Resolving closely spaced normal modes of Earth	MUSIC, Max-Entropy
Communications	Frequency estimation of narrowband interferers	ESPRIT, MUSIC
NMR / MRI	Chemical shift estimation from short free induction decay signals	MUSIC, LP (linear prediction)

Practical Considerations

Parametric methods have several practical pitfalls. First, they all require the model order as input — either the AR order p or the number of sinusoids M. Overestimating M in MUSIC creates spurious peaks; underestimating misses real ones. Second, the covariance matrix estimate becomes unreliable when the number of snapshots (data blocks) is small relative to the array size. Third, coherent signals (two signals with a fixed phase relationship, as in multipath propagation) cause the covariance matrix to lose rank, degrading MUSIC performance. Spatial smoothing pre-processing restores full rank in these cases.

For real-valued signals, efficient implementations use the real autocorrelation structure. The Levinson-Durbin recursion solves the Yule-Walker equations in O(p²) operations for any AR order, making AR estimation computationally cheap even at high orders.

Key Takeaways

Parametric methods fit a signal model (AR, ARMA, or subspace) to data before estimating the PSD — enabling super-resolution beyond the 1/(NT) DFT limit.
The AR(p) model is an all-pole filter driven by white noise; its PSD is σ²/|A(ω)|² — smooth peaks at pole locations.
The Yule-Walker equations relate AR coefficients to the autocorrelation sequence and are solved by the O(p²) Levinson-Durbin recursion.
Model order p must be chosen carefully — AIC, MDL, and FPE criteria provide principled guidance.
MUSIC and ESPRIT are subspace methods for resolving sinusoids: they exploit the orthogonality of signal and noise eigenvectors of the covariance matrix.
Parametric methods work best with short records, high SNR, and line spectra; non-parametric methods are preferred for broadband or unknown spectra.
Coherent signals (multipath) degrade subspace methods — spatial smoothing restores rank.