如何重新取样/重新组合频谱？

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

这很可爱并且大多数都可以工作,但是箱子中心不是很正确,并且它不能智能地处理低频区域,其中频率网格可以变得更精细地采样.

其他想法:

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

Answer 1

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

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