scipy.signal.spectrogram() - how to handle gaps in the timeseries data

Question

Looking at the magnitude of variable stars - the dataset is from here:

https://dogwood.physics.mcmaster.ca/Cepheid/URL/MW/BD-10d4669.html

The magnitude plot:

import pandas as pd
from matplotlib import pyplot as plt
from scipy import signal
import numpy as np
plt.rcParams["figure.figsize"] = (16, 9)

df = pd.read_table('BD-10d4669.p.1',
                   sep='   ',
                   engine='python',
                   names=['time', 'mag'])
df.plot(x='time', y='mag');

Clearly there are gaps in the timeline.

But how do I tell signal.spectrogram() that the gaps exist? It seems to me like that function assumes the timeseries could never have any data gaps.

This is the naive spectrogram plot that does not take into account the gaps:

sig = np.array(df['mag'].tolist())
nseg = 20
f, t, Sxx = signal.spectrogram(sig, 1, nperseg=nseg, noverlap=nseg-1)
plt.pcolormesh(t, f, np.log10(Sxx), shading='auto');

I would like the timeline of the spectrogram to match at least partially the timeline of the plot.

I am aware there might be issues due to the bucket size (nperseg), and I am not sure how to handle that either.

Answer 1

OP's time vector is

What I'd do:

1. Treat it as piecewise-lienar, i.e. ignore that the time vector isn't uniformly spaced except for jumps. This should work reasonably - but if greater accuracy is desired, there's a related inquiry.
2. Define "jump" threshold that separates each "segment"
3. Pad each such segment from each side - that is, have left segment's pad meet right segment's pad mid-way -- I recommend reflect, related post
4. Define "discontinuity" or "big jump" threshold, and reflect-zero pad it - that is, pad part of it with reflect, and other with zero. That'd be the second-to-last jump in above pic. The idea is we don't want to impute too much.

These are all "heuristics", but there isn't "the correct" answer for this engineering problem. It's about dealing with missing information.

1 extended

The greatest step size ratio I'd permit is x0.5/x2 of largest step within each piecewise-linear segment. OP's time instances sometimes vary by much more, which invalidates the approximation; below's for the third segment:

With this much variability I'd opt for another method entirely, but if we insist, it becomes an interpolation problem: upsample or downsample until the time vector is sufficiently uniform.

I'd also work with a narrow window, which requires less points to conform to uniformity for a given STFT snapshot. It can then plot as a single 2D representation, but each timestep may encode a different physical time scale, so caution in interpreting is due.

Answer 2

Well it's not just the gaps; your data is also non-uniformly sampled.

Use index_col to use the time column as the index to your dataframe:

df = pd.read_table('BD-10d4669.p.1',
                   sep='   ',
                   engine='python',
                   names=['time', 'mag'],
                   index_col=0)
df.plot(y='mag');

             mag
time            
38964.297  10.92
38965.349  10.86
38966.364  10.52
38968.293  11.38
38972.293  11.09
         ...
48098.300  11.26
48098.399  11.32
48099.308  11.38
48100.325  10.57
48179.272  10.97

[221 rows x 1 columns]

Then you can see the spacing between samples is not exactly the same:

np.diff(df.index)[0:2]
Out[26]: array([1.052, 1.015])

signal.spectrogram only operates on uniformly-sampled data, and you're throwing away the time information with df['mag'].tolist().

You can use signal.lombscargle to see a spectrum of the non-uniformly sampled data, following the example there:

The .note file says:

Columns: 1) HJD-2,400,000
     2) Bpg

P0 = 4.84125 days
P1 = 3.38530 days

HJD is Heliocentric Julian Date, so the sampling times are in days and we want to measure periods of several days.

P0 = 4.84125  # days / cycle
P1 = 3.38530  # days / cycle

f0 = 1/P0  # cycles / day
f1 = 1/P1  # cycles / day

# But lombscargle uses angular frequency, so

w0 = f0 * 2*np.pi  # = 1.29784 radians / day
w1 = f1 * 2*np.pi  # = 1.85602 radians / day

w = np.linspace(0.01, 3, 100000)  # radians / day
pgram = signal.lombscargle(df.index, df['mag'], w, normalize=True)
plt.subplot(2, 1, 1)
plt.plot(df.index, df['mag'], 'b+')
plt.subplot(2, 1, 2)
plt.plot(w, pgram)
plt.show()

This isn't showing the expected peak at angular frequency of ω0 = 1.29784 rad/day, so I'm confused:

I don't know why the units aren't working out. In the example, the time variable is in seconds, and the peak of the periodogram is at 1, which represents 1 rad/sec.

So if the time variable here is in days, I would expect the w variable to be in radians/day.

Anyway, after figuring out why the units don't work (maybe my mistake somewhere), you can try to interpolate the original data to regularly-spaced timestamps with NaN values for the missing days:

df.index = pd.to_datetime(df.index + 2_400_000, unit='D', origin='julian')
df.resample('1D')

But there's so little data with such large gaps between them that I don't think you can get any meaningful spectrogram out of it.

Answer 3

I am experiencing similar issues, hence a recent question on Running window design for irregular or nonuniform time series. One possibility is to invest on uneven or non-uniformly sampled spectrogram, using tools for "least-squares spectral analysis", like Lomb-Scargle periodograms in a moving-window or running frame fashion (which I am redevelopping in Matlab). Yet I still don't know of a canonical response to: how to deal with arbitrary gaps, how to balance power on uneven time-frames.