如何对频谱进行重新采样/重新分类？
在 Matlab 中，我经常使用 Welch 方法 (pwelch) 计算功率谱，然后将其显示在对数图上。 pwelch 估计的频率是等距的，但对数间隔的点更适合双对数图。特别是，将绘图保存为 PDF 文件时，由于高频点过多，会导致文件大小过大。
从线性间隔频率到对数间隔频率重新采样（重新分类）频谱的有效方案是什么？或者，在 PDF 文件中包含高分辨率频谱而不生成过大文件的方法是什么？尺寸？
显而易见的做法是简单地使用 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));
然而，这是不可取的，因为它不会节省频谱功率。例如，如果两个新频率仓之间存在较大的谱线，则它将简单地从生成的对数采样频谱中排除。
为了解决这个问题，我们可以对功率谱的积分进行插值：
df = f(2) - f(1);
intPxx = cumsum(Pxx) * df; % integrate
intPxx2 = interp1(f, intPxx, f2); % interpolate
Pxx2 = diff([0 intPxx2]) ./ diff([0 F]); % difference
这很可爱并且大部分都有效，但是箱中心不太正确，并且它不能智能地处理低频区域，其中频率网格可以变得更精细采样。
其他想法：
- 编写一个函数来确定新的频率分级，然后使用 Accumarray 进行重新分级。
- 在进行插值之前对频谱应用平滑滤波器。问题：平滑核大小必须适应所需的对数平滑。
pwelch函数接受频率向量参数f，在这种情况下，它使用 Goetzel 算法计算所需频率的 PSD。也许首先使用对数间隔频率向量调用 pwelch 就足够了。 （这效率更高还是更低？）- 对于 PDF 文件大小问题：包含光谱的位图图像（看起来很笨拙 - 我想要漂亮的矢量图形！）；
- 或者可能显示一个区域（多边形/置信区间），而不是简单地用分段线来指示频谱。
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
accumarrayto 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
pwelchfunction accepts a frequency-vector argumentf, in which case it computes the PSD at the desired frequencies using the Goetzel algorithm. Maybe just callingpwelchwith 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.
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（3）
我会让它为我工作并从一开始就给它频率。该文档指出您指定的频率将四舍五入到最近的 DFT bin。这应该不是问题，因为您正在使用结果来绘图。如果您担心运行时间，我会尝试并计时。
如果您想自己重新装箱，我认为您最好编写自己的函数来对每个新箱进行集成。如果您想让您的生活更轻松，您可以做他们所做的事情，并确保您的日志箱与线性日志箱共享边界。
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.
找到的解决方案： https://dsp.stackexchange.com/a/2098/64
简而言之，这是一个解决方案问题是使用频率相关的变换长度执行韦尔奇方法。上面的链接是 dsp.SE 答案，其中包含论文引用和示例实现。此技术的缺点是无法使用 FFT，但由于计算的 DFT 点数大大减少，因此这不是一个严重的问题。
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.
如果您想以可变速率（对数）重新采样 FFT，则平滑或低通滤波器内核也需要具有可变宽度，以避免混叠（采样点丢失）。只需对每个绘图点使用不同宽度的同步插值内核（同步宽度大约是本地采样率的倒数）。
If you want to resample an FFT at a variable rate (logarithmically), then the smoothing or low pass filter kernel will need to be variable width as well to avoid aliasing (loss of sample points). Just use a different width Sync interpolation kernel for each plot point (Sync width approximately the reciprocal of the local sampling rate).