How to resample/rebin a spectrum?

Question

In Matlab, I frequently compute power spectra using Welch's method (pwelch), which I then display on a log-log plot. The frequencies estimated by pwelch are equally spaced, yet logarithmically spaced points would be more appropriate for the log-log plot. In particular, when saving the plot to a PDF file, this results in a huge file size because of the excess of points at high frequency.

What is an effective scheme to resample (rebin) the spectrum, from linearly spaced frequencies to log-spaced frequencies? Or, what is a way to include high-resolution spectra in PDF files without generating excessively large files sizes?

The obvious thing to do is to simply use interp1:

  rate = 16384; %# sample rate (samples/sec)  
  nfft = 16384; %# number of points in the fft

  [Pxx, f] =  pwelch(detrend(data), hanning(nfft), nfft/2, nfft, rate);

  f2 = logspace(log10(f(2)), log10(f(end)), 300);
  Pxx2 = interp1(f, Pxx, f2);

  loglog(f2, sqrt(Pxx2)); 

However, this is undesirable because it does not conserve power in the spectrum. For example, if there is a big spectral line between two of the new frequency bins, it will simply be excluded from the resulting log-sampled spectrum.

To fix this, we can instead interpolate the integral of the power spectrum:

  df = f(2) - f(1);
  intPxx = cumsum(Pxx) * df;                     % integrate
  intPxx2 = interp1(f, intPxx, f2);              % interpolate
  Pxx2 = diff([0 intPxx2]) ./ diff([0 F]);       % difference

This is cute and mostly works, but the bin centers aren't quite right, and it doesn't intelligently handle the low-frequency region, where the frequency grid may become more finely sampled.

Other ideas:

• write a function that determines the new frequency binning and then uses accumarray to do the rebinning.
• Apply a smoothing filter to the spectrum before doing interpolation. Problem: the smoothing kernel size would have to be adaptive to the desired logarithmic smoothing.
• The pwelch function accepts a frequency-vector argument f, in which case it computes the PSD at the desired frequencies using the Goetzel algorithm. Maybe just calling pwelch with a log-spaced frequency vector in the first place would be adequate. (Is this more or less efficient?)
• For the PDF file-size problem: include a bitmap image of the spectrum (seems kludgy--I want nice vector graphics!);
• or perhaps display a region (polygon/confidence interval) instead of simply a segmented line to indicate the spectrum.

Answer 1

I would let it do the work for me and give it the frequencies from the start. The doc states the freqs you specify will be rounded to the nearest DFT bin. That shouldn't be a problem since you are using the results to plot. If you are concerned about the runtime, I'd just try it and time it.

If you want to rebin it yourself, I think you're better off just writing your own function to do the integration over each of your new bins. If you want to make your life easier, you can do what they do and make sure your log bins share boundaries with your linear ones.

Answer 2

Solution found: https://dsp.stackexchange.com/a/2098/64

Briefly, one solution to this problem is to perform Welch's method with a frequency-dependent transform length. The above link is to a dsp.SE answer containing a paper citation and sample implementation. A disadvantage of this technique is that you can't use the FFT, but because the number of DFT points being computed is greatly reduced, this is not a severe problem.