CODEWITHSUNDEEP
pythonmatlabnumpyspectral-density
Daniel.Reardon
Daniel.Reardon

Reputation: 41

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

I am trying to find the power spectral density of a signal measured at uneven times. The data looks something like this:

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

where the first column is the time since the first measurement (in seconds) and the second column is the value of the measurement.

Currently, using the periodogram function in Matlab, I have been able to estimate the power spectral density by using:

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

Now at the moment, to plot this I have been using

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

Which clearly does not place the correct x-axis on the power spectral density (and yields something like the plot below), which needs to be in frequency space. I am confused about how to go about this given the uneven sampling of the data.

example

What is the correct way to convert to (and then plot in) frequency space when estimating the power spectral density for data that has been unevenly sampled? I am also interested in tackling this from a python/numpy/scipy perspective but have so far only looked at the Matlab function.

Upvotes: 4

Views: 3762

Answers (1)

outdoor_guy
outdoor_guy

Reputation: 146

I am not aware of any functions that calculate a PSD from irregulary sampled data, so you need to convert the data to a uniform sample rate first. So the first step is to use interp1 to resample at regular time intervals.

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);

I always use pwelch for calculating PSD's, here is how I would go about it

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

You will definitely want to experiment with nfft, larger numbers will give you more frequency resolution (smaller spacing between frequencies), but the PSD will be noisier. One trick you can do to get fine resolution and low noise is to make the window smaller than nfft.

Upvotes: 0

Related Questions