How can I compute a log-spaced power spectrum?
I found a paper that addresses this question directly:


Michael Tröbs and Gerhard Heinzel, "Improved spectrum estimation
 from digitized time series on a logarithmic frequency axis," in
 Measurement, vol 39 (2006), pp 120-129. (Or free pre-print.)


The first few figures in the paper nicely illustrate the problem that this algorithm solves, and the references contain a useful bibliography of other approaches (constant-Q transform, tempered Fourier transform, a survey article, etc).

Their approach is not to re-bin the output of an existing FFT-based power spectrum estimation, but to only compute the discrete Fourier transform at the (logarithmically-spaced) frequencies of interest.  For each frequency to be estimated, they basically implement Welch's algorithm, but with a transform length (and hence also, number of averages) specifically chosen for each frequency. The computation of each frequency bin uses the entire time series, but segmented differently. The results have the desirable property that the resolution (bin width) is a smooth function of frequency, and the results can be calibrated as either a power spectral density or a power spectrum.

Matlab implementation here: https://github.com/tobin/lpsd


Disclosure: The authors of this paper are at the same institution as me.