Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search


Power Spectral Density Estimation

Welch's method [87] (or the periodogram method [21]) 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

$\displaystyle {\hat R}_x(\omega_k) = \frac{1}{M}\sum_{m=0}^{M-1}\left\vert DFT_k(x_m)\right\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\;\leftrightarrow\;x\star x$ 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 $ x_m$ above by $ \hbox{\sc CausalZeroPad}_{2N-1}(x_m) =
[x_m,0,\ldots,0]$ .8.12Although 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.13it 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 operator7.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.

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


Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

[How to cite this work]  [Order a printed hardcopy]  [Comment on this page via email]

``Mathematics of the Discrete Fourier Transform (DFT), with Audio Applications --- Second Edition'', by Julius O. Smith III, W3K Publishing, 2007, ISBN 978-0-9745607-4-8.
Copyright © 2014-04-06 by Julius O. Smith III
Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
CCRMA