Spectrograms and Time-Frequency Tracks¶
One of the most common forms of data visualization for gravitational wave data is the spectrogram, which shows a time-frequency representation of the data. These are often useful for getting a rough sense of the signals characteristics, and are also useful for identifying any detector characterization issues which may overlap the signal. You can also add predictions of the signal's time-frequency track to these, so you can see how the prediction from the inferred posterior lines up with the data you see. This tutorial shows how to make a spectrogram from the bilby IFO data, and add a signal TF track to it.
Preliminaries¶
# For the purposes of notebook readability, we'll surpress the various logging messages from gwpy, lal, bilby, etc
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import logging
logging.disable(logging.CRITICAL)
import bilby
import ptfp
import numpy as np
from gwpy.timeseries import TimeSeries
import matplotlib.pyplot as plt
ptfp.configuration.set_better_matplotlib_defaults(fontsize=12)
%matplotlib inline
Loading in Bilby Results¶
Once a bilby_pipe run has completed, the posteriors are easily accessible with bilby.core.result.read_in_result(), but other parts of the analyses may be more difficult to load in.
When all you want to make are corner plots this is fine, but we may often be interested in more complicated sorts of figures using the data (that's what ptfp is for!).
So, we want to load everything we'll need in, something which ptfp lets you do in one line!
The catch is that this requires putting in the original run directory, which may be hard to find.
This is something worth improving on in the future, but due to how bilby_pipe works it's a much more complicated task to rebuild from the metadata alone.
On CIT, I've done an example bilby_pipe run on GW150914 to use for these example pages.
# If these were injections and we wanted the Nth injection (i.e. "_data{N}", we could pass "data_case=N")
# This defaults to using the IFOs H1L1 - this is used to search for the correct file names
# For 3 detector events do H1L1V1, or for a single detector do e.g. H1
result, data_dump, analysis, likelihood = ptfp.bilby_pipe_interface.get_all_analysis_components("/home/rhiannon.udall/PlottingToolsForPEExampleRun/GW150914_bilby_output")
Making the Spectrogram¶
Spectrograms require timeseries data. There will generally be two ways to obtain this: reading the frame files directly, or reading it out of bilby's data dump. While in principle these are essentially the same, bilby does do some peculiar data conditioning that may impact gwpy's normal operations, so it's broadly advisable to load the data directly. We'll do this using data copied directly to the cluster.
lho_timeseries = TimeSeries.read("/home/rhiannon.udall/PlottingToolsForPEExampleRun/H-H1_GWOSC_4KHZ_R1-1126259447-32.gwf", channel="H1:GWOSC-4KHZ_R1_STRAIN")
llo_timeseries = TimeSeries.read("/home/rhiannon.udall/PlottingToolsForPEExampleRun/L-L1_GWOSC_4KHZ_R1-1126259447-32.gwf", channel="L1:GWOSC-4KHZ_R1_STRAIN")
Now, the actual spectrogram; we do this with a wrapper around gwpy methods.
fig, (lho_ax, llo_ax) = plt.subplots(2, 1, tight_layout=True)
ptfp.spectrograms.plot_spectrogram_from_timeseries(
lho_timeseries,
q_value=16,
f_max=400,
frequency_resolution=0.01,
time_resolution=0.001,
fig_and_ax=(fig, lho_ax),
start_time=result.start_time + 5.5,
end_time=result.start_time + 6.5,
colormap='binary'
)
lho_ax.set_title("Hanford")
ptfp.spectrograms.plot_spectrogram_from_timeseries(
llo_timeseries,
q_value=16,
f_max=400,
frequency_resolution=0.01,
time_resolution=0.001,
fig_and_ax=(fig, llo_ax),
start_time=result.start_time + 5.5,
end_time=result.start_time + 6.5,
colormap='binary'
)
llo_ax.set_title("Livingston")
Text(0.5, 1.0, 'Livingston')
As you can see, many formatting headaches are taken care of automatically with this function.
However, there are still a few less-than-perfect things.
One that stands out are the times, with a somewhat weird reference including the epoch of the segment.
ptfp always maintains times at the true values for consistency, but provides a trick to overwrite the times that appear on an axis:
ptfp.time_domain_plots.force_axis_time_offset(
lho_ax,
new_ticks=np.linspace(-0.5, 0.5, 11),
tick_offset=result.start_time+6
)
ptfp.time_domain_plots.force_axis_time_offset(
llo_ax,
new_ticks=np.linspace(-0.5, 0.5, 11),
tick_offset=result.start_time+6
)
<Axes: title={'center': 'Livingston'}, xlabel='Time [s]', ylabel='Frequency [Hz]'>
fig
Making the Time-Frequency Track¶
For ease of use, ptfp splits up the process into two steps.
First, we invoke a function which draws $N$ samples from a posterior, then computes $N$ $f(t)$ series for it.
Then, we use another function to plot these onto the spectrogram.
The first function can take a bit of time (around a minute, depending on the number of samples and the waveform), and so it's convenient to only perform that step once, then play with the plot settings separately.
lho_time_frequency_posterior = ptfp.time_frequency_tracks.time_frequency_posterior_from_bilby_components(
result.posterior,
data_dump.interferometers[0],
likelihood.waveform_generator,
relative_time_of_track_start=5.5,
relative_time_of_track_end=6.05,
number_of_samples=100
)
llo_time_frequency_posterior = ptfp.time_frequency_tracks.time_frequency_posterior_from_bilby_components(
result.posterior,
data_dump.interferometers[1],
likelihood.waveform_generator,
relative_time_of_track_start=5.5,
relative_time_of_track_end=6.04,
number_of_samples=100
)
ptfp.time_frequency_tracks.plot_time_frequency_posterior(
lho_time_frequency_posterior,
data_dump.interferometers[0].time_array,
color='C0',
time_frequency_label='Signal',
fig_and_ax = (fig, lho_ax)
)
ptfp.time_frequency_tracks.plot_time_frequency_posterior(
llo_time_frequency_posterior,
data_dump.interferometers[1].time_array,
color='C0',
time_frequency_label='Signal',
fig_and_ax = (fig, llo_ax)
)
/home/rhiannon.udall/.conda/envs/bilby-MRs-testing/lib/python3.11/site-packages/ptfp/time_frequency_tracks.py:380: RuntimeWarning: All-NaN slice encountered np.nanmedian(time_frequency_posterior, axis=0), /home/rhiannon.udall/.conda/envs/bilby-MRs-testing/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1563: RuntimeWarning: All-NaN slice encountered return function_base._ureduce(a,
<Axes: title={'center': 'Livingston'}, xlabel='Time [s]', ylabel='Frequency [Hz]'>
fig
