math - 如何重新采样/重新组合光谱？

Question

在 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 文件大小问题:包括一个光谱位图图像(看起来很笨拙——我想要漂亮的矢量图形!)；
• 或者可能显示一个区域(多边形/置信区间)而不是简单的分段线来指示光谱。

Answer 1

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

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