This tutorial will show you how to make a few common plots with time series data, including a Fourier domian representation, an amplitude spectral density, and a spectrogram.
Create a file called lotsofplots.py
and start with
the following code:
import numpy as np
import h5py
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import readligo as rl
Notice that readligo
just refers to the
sample reader functions
described in step 4, which
read data from a LIGO data file. Of course, you can also
write your own code for this.
We also added the mlab
module from matplotlib
,
which includes a function to calculate a power spectral density (PSD).
To start, we'll read the strain data and data quality information
from the
LIGO data file
that you downloaded in
step 1 of the introductory
tutorial:
fileName = 'fileName = 'H-H1_LOSC_4_V1-815411200-4096.hdf5'
strain, time, channel_dict = rl.loaddata(fileName)
ts = time[1] - time[0] #-- Time between samples
fs = int(1.0 / ts) #-- Sampling frequency
To make some example plots, let's pick a segment of "good data", and then
work with just the first 16 seconds of data in the segment:
segList = rl.dq_channel_to_seglist(channel_dict['DEFAULT'], fs)
length = 16 # seconds
strain_seg = strain[segList[0]][0:(length*fs)]
time_seg = time[segList[0]][0:(length*fs)]
In the example, we'll use pyplot.suplot
so that all of the
plots will be on one figure. You may prefer to use pyplot.figure()
instead to make separate figures. To set up a grid of
6 figures, with a little extra space for axis labels:
fig = plt.figure(figsize=(10,10))
fig.subplots_adjust(wspace=0.3, hspace=0.3)
plt.subplot(321)
plt.plot(time_seg - time_seg[0], strain_seg)
plt.xlabel('Time since GPS ' + str(time_seg[0]))
plt.ylabel('Strain')
If you don't see the plot, try the command plt.show()
.
window = np.blackman(strain_seg.size)
windowed_strain = strain_seg*window
freq_domain = np.fft.rfft(windowed_strain) / fs
freq = np.fft.rfftfreq(len(windowed_strain))*fs
The vector returned in freq_domain
is complex.
Let's just plot the magnitude of this vector:
plt.subplot(322)
plt.loglog( freq, abs(freq_domain))
plt.axis([10, fs/2.0, 1e-24, 1e-18])
plt.grid('on')
plt.xlabel('Freq (Hz)')
plt.ylabel('Strain / Hz')
Pxx, freqs = mlab.psd(strain_seg, Fs = fs, NFFT=fs)
plt.subplot(323)
plt.loglog(freqs, Pxx)
plt.axis([10, 2000, 1e-46, 1e-36])
plt.grid('on')
plt.ylabel('PSD')
plt.xlabel('Freq (Hz)')
You can also plot an amplitude spectra density the same way,
just by taking the square root, np.sqrt()
, of Pxx
.
NFFT = 1024
window = np.blackman(NFFT)
plt.subplot(325)
spec_power, freqs, bins, im = plt.specgram(strain_seg, NFFT=NFFT, Fs=fs,
window=window)
Mainly, this spectrogram shows what we already knew from the PSD:
there is a lot more power at very low frequencies than high frequencies.
If we are looking for time variation in the data, it can be
helpful to normalize each frequency bin by the typical power
at that frequency. We can do this, and then re-plot the
spectrogram, using the colormesh
function:
med_power = np.zeros(freqs.shape)
norm_spec_power = np.zeros(spec_power.shape)
index = 0
for row in spec_power:
med_power[index] = np.median(row)
norm_spec_power[index] = row / med_power[index]
index += 1
ax = plt.subplot(326)
ax.pcolormesh(bins, freqs, np.log10(norm_spec_power))