I would like to compute a power spectrum in which the frequencies are logarithmically spaced.
In Welch's method there is a trade-off between the frequency resolution of the resulting power spectrum and the number of averages (i.e. error in the result). I would like this trade-off to be dynamic, i.e. do fewer averages for low-frequency points in order to have a finer resolution at low frequency.
Is there a standard way to do this?
I suppose one way would be to initially do pwelch with a very high resolution (low number of averages), and then rebin the resulting spectrum using logarithmic binning.
Answer
I found a paper that addresses this question directly:
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.
No comments:
Post a Comment