How to reprogram / rebuild the spectrum?

Question

In Matlab, I often compute power spectra using the pwelch method, which I then show on a graph in a logarithmic log. The frequencies evaluated on pwelch are evenly distributed, but log-spaced points will be more suitable for a log-log plot. In particular, when saving a graph to a PDF file, this leads to a huge file size due to excess points at a high frequency.

What is an efficient scheme for re-sampling a (re) spectrum from linearly spaced frequencies to a frequency with an interval in the log? Or, what is the way to include high-resolution spectra in PDF files without generating too large file sizes?

The obvious thing to do is just 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 since it does not retain strength in the spectrum. For example, if there is a large spectral line between two new frequency cells, it will simply be excluded from the resulting spectrum with a discrete spectrum.

To fix this, we can instead interpolate the integral from 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 nice and basically works, but the centers of the bins are not quite right, and he does not know how to handle the low-frequency region, where the frequency grid can be more finely rejected.

Other ideas:

• write a function that defines a new binary frequency value, and then uses accumarray to replay.
• Apply a smoothing filter to the spectrum before interpolation. Problem: The size of the smoothing kernel must be adaptive to the desired logarithmic smoothing.
• The pwelch function takes the argument of the frequency vector f , in which case it calculates the PSD at the required frequencies using the Goetzel algorithm. Maybe just calling pwelch with a frequency vector with an interval in the log space will be primarily adequate. (Is it more or less effective?)
• For a problem with the PDF file size: include a spectrum bitmap (it seems kludgy - I need good vector graphics!);
• or perhaps display an area (polygon / confidence interval) instead of just a segmented line to indicate the spectrum.

Answer 1

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

In short, one solution to this problem is to implement the Welch method with a frequency-dependent conversion length. The above link refers to dsp.SE's answer, containing a citation and selection in the form of paper. The disadvantage of this method is that you cannot use FFT, but since the number of calculated DFT points is significantly reduced, this is not a serious problem.

Answer 2

I would let him do the work for me and give her frequencies from the start. The document states that the frequencies you specify will be rounded to the nearest DFT bean. This should not be a problem as you use the results to build. If you are concerned about the runtime, I will just try the time.

If you want to rebuild it yourself, I think you better just write your own function to do the integration for each of your new bins. If you want to make your life easier, you can do what they do and make sure that your log files have borders with your linear ones.

Answer 3

If you want to reprogram the FFT with a variable speed (logarithmically), then the smoothing core or low-pass filter should also be of variable width to avoid smoothing (loss of sample points). Just use a different Sync interpolation width for each point in the graph (the synchronization width is approximately equal to the inverse sampling frequency).