Free Books    Mathematics of the DFT  

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 $ 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_... ...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\;\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 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.

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