Getting started with timeseries

This notebook will walk you through some of the basic methods in PTSA’s TimeSeries class. We’ll see how to construct an instance of TimeSeries, how to save and load data, as well as resampling and filtering TimeSeries.

# relevant imports
import numpy as np
from matplotlib import pyplot as plt
%matplotlib notebook
from ptsa.data.timeseries import TimeSeries

Let’s start with creating some data

In real applications, you will most likely have your own timeseries data for analysis. For the purpose of illustrating the functionalities of the timeseries class, we will construct sinusoids as our timeseries data. Let’s create an array of 5000 data points, or samples. Suppose the sampling rate is 10Hz, this means that our timeseires is 5000/10=500 seconds long.

num_points = 5000
sample_rate = 10.

# We can specify the timestamps for each data point, from 0s to 500s.
t = np.linspace(1, num_points, num_points) / sample_rate

# Let's create two noisy sinusoids with different frequencies.
frequency1 = .5 # 1 cycle every 2 seconds
frequency2 = .1 # 1 cycle every 10 seconds
data1 = np.sin(2*np.pi*frequency1*t) + np.random.uniform(-0.5, 0.5, num_points)
data2 = np.sin(2*np.pi*frequency2*t) + np.random.uniform(-0.5, 0.5, num_points)

Let’s check our timepoints.

print ('First 5 timestamps: ', t[:5])
print ('Last 5 timestamps: ', t[-5:])

First 5 timestamps:  [ 0.1  0.2  0.3  0.4  0.5]
Last 5 timestamps:  [ 499.6  499.7  499.8  499.9  500. ]

We can also visualize the timeseries using matplotlib.

plt.figure(figsize=[10,2])
plt.plot(t, data1, label='%sHz'%str(frequency1))
plt.plot(t, data2, label='%sHz'%str(frequency2))

plt.legend()

<matplotlib.legend.Legend at 0x2b7884cc0518>

As we zoom in to the data array, the random noise we added to the sinusoids becomes clear.

plt.figure(figsize=[10, 2])
plt.plot(t[500:1000], data1[500:1000], label='%sHz'%str(frequency1))
plt.plot(t[500:1000], data2[500:1000], label='%sHz'%str(frequency2))
plt.legend()

<matplotlib.legend.Legend at 0x2b78be4cc908>

Create a TimeSeries object

The TimeSeries class is a convenient wrapper of xarray that offers basic functionalities for timeseries analysis. Although we focus our analysis here in the context of timeseries data, many of the following examples can be easily extended to non-timeseries, multidimensional data. To create a TimeSeries object, we simply need to construct dimensions and the corresponding coordinates in each dimension.

# Let's stack the two time-series data arrays.
data = np.vstack((data1, data2))

# and construct the TimeSeries object
ts = TimeSeries(data,
                dims=('frequency', 'time'),
                coords={'frequency':[frequency1, frequency2],
                        'time':t,
                        'samplerate':sample_rate})
print (ts)

<xarray.TimeSeries (frequency: 2, time: 5000)>
array([[-0.02586 ,  0.497876,  0.336389, ..., -0.680535, -0.378594,  0.029344],
       [-0.016682, -0.123327, -0.285444, ..., -0.406314,  0.222235, -0.272452]])
Coordinates:
  * frequency   (frequency) float64 0.5 0.1
  * time        (time) float64 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 ...
    samplerate  float64 10.0

TimeSeries also has a convenient plotting method, inherited from xarray.

Note: since the frequency dimension has float coordinates, we want to be careful for exact float comparisons. Thus, instead of using ts.sel(frequency=frequency1), we use ts.sel(frequency=ts.frequency[0]). See more about the .sel() method in later sections.

plt.figure(figsize=[10,2])
ts.sel(frequency=ts.frequency[0]).plot()

[<matplotlib.lines.Line2D at 0x2b78be58e320>]

Saving and loading your data

TimeSeries objects can be easily saved and loaded in HDF5 format.

Note: You will need to install h5py to use these methods.

ts

<xarray.TimeSeries (frequency: 2, time: 5000)>
array([[-0.02586 ,  0.497876,  0.336389, ..., -0.680535, -0.378594,  0.029344],
       [-0.016682, -0.123327, -0.285444, ..., -0.406314,  0.222235, -0.272452]])
Coordinates:
  * frequency   (frequency) float64 0.5 0.1
  * time        (time) float64 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 ...
    samplerate  float64 10.0

# use to_hdf to save data
fname = 'my_ts_data.h5'
ts.to_hdf(fname)

# use from_hdf to load data
ts = TimeSeries.from_hdf(fname)
print (ts)

<xarray.TimeSeries (frequency: 2, time: 5000)>
array([[-0.02586 ,  0.497876,  0.336389, ..., -0.680535, -0.378594,  0.029344],
       [-0.016682, -0.123327, -0.285444, ..., -0.406314,  0.222235, -0.272452]])
Coordinates:
  * frequency   (frequency) float64 0.5 0.1
  * time        (time) float64 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 ...
    samplerate  float64 10.0

Indexing your data

We can used xarray’s .sel() method to index our data, or give the dimension name instead of the dimension index to other methods.

Click here to learn more about xarray methods.

# Select the data from 100s to 200s
ts.sel(time=(ts.time>100)&(ts.time<200))

<xarray.TimeSeries (frequency: 2, time: 999)>
array([[-0.111702,  0.824685,  0.499759, ..., -0.377922, -0.517297,  0.072835],
       [ 0.368357, -0.372158,  0.431717, ..., -0.62199 , -0.475753,  0.178974]])
Coordinates:
  * frequency   (frequency) float64 0.5 0.1
  * time        (time) float64 100.1 100.2 100.3 100.4 100.5 100.6 100.7 ...
    samplerate  float64 10.0

# mean over the time dimension
ts.mean('time')

<xarray.TimeSeries (frequency: 2)>
array([-0.003609,  0.005354])
Coordinates:
  * frequency   (frequency) float64 0.5 0.1
    samplerate  float64 10.0

Matching coordinates

We can also place the time dimension as the first dimension. The coordinates in xarray are very useful and are used to keep track of the dimensions of the data. Therefore, the exact shape of the data or the order of the dimensions can vary.

ts_transpose = ts.T

print (ts_transpose)

<xarray.TimeSeries (time: 5000, frequency: 2)>
array([[-0.02586 , -0.016682],
       [ 0.497876, -0.123327],
       [ 0.336389, -0.285444],
       ...,
       [-0.680535, -0.406314],
       [-0.378594,  0.222235],
       [ 0.029344, -0.272452]])
Coordinates:
  * frequency   (frequency) float64 0.5 0.1
  * time        (time) float64 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 ...
    samplerate  float64 10.0

We can try adding the original timeseries and its “transposed” version together. Because their coordinates match, operations can be performed between them even though the data arrays are technically of different shapes. Note that because the transposed version of the data has the same value at every coordinate, by adding them together we’re effectively doubling the values.

results = ts + ts_transpose

plt.figure(figsize=[10, 3])

# we'll only plot the first 250 samples to see things clearly
plt.plot(ts.time[:250], ts.sel(frequency=ts.frequency[0])[:250], label='original')
plt.plot(results.time[:250], results.sel(frequency=ts.frequency[0])[:250], label='original+transposed')
plt.ylim([-4,4])
plt.legend(ncol=2)

<matplotlib.legend.Legend at 0x2b78be5e9a20>

Without the coordinates, the data cannot be added because their shapes wouldn’t match.

try:
    print (data+data.T)
except Exception as e:
    print ('Error: ' + str(e))

Error: operands could not be broadcast together with shapes (2,5000) (5000,2)

Note that sometimes, this feature can produce undesired results. Suppose that we want to subtract a subset of our data from another subset of our data, because xarray matches coordinates, it will drop all unmatched coordinates and return 0 at each matched coordinate.

# Recall that our data contains 5000 points.
# Let's select the first 3000 points and the last 3000 points and subtract one from the other.
# The result is all zeros because unmatched coordinates are dropped,
# and the value at each of the 1000 overlapping coordinates is subtracted from itself.

early_ts = ts.sel(time=ts.time[:3000])
late_ts = ts.sel(time=ts.time[2000:])

print (late_ts - early_ts)

<xarray.TimeSeries (frequency: 2, time: 1000)>
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])
Coordinates:
  * time        (time) float64 200.1 200.2 200.3 200.4 200.5 200.6 200.7 ...
  * frequency   (frequency) float64 0.5 0.1
    samplerate  float64 10.0

# However, if we only subtract the underlying data arrays,
# it will return the actual difference of the early subset and the late subset of our timeseries.

print (late_ts.values - early_ts.values)

[[ 0.28978477  0.07390804  0.070932   ..., -0.17307536 -0.16645912
  -0.244184  ]
 [ 0.15750319  0.59501017  0.48198941 ...,  0.04690397  0.33666706
  -0.21464846]]

Resampling your data

We can resample the data to a specific samplerate.

original = ts.sel(frequency=ts.frequency[0]).sel(time=ts.time<50.0)
downsampled = original.resampled(resampled_rate=2.0)

plt.figure(figsize=[10, 4])
plt.plot(original.time, original, color='C0', label='10.0Hz samplerate (original)')
plt.plot(downsampled.time, downsampled, color='C6', label='2.0Hz samplerate (downsampled)')
plt.legend()

<matplotlib.legend.Legend at 0x2b78befd2d68>

However, as we downsample to lower frequencies, we lose precision of the signal.

Read more about the relationship between sampling rate and observed frequencies here.

original = ts.sel(frequency=ts.frequency[0]).sel(time=ts.time<50.0)
downsampled = original.resampled(resampled_rate=1.0)

plt.figure(figsize=[10, 4])
plt.plot(original.time, original, color='C0', label='10.0Hz samplerate (original)')
plt.plot(downsampled.time, downsampled, color='C6', label='1.0Hz samplerate (downsampled)')
plt.legend()

<matplotlib.legend.Legend at 0x2b78bfcc0b00>

# TimeSeries updates the samplerate for you
print (original.samplerate)
print (downsampled.samplerate)

<xarray.TimeSeries 'samplerate' ()>
array(10.0)
Coordinates:
    frequency   float64 0.5
    samplerate  float64 10.0
<xarray.TimeSeries 'samplerate' ()>
array(1.0)
Coordinates:
    samplerate  float64 1.0

Filtering your data

Let’s create a new TimeSeries object composed of sinunoids at three different frequencies. We’ll show how to manipulate these data using different filtering methods.

freq1 = 3.2
freq2 = 1.6
freq3 = 0.2

data1 = np.sin(2*np.pi*freq1*t)
data2 = np.sin(2*np.pi*freq2*t)
data3 = np.sin(2*np.pi*freq3*t)

# our data are simply the sum of the three sinusoids
data = data1 + data2 + data3

ts = TimeSeries(data, dims=('time'), coords={'time':t, 'samplerate':sample_rate})
print (ts)

<xarray.TimeSeries (time: 5000)>
array([  1.874488e+00,   3.830037e-01,   2.447679e-01, ...,  -3.830037e-01,
        -1.874488e+00,   9.821934e-14])
Coordinates:
  * time        (time) float64 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 ...
    samplerate  float64 10.0

# Let's plot the first 200 samples of the data
fig, ax = plt.subplots(4, figsize=[10, 4], sharex=True, sharey=True)
ax[0].plot(t[:200], data1[:200])
ax[1].plot(t[:200], data2[:200])
ax[2].plot(t[:200], data3[:200])
ax[3].plot(t[:200], ts[:200])

ax[0].set_ylabel('3.2Hz')
ax[1].set_ylabel('1.6Hz')
ax[2].set_ylabel('0.2Hz')
ax[3].set_ylabel('sum')
ax[3].set_xlabel('Time(s)')

<matplotlib.text.Text at 0x2b78bfeb2ba8>

We will use three different filters to filter out each component.

1. To filter out the component with the highest frequency (3.2Hz), we’ll use a lowpass filter. A lowpass filter perserves any frequency that is lower than the given frequency.

2. To filter out the component with the lowest frequency (0.2Hz), we’ll use a highpass filter. A highpass filter perserves any frequency that is higher than the given frequency.

3. To filter out the component with the middle frequency (1.6Hz), we’ll use a bandstop filter. A bandstop filter perserves any frequency that is outside of the given frequency range.

Note that these filters suffer from edge effects at both ends of the timeseries.

# lowpass filter
filtered_data = ts.filtered(3.0, filt_type='lowpass', order=4)

fig, ax = plt.subplots(4, figsize=[10, 5], sharex=True, sharey=True)

ax[0].plot(t[:200], ts[:200]) # origianl timeseries
ax[1].plot(t[:200], filtered_data[:200]) # losspass filtered
ax[2].plot(t[:200], (data2+data3)[:200]) # what we should get (mid + low frequencies)
ax[3].plot(t[:200], (data2+data3-filtered_data)[:200]) # the difference between what we should get and what we got should be close to zero

ax[0].set_ylabel('unfiltered')
ax[1].set_ylabel('filtered')
ax[2].set_ylabel('mid + low')
ax[3].set_ylabel('(mid+low)-filtered')
ax[3].set_xlabel('Time(s)')

/home1/liyuxuan/anaconda3/envs/greatest_environment/lib/python3.6/site-packages/ptsa/data/timeseries.py:416: UserWarning: The filtered method is not very flexible. Consider using filters in ptsa.data.filters instead.
  warnings.warn("The filtered method is not very flexible. "

<matplotlib.text.Text at 0x2b78c8143400>

# highpass filter
filtered_data = ts.filtered(0.5, filt_type='highpass', order=4)

fig, ax = plt.subplots(4, figsize=[10, 5], sharex=True, sharey=True)

ax[0].plot(t[:200], ts[:200]) # origianl timeserids
ax[1].plot(t[:200], filtered_data[:200]) # highpass filtered
ax[2].plot(t[:200], (data2+data1)[:200]) # what we should get (mid + high frequencies)
ax[3].plot(t[:200], (data2+data1-filtered_data)[:200]) # the difference between what we should get and what we got should be close to zero

ax[0].set_ylabel('unfiltered')
ax[1].set_ylabel('filtered')
ax[2].set_ylabel('mid + high')
ax[3].set_ylabel('(mid+high)-filtered')
ax[3].set_xlabel('Time(s)')

/home1/liyuxuan/anaconda3/envs/greatest_environment/lib/python3.6/site-packages/ptsa/data/timeseries.py:416: UserWarning: The filtered method is not very flexible. Consider using filters in ptsa.data.filters instead.
  warnings.warn("The filtered method is not very flexible. "

<matplotlib.text.Text at 0x2b78c92c8940>

# bandstop filter
filtered_data = ts.filtered([1.4, 1.8], filt_type='stop', order=4)

fig, ax = plt.subplots(4, figsize=[10, 5], sharex=True, sharey=True)

ax[0].plot(t[:200], ts[:200]) # origianl timeserids
ax[1].plot(t[:200], filtered_data[:200]) # bandstop filtered
ax[2].plot(t[:200], (data1+data3)[:200]) # what we should get (high + low frequencies)
ax[3].plot(t[:200], (data1+data3-filtered_data)[:200]) # the difference between what we should get and what we got should be close to zero

ax[0].set_ylabel('unfiltered')
ax[1].set_ylabel('filtered')
ax[2].set_ylabel('high + low')
ax[3].set_ylabel('(high+low)-filtered')
ax[3].set_xlabel('Time(s)')

/home1/liyuxuan/anaconda3/envs/greatest_environment/lib/python3.6/site-packages/ptsa/data/timeseries.py:416: UserWarning: The filtered method is not very flexible. Consider using filters in ptsa.data.filters instead.
  warnings.warn("The filtered method is not very flexible. "

<matplotlib.text.Text at 0x2b78c94bc9e8>

There you have it! You’re now ready to get rolling with PTSA’s timeseries.