Matlab/Python: спектральная плотность мощности неоднородных временных рядов

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 для вычисления СПМ, вот как я бы об этом

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.