如何重新采样/重新组合光谱？答案

【问题标题】：How to resample/rebin a spectrum?如何重新采样/重新组合光谱？
【发布时间】：2011-10-03 22:56:31
【问题描述】：

在 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

这很可爱，而且大部分都有效，但是 bin 中心不太正确，并且它不能智能地处理低频区域，频率网格可能会变得更精细地采样。

其他想法：

编写一个函数来确定新的频率分级，然后使用accumarray 进行重新分级。
在进行插值之前对频谱应用平滑滤波器。问题：平滑内核大小必须适应所需的对数平滑。
pwelch 函数接受频率矢量参数f，在这种情况下，它使用 Goetzel 算法计算所需频率的 PSD。也许首先使用对数间隔的频率向量调用pwelch 就足够了。（这或多或少有效率？）
对于 PDF 文件大小问题：包括光谱的位图图像（看起来很笨拙——我想要漂亮的矢量图形！）；
或者可能显示一个区域（多边形/置信区间）而不是简单的分段线来指示光谱。

【问题讨论】：

标签： math matlab pdf fft spectrum

【解决方案1】：

我会让它为我完成工作，并从一开始就给它频率。该文档指出您指定的频率将四舍五入到最近的 DFT bin。这应该不是问题，因为您正在使用结果进行绘图。如果您担心运行时，我会尝试并计时。

如果您想自己重新装箱，我认为您最好编写自己的函数来对每个新箱进行集成。如果您想让您的生活更轻松，您可以做他们所做的事情，并确保您的日志箱与您的线性日志箱共享边界。

【讨论】：

这种方法的问题是平均数是由最低频率（最高分辨率）的bin设置的。我们在高频下得到的分数更少（如所愿），但我们没有在高频下获得更多平均值的优势。

【解决方案2】：

找到的解决方案：https://dsp.stackexchange.com/a/2098/64

简而言之，此问题的一种解决方案是使用与频率相关的变换长度执行 Welch 方法。上面的链接是一个 dsp.SE 答案，其中包含论文引用和示例实现。这种技术的一个缺点是您不能使用 FFT，但是由于要计算的 DFT 点的数量大大减少，所以这不是一个严重的问题。

【讨论】：

【解决方案3】：

如果您想以可变速率（对数）重新采样 FFT，则平滑或低通滤波器内核也需要可变宽度以避免混叠（样本点丢失）。只需为每个绘图点使用不同宽度的同步插值内核（同步宽度大约为本地采样率的倒数）。

【讨论】：