Note

このページは Jupyter Notebook から生成されました。 ノートブックをダウンロード (.ipynb)

# Skipped in CI: Colab/bootstrap dependency install cell.

ASD 解析: パイプライン

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

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

import warnings

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)

import warnings

with warnings.catch_warnings():
    warnings.simplefilter('ignore')

    import numpy as np
    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

切り出し (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.425641	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
1	(16, 32)	1.395076	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
2	(32, 48)	1.382460	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>
3	(48, 64)	1.437775	<timeseries: 1024 samples>	<timeseries: 1024 samples>	<frequencyseries: 65 bins>