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

Commissioner：GUI から保存可能な解析へ#

使い慣れたコミッショニングの手順を Python で再現します。 チャンネルの読込、時間区間の選択、ASD とコヒーレンスの比較、解析条件と図の保存まで進めます。 DiagGUI、ndscope、Virgo dataDisplay の GUI でチャンネルやスペクトル設定を選んでいる方を対象とします。 以下のコマンドを実行する以外の Python 経験は前提にしていません。 構文が必要なら最初の解析で確認できます。

前提：GWexpy のインストール。 主な例には合成データと HDF5 サポートを含む標準の依存パッケージを使い、検出器への接続は必要ありません。 学習時間の目安は 20〜30 分です。 ノート PC で数秒程度の実行を目標としますが、環境によって異なります。

GUI の選択をコードに対応付ける#

GUI の操作や設定

スクリプトでの対応

チャンネルを選択する

TimeSeriesDict のキー

ndscope の記録を読み込む

TimeSeriesDict.read(path, format="hdf.ndscope")

時間区間を選択する

channels.copy().crop(start, end)。 記録した GW データには GPS 秒を使います

ASD の FFT 長とオーバーラップを選ぶ

channels.asd(fftlength=2, overlap=1, window="hann", method="welch")

参照チャンネルを選ぶ

sensor.coherence(reference, fftlength=2, overlap=1, window="hann")

トレースを図として保存する

plot.savefig("asd.png")

解析条件を記録する

チャンネル名、時刻、スペクトル設定を保存する JSON ファイル

Virgo dataDisplay は、読者が使い慣れた解析環境として挙げています。 このチュートリアルでは以下の ndscope と DiagGUI の形式を明示して使用します。 dataDisplay 専用の直接読込機能を示すものではありません。

ローカルで一連の解析を実行する#

commissioner.py を作業フォルダにダウンロードします。 ターミナルでインストールガイドの環境を有効にし、次を実行します。

python commissioner.py

スクリプトは commissioner-output/ に channels.hdf5、asd.png、coherence.png、analysis-parameters.json を作成します。 再実行すると、これらの練習用出力を置き換えます。 PNG は画像ビューア、JSON はテキストエディタで開いてください。

"""Save synthetic ndscope data, analyze a segment, and record the settings."""

import json
import platform
from importlib.metadata import version
from pathlib import Path

from gwexpy.noise.wave import gaussian, sine
from gwexpy.timeseries import TimeSeriesDict

# settings-begin
output = Path("commissioner-output")
output.mkdir(exist_ok=True)
parameters = {
    "sample_rate_hz": 512,
    "t0_gps_s": 1400000000,
    "duration_s": 32,
    "unit": "V",
    "tone_hz": 40,
    "noise_std_v": {"X1:REFERENCE": 0.3, "X1:SENSOR": 0.8},
    "noise_seeds": {"X1:REFERENCE": 10, "X1:SENSOR": 20},
    "crop_offset_s": [4, 28],
    "fftlength_s": 2,
    "overlap_s": 1,
    "window": "hann",
    "asd_method": "welch",
    "reference_channel": "X1:REFERENCE",
    "sensor_channel": "X1:SENSOR",
}
# settings-end

# data-begin
settings = dict(
    duration=parameters["duration_s"],
    sample_rate=parameters["sample_rate_hz"],
    t0=parameters["t0_gps_s"],
    unit=parameters["unit"],
)
tone = sine(frequency=parameters["tone_hz"], **settings)
channels = TimeSeriesDict(
    {
        name: tone
        + gaussian(std=parameters["noise_std_v"][name], seed=seed, **settings)
        for name, seed in parameters["noise_seeds"].items()
    }
)
data_path = output / "channels.hdf5"
channels.write(data_path, format="hdf.ndscope", overwrite=True)
loaded = TimeSeriesDict.read(data_path, format="hdf.ndscope")
print("Loaded channels:", list(loaded))
# data-end

# analysis-begin
t0 = parameters["t0_gps_s"]
start = t0 + parameters["crop_offset_s"][0]
end = t0 + parameters["crop_offset_s"][1]
segment = loaded.copy().crop(start, end)
spectral_settings = dict(
    fftlength=parameters["fftlength_s"],
    overlap=parameters["overlap_s"],
    window=parameters["window"],
)
spectra = segment.asd(method=parameters["asd_method"], **spectral_settings)
asd_plot = spectra.plot(xlim=(1, 256), ylabel=r"ASD [V/$\sqrt{\mathrm{Hz}}$]")
asd_plot.gca().legend()
asd_plot.savefig(output / "asd.png")

reference = segment[parameters["reference_channel"]]
sensor = segment[parameters["sensor_channel"]]
coherence = sensor.coherence(reference, **spectral_settings)
coherence_plot = coherence.plot(
    xlim=(1, 256), ylim=(0, 1), yscale="linear", ylabel="Magnitude-squared coherence"
)
coherence_plot.savefig(output / "coherence.png")
# analysis-end

# save-begin
parameters["crop_start_gps_s"] = start
parameters["crop_end_gps_s"] = end
parameters["input_file"] = str(data_path)
parameters["versions"] = {
    package: version(package)
    for package in ("gwexpy", "gwpy", "numpy", "scipy", "astropy", "h5py")
}
parameters["versions"]["python"] = platform.python_version()
(output / "analysis-parameters.json").write_text(
    json.dumps(parameters, indent=2) + "\n", encoding="utf-8"
)
print("Saved data, figures, and analysis-parameters.json in", output)
# save-end

保存データと時間選択を確認する#

スクリプトは、GPS 1400000000 を開始時刻とする長さ 32 秒、512 Hz サンプリングの電圧チャンネルを 2 つ作ります。 両方に 40 Hz の正弦波と、独立したシードのガウスノイズを含めます。 チャンネル名は合成データ用の例です。

channels.write(..., format="hdf.ndscope") は ndscope 形式の HDF5 ファイルを作ります。 TimeSeriesDict.read(..., format="hdf.ndscope") はチャンネル値とメタデータを読み込みます。 表示される一覧には X1:REFERENCE と X1:SENSOR が含まれます。 これらの公開 API は必要な I/O ハンドラを呼び出し時に読み込みます。

開始時刻から 4 秒後から 28 秒後まで、つまり GPS 1400000004 以上、GPS 1400000028 未満を切り出します。 両チャンネルに同じ 24 秒間を選択します。 TimeSeriesDict.crop() はコレクションを更新するため、スクリプトでは先にコピーしています。 実際の記録では、区間と参照を選ぶ前にチャンネル名、開始時刻、サンプルレート、単位を確認してください。

2 枚の図を読み取る#

ASD の図では、両チャンネルの 40 Hz 付近に線が現れ、X1:SENSOR の広帯域ノイズフロアが高くなります。 単位は V/√Hz です。 Hann 窓、2 秒の FFT 区間、1 秒のオーバーラップ、Welch 平均を使い、周波数ビンの間隔は 0.5 Hz になります。

2 枚目は二乗コヒーレンスです。 各周波数での線形な関連を 0 から 1 の無次元量で示します。 共通の 40 Hz の正弦波付近で高くなるはずです。 その周囲では独立ノイズにより小さく変動する推定値になり、有限回の平均では正確にゼロにはなりません。 高い値は共通のスペクトル成分を示しますが、それだけで物理的な結合の方向は判断できません。

analysis-parameters.json は入力ファイル、チャンネル選択、絶対時刻の切出し範囲、FFT 設定、合成データのシード、Python とパッケージのバージョンを記録します。 計算を再実行できるように、図とともに保存してください。

DiagGUI の時系列出力を読み込む#

この節では追加の依存パッケージ dttxml を同じ環境にインストールします。

python -m pip install dttxml

小さな合成サンプル commissioner.xml を作業フォルダにダウンロードします。 以下を同じ場所に read_diaggui.py として保存し、python read_diaggui.py を実行します。

from gwexpy.timeseries import TimeSeriesDict

diaggui = TimeSeriesDict.read(
    "commissioner.xml", format="xml.diaggui", products="TS", unit="V"
)
print(list(diaggui))
for name, channel in diaggui.items():
    print(name, channel.t0, channel.sample_rate, channel.unit)
plot = diaggui.plot()
plot.savefig("diaggui-timeseries.png")

products="TS" は保存された時系列データを選択します。 合成サンプルには TEST:SYNTHETIC_INPUT という名前で、4 Hz の電圧サンプルを 4 点含めています。 この時系列アダプタは XML から物理単位を復元しないため、例では unit="V" を明示しています。 この単位は、付属サンプルについて既知の校正上の前提です。 実際の出力には、その校正に基づく単位を指定してください。 DiagGUI は周波数領域の結果も保存できますが、それらは別のスペクトルとして、対応する周波数系列リーダーで読み込みます。 この例には時系列データを含む出力が必要です。 保存したスペクトル結果を読む場合の入口は I/O 形式を参照してください。

さらに読む#