Yule–Walker Estimation (Kay / Proakis Convention)

Theory

As per standard signal processing literature (e.g., Kay, Proakis), many signals can be modeled as autoregressive (AR) process i.e. a linear combination of past samples plus white noise. The spectrum of these AR modelled signal can be estimated using a parametric spectral estimation approach, where the power spectral density (PSD) can be directly evaluated using the AR model parameters which describes the signal. This is a parametric spectral estimation approach, where the power spectral density (PSD) is derived from the model parameters.

We assume that $x[n]$ is a zero-mean, wide-sense stationary (WSS) random process generated by passing white noise $e[n]$ through an all-pole system.

An all-pole system is a system whose transfer function has only poles (no finite zeros), meaning its behavior is entirely determined by the denominator polynomial. In an AR model, white noise is passed through this all-pole filter, and the AR coefficients define the pole locations. These poles shape the frequency response and thus determine the power spectral density of the signal.

Autoregressive (AR) Model

The AR model of order $p$ is written as:

$$ x[n] = \sum_{k=1}^{p} \phi_k x[n-k] + e[n] $$

• $\phi_k$ = AR coefficients
• $e[n]$ = white noise with $E[e[n]] = 0$ and variance $\sigma_e^2$
• $p$ = model order

Autocorrelation Function

The autocorrelation function is defined as:

$$ R[k] = E[x[n] x^*[n-k]] $$

For real-valued signals:

$$ R[k] = E[x[n] x[n-k]] $$

For a WSS process, $R[k]$ depends only on the lag $k$.

Yule–Walker Equations

To derive the Yule–Walker equations, we multiply the AR model by $x[n-m]$ and take expectations:

$$ E[x[n]x[n-m]] = \sum_{k=1}^{p} \phi_k E[x[n-k]x[n-m]] + E[e[n]x[n-m]] $$

Since $e[n]$ is white noise and uncorrelated with past samples:

$$ E[e[n]x[n-m]] = 0 \quad \text{for } m \ge 1 $$

Using the definition of autocorrelation, we obtain:

$$ R[m] = \sum_{k=1}^{p} \phi_k R[m-k], \quad m = 1,2,\dots,p $$

These equations can be written in matrix form as:

$$ \begin{bmatrix} R[0] & R[1] & \cdots & R[p-1] \\ R[1] & R[0] & \cdots & R[p-2] \\ \vdots & \vdots & \ddots & \vdots \\ R[p-1] & R[p-2] & \cdots & R[0] \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \vdots \\ \phi_p \end{bmatrix} = \begin{bmatrix} R[1] \\ R[2] \\ \vdots \\ R[p] \end{bmatrix} $$

The autocorrelation matrix is Toeplitz because its entries depend only on the difference between the row and column indices (i.e., the lag). This structured form significantly reduces the computational complexity of solving the Yule–Walker equations through efficient recursive methods such as the Levinson–Durbin algorithm.

Noise Variance

Once the AR coefficients are estimated, the variance prediction error is given by:

$$ \sigma_e^2 = R[0] - \sum_{k=1}^{p} \phi_k R[k] $$

This represents the portion of the signal power that cannot be explained by the AR model. In other words, even after using past samples to estimate the current value, there is still some remaining randomness or error. That leftover uncertainty measures error variance.

Power Spectral Density (PSD)

The PSD of the AR process is:

$$ P_x(f) = \frac{\sigma_e^2}{\left|1 - \sum_{k=1}^{p} \phi_k e^{-j 2 \pi f k}\right|^2} $$

The denominator corresponds to the frequency response of the all-pole system. The locations of the poles determine the spectral characteristics of the process. Frequencies at which the denominator has a small magnitude correspond to resonances of the filter, resulting in peaks in the estimated power spectrum. In summary, the AR coefficients control the shape and location of the spectral peaks.

$$ -\tfrac{1}{2} \le f \le \tfrac{1}{2}, \quad \omega = 2\pi f $$

In discrete-time signal processing, frequency is periodic, so values outside a fundamental interval represent repeated information. This is similar to angular frequency in continuous-time systems, which is periodic with period $2\pi$.

Due to this periodicity, normalized frequency $f$ is typically restricted to the principal interval $-\tfrac{1}{2} \le f \le \tfrac{1}{2}$, which captures all unique spectral content.

The relationship $\omega = 2\pi f$ converts normalized frequency into angular frequency in radians/sample, scaling the interval so that one full cycle corresponds to $2\pi$ radians.

Advantages

• Provides high spectral resolution, particularly for short data records.
• Produces smooth PSD estimates without the leakage commonly observed in periodograms.
• Requires estimation of only a small number of model parameters.
• Exploits the Toeplitz structure of the autocorrelation matrix for computationally efficient implementation.
• Is especially effective for signals that can be accurately represented by an autoregressive model.

In summary, Yule–Walker estimation determines the AR parameters from autocorrelation values and then computes the spectrum analytically from the resulting model.