# 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
python library
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))
```