python - Matlab/Python : Power spectral density of non-uniform time series

Question

我正在尝试查找在不均匀时间测量的信号的功率谱密度。数据看起来像这样:

0 1.55
755 1.58
2412256 2.42
2413137 0.32
2497761 1.19
...

其中第一列是自第一次测量以来的时间(以秒为单位)，第二列是测量值。

目前，使用 Matlab 中的周期图函数，我已经能够通过以下方式估计功率谱密度:

nfft = length(data(:,2));
pxx = periodogram(data(:,2),[],nfft);

现在，为了绘制这个我一直在使用

len = length(pxx);
num = 1:1:len;
plot(num,pxx)

这显然没有将正确的 x 轴放置在功率谱密度上(并产生如下图所示的结果)，而功率谱密度需要位于频率空间中。鉴于数据采样不均匀，我对如何解决这个问题感到困惑。

在估计不均匀采样数据的功率谱密度时，转换到(然后在频率空间中绘制)的正确方法是什么？我也有兴趣从 python/numpy/scipy 的角度解决这个问题，但到目前为止只研究了 Matlab 函数。

Answer 1

我不知道有任何函数可以根据不规则采样数据计算 PSD，因此您需要首先将数据转换为统一采样率。所以第一步是使用 interp1 定期重新采样。

avg_fs = 1/mean(diff(data(:, 1)));
min_time = min(data(:, 1));
max_time = max(data(:, 1));
num_pts = floor((max_time - min_time) * avg_fs);
new_time = (1:num_pts)' / avg_fs;
new_time = new_time - new_time(1) + min_time;
new_x = interp1(data(:, 1), data(:, 2), new_time);

我总是使用 pwelch 来计算 PSD，这是我的方法

nfft = 512; % play with this to change your frequency resolution
noverlap = round(nfft * 0.75); % 75% overlap
window = hanning(nfft);
[Pxx,F] = pwelch(new_x, window, noverlap, nfft, avg_fs);
plot(F, Pxx)
xlabel('Frequency (Hz)')
grid on

您肯定想尝试 nfft，较大的数字将为您提供更高的频率分辨率(频率之间的间距更小)，但 PSD 的噪声会更大。获得高分辨率和低噪声的一个技巧是使窗口小于 nfft。