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