Power Spectral Density Estimation

Welch's method [83] (or the periodogram method [18]) for estimating power spectral densities (PSD) is carried out by dividing the time signal into successive blocks, and averaging squared-magnitude DFTs of the signal blocks. Let $x_m(n)=x(n+mN)$, $n=0,1,\dots,N-1$, denote the $m$th block of the signal $x\in{\bf C}^{MN}$, with $M$ denoting the number of blocks. Then the Welch PSD estimate is given by

$${\hat R}_x(\omega_k) = \frac{1}{M}\sum_{m=0}^{M-1}\left\vert DFT_... ...t\vert^2 \isdef \left\{\left\vert X_m(\omega_k)^2\right\vert\right\}_m \protect$$ (8.3)

where `` $\{\cdot\}_m$'' denotes time averaging across blocks (or ``frames'') of data indexed by $m$. The function pwelch implements Welch's method in Octave (Octave-Forge collection) and Matlab (Signal Processing Toolbox).

Recall that $\left\vert X_m\right\vert^2\;\longleftrightarrow\;x\star x$ which is circular (cyclic) autocorrelation. To obtain an acyclic autocorrelation instead, we may use zero padding in the time domain. That is, we can replace $x_m$ above by $\hbox{\sc ZeroPad}_{2N}(x_m) = [x_m,0,\ldots,0]$.^8.8Although this fixes the ``wrap-around problem'', the estimator is still biased because its expected value is the true autocorrelation $r_x(l)$ weighted by $N-\vert l\vert$. This bias is equivalent to multiplying the correlation in the ``lag domain'' by a triangular window (also called a ``Bartlett window''). The bias can be removed by simply dividing it out, as in Eq. (8.2), but it is common to retain the Bartlett weighting since it merely corresponds to smoothing the power spectrum (or cross-spectrum) with a sinc$^2$ kernel;^8.9it also down-weights the less reliable large-lag estimates, weighting each lag by the number of lagged products that were summed.

Since $\vert X_m(\omega_k)\vert^2=N\cdot\hbox{\sc DFT}_k({\hat r}_{x_m})$, and since the DFT is a linear operator (§7.4.1), averaging magnitude-squared DFTs $\vert X_m(\omega_k)\vert^2$ is equivalent, in principle, to estimating block autocorrelations ${\hat r}_{x_m}$, averaging them, and taking a DFT of the average. However, this would normally be slower.