開発版ドキュメント · 0.2.3 30c2f8ba · 入門例の検証対象 0.2.3 · 版情報 · 既知の制限

ASD 解析: パイプライン#

このチュートリアルでは、SegmentTable を使用して、複数セグメントにわたる一括 ASD（振幅スペクトル密度）解析パイプラインを構築する方法を学びます。各セグメントの切り出し（crop）、ASD計算、可視化までを流れるように実行します。

ここで使う区間は gwpy.segments.Segment、各行の波形データは GWpy の TimeSeries です。gwexpy はそれらの GWpy 基盤クラスを SegmentTable に載せ、crop() / asd() で表単位の一括処理を行う部分を拡張しています。セグメント解析の基礎関係は SegmentTable: 基本を参照してください。

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) と ASD 計算#

crop() と asd() という便利な Sugar API を使って、全セグメントを一括処理できます。

# 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>

解析結果の集計#

map() を使うことで、独自関数（例：特定の周波数帯域の RMS 計算）を適用してテーブルに結果を統合できます。

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>