## Power Spectral Density Estimation

Welch's method [85] (or the *periodogram method*
[20]) 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
,
, denote the th block of the
signal
, with denoting the number of blocks.
Then the Welch PSD estimate is given by

where `` '' denotes

*time averaging*across blocks (or ``frames'') of data indexed by . The function

`pwelch`implements Welch's method in Octave (Octave-Forge collection) and Matlab (Signal Processing Toolbox).

Recall that
which is
*circular* (cyclic) autocorrelation. To obtain an *acyclic*
autocorrelation instead, we may use *zero padding* in the time
domain, as described in §8.4.2.
That is, we can replace above by
.^{8.12}Although this fixes the ``wrap-around problem'', the estimator is
still *biased* because its expected value is the true
autocorrelation weighted by . 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 kernel;^{8.13}it also down-weights the less reliable large-lag
estimates, weighting each lag by the number of lagged products that
were summed.

Since
, and since the DFT
is a *linear operator* (§7.4.1), averaging
magnitude-squared DFTs
is *equivalent*, in
principle, to estimating block autocorrelations
, averaging
them, and taking a DFT of the average. However, this would normally
be slower.

We return to power spectral density estimation in Book IV [70] of the music signal processing series.

**Next Section:**

Coherence Function

**Previous Section:**

Correlation Analysis