# Tutorial: Lots of Plots

This tutorial assumes that you already know how to use Python to read LIGO data files. If you don't, you may want to take a look at the Introductory Tutorial before you go any further.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)
```

## Plot a Time Series

You can plot the time series data like this:```
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()`

.

## Fourier Transform (FFT)

Next, let's try plotting the data in the Fourier domain, using a Fast Fourier Transform (FFT). Because GWOSC data has a great deal of excess noise at low frequencies, spectral leakage is often an issue. To mitigate this, we'll apply a Blackman Window before we take the FFT:```
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')
```

## Power Spectral Density and Amplitude Spectral Density

Plotting in the Fourier domain gives us an idea of the frequency content of the data. Another way to visualize the frequency content of the data is to plot the power spectral density:```
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`

.
## Spectrograms

A spectrogram shows the power spectral density of a signal in a series of time bins. Pyplot has a convienient function for making spectrograms:```
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))
```