Development documentation · 0.2.3 30c2f8ba · Intro examples tested with 0.2.3 · Version details · Known limitations

ASD Analysis: Pipeline#

This tutorial demonstrates how to use SegmentTable for a batch ASD analysis pipeline. We will crop data for each segment, compute ASDs, and visualize the variation.

The segment windows are still gwpy.segments.Segment, and the waveform payload in each row is a GWpy TimeSeries. gwexpy extends that GWpy base-class setup by placing the rows in SegmentTable and applying crop() and asd() across the whole table at once. For the core relationship between GWpy base classes and gwexpy extensions, see SegmentTable: Basics.

import warnings


import warnings

with warnings.catch_warnings():

    import numpy as np
    np.random.seed(42)
    from gwpy.segments import Segment
    from gwpy.timeseries import TimeSeries

    from gwexpy.table import SegmentTable

    def get_synthetic_data(t0):
        return TimeSeries(np.random.randn(1024), sample_rate=64, t0=t0)

    segs = [Segment(i*16, i*16+16) for i in range(4)]
    st = SegmentTable.from_segments(segs)
    st.add_series_column("raw", data=[get_synthetic_data(seg[0]) for seg in segs], kind="timeseries")
    st

/home/runner/micromamba/envs/gwexpy/lib/python3.11/site-packages/gwpy/time/_ligotimegps.py:42: UserWarning: Wswiglal-redir-stdio:

SWIGLAL standard output/error redirection is enabled in IPython.
This may lead to performance penalties. To disable locally, use:

with lal.no_swig_redirect_standard_output_error():
    ...

To disable globally, use:

lal.swig_redirect_standard_output_error(False)

Note however that this will likely lead to error messages from
LAL functions being either misdirected or lost when called from
Jupyter notebooks.

To suppress this warning, use:

import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal

  from lal import LIGOTimeGPS

Crop and ASD#

We can use sugar APIs like crop() and asd() to process all segments at once.

# Cut the data to each segment span so every ASD is computed from intervals that share the same operating state.
st_cropped = st.crop("raw", out_col="cropped")
# Estimate ASD per segment to compare stationary noise floors without smearing glitches or state changes across the whole run.
st_asd = st_cropped.asd("cropped", out_col="asd", fftlength=2.0)
st_asd.display()

	span	raw	cropped	asd
0	(0, 16)	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
1	(16, 32)	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
2	(32, 48)	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
3	(48, 64)	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>

Multi-channel Summary#

You can map custom functions (like calculating band RMS) using map().

def calc_rms(fs):
    return np.sqrt(np.sum(fs.value**2))  # Simplified band-integrated amplitude proxy for comparing segment-to-segment loudness.

st_asd.map("asd", calc_rms, out_col="band_rms", inplace=True)
st_asd.display()

	span	band_rms	raw	cropped	asd
0	(0, 16)	1.410130	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
1	(16, 32)	1.391752	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
2	(32, 48)	1.404694	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
3	(48, 64)	1.458929	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>