• KAGRA
• Subgroups
• PEM
• PythonMemoJP
• gwpy_TimeSeries

gwpy.timeseries.TimeSeries

TimeSeries とは？

gwpy.timeseries.TimeSeries は、GWpyの中で最も頻繁に使うクラスです。 単純な「数値の配列 (Numpy Array)」に、時間情報（開始時刻、サンプリングレート） や 単位 などの付加情報をセットにしたものです。

• データ本体: 電圧値やひずみなどの時系列データ

• t0: 開始時刻 (GPS時刻)

• dt: サンプリング間隔 (1/sample_rate)

• unit: データの単位 (V, m, strainなど)

• name: チャンネル名

これらがひとまとまりになっているため、「横軸を作成する手間」や「グラフのラベル付け」をGWpyが自動でやってくれます。

データの取得 (Fetch)

ローカルにあるファイル(gwf)を読むのではなく、サーバーから直接データを取得する方法です。

1. 公開データ (Open Data) を取得する場合

LIGO/Virgo/KAGRAの公開データサーバー (GWOSC) からデータを取得します。VPN不要で誰でも動かせます。

from gwpy.timeseries import TimeSeries

# LIGO Hanford (H1) の GW170817 前後のデータを32秒間取得
data = TimeSeries.fetch_open_data('H1', 1187008882, 1187008882+32)
print(data)

2. KAGRAのサーバー (NDS) から取得する場合

フレームファイルの時系列信号を読むチュートリアル (要JGWdocアカウント)

データのプロット (Plot)

TimeSeries オブジェクトは .plot() メソッドを持っており、一瞬でグラフを描画できます。 横軸（時刻）や単位のラベルは自動で付与されます。

plot = data.plot()
plot.show()

グラフのカスタマイズ

タイトルや軸の範囲を変えたい場合は、matplotlib の流儀で操作します。

import matplotlib.pyplot as plt

plot = data.plot(title='My First Plot', ylabel='Strain amplitude')
ax = plot.gca()  # get current axes (軸オブジェクトを取得)

ax.set_xlim(1187008882, 1187008882+5)  # 横軸の範囲指定
plot.show()

基本的なデータ操作

取得したデータに対して、よく行う操作を紹介します。

切り出し (Crop)

長いデータの一部だけを使いたい場合に使います。

# 最初の10秒間だけ切り出す
crop_data = data.crop(start=1187008882, end=1187008892)

リサンプリング (Resample)

データが重すぎる場合や、低い周波数を見たい場合に、サンプリングレートを下げます。

# 16384 Hz -> 4096 Hz にダウンサンプリング
resampled_data = data.resample(4096)

フィルタリング (Filter)

特定の周波数帯域を取り出します。

# 50Hz ハイパスフィルタ (50Hz以下の低周波ノイズをカット)
hp_data = data.highpass(50)

# 50Hz ~ 250Hz バンドパスフィルタ
bp_data = data.bandpass(50, 250)

# 60Hz ノッチフィルタ (商用電源ノイズをカット)
notch_data = data.notch(60)

白色化 (Whiten)

重力波データ解析で重要な処理です。感度の悪い（ノイズの大きい）周波数帯の重みを下げ、平坦化します。 微小な信号を見つけるためによく使われます。

white_data = data.whiten(fftlength=4, overlap=2)

周波数領域への変換 (FFT / スペクトル解析)

時系列データ (TimeSeries) から、周波数領域のデータ (FrequencySeries) を計算する方法です。 単なる FFT (numpy.fft) ではなく、データセットを分割して平均化する Welch法 がデフォルトで使われるため、安定したスペクトル評価が可能です。

ASD と PSD (単一チャンネル解析)

データのノイズレベルを評価する基本です。

• ASD (振幅スペクトル密度): 単位は 1/√Hz (ルートヘルツ分の1)。感度曲線のプロットには通常こちらを使います。

• PSD (パワースペクトル密度): 単位は 1/Hz (ヘルツ分の1)。ASDの2乗です。

gwpy.timeseries.TimeSeries.asd メソッドを使用します。

# ASDの計算
# fftlength: FFTを行う窓の長さ(秒)。周波数分解能 df = 1/fftlength になります。
# overlap: オーバーラップする長さ(秒)。通常は fftlength / 2 に設定します。
asd = data.asd(fftlength=4, overlap=2, method='welch')

# プロット
plot = asd.plot()
ax = plot.gca()
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel('Amplitude [1/sqrt(Hz)]')
plot.show()

CSD と コヒーレンス (2チャンネル相関解析)

2つのチャンネル間の相関を調べます。ノイズハンティング（環境モニターと重力波チャンネルの相関調査など）で頻出します。 ここでは、data1 と data2 という2つの時系列データが用意されていると仮定します。

CSD (クロススペクトル密度)

gwpy.timeseries.TimeSeries.csd 2つの信号の共通成分の強さを表します。

csd = data1.csd(data2, fftlength=4, overlap=2)

Coherence (コヒーレンス)

gwpy.timeseries.TimeSeries.coherence 2つの信号の相関の強さを 0 から 1 で正規化したものです。1に近いほど相関が強く、線形な関係にあります。

# コヒーレンスの計算
coh = data1.coherence(data2, fftlength=4, overlap=2)

# プロット (コヒーレンスは線形軸で見ることが多いです)
plot = coh.plot()
ax = plot.gca()
ax.set_xscale('log')
ax.set_ylabel('Coherence')
ax.set_ylim(0, 1)  # 0~1の範囲に固定
plot.show()

伝達関数 (Transfer Function)

gwpy.timeseries.TimeSeries.transfer_function ある信号が、システムを通してどのように別の信号へ伝わったか（ゲインと位相）を計算します。 数式的には H(f) = CSD_12 / PSD_11 で計算されます。

戻り値は複素数の FrequencySeries なので、ボード線図（ゲイン・位相）を描くには絶対値と偏角を取り出す必要があります。

# 伝達関数の計算 (data1 から data2 への伝達関数)
tf = data1.transfer_function(data2, fftlength=4, overlap=2)

# --- ボード線図のプロット ---
import numpy as np
import matplotlib.pyplot as plt

# 振幅 (Gain) と 位相 (Phase) の計算
mag = np.abs(tf.value)
phase = np.angle(tf.value, deg=True)
freq = tf.frequencies.value

# 描画
fig, (ax_mag, ax_phase) = plt.subplots(2, 1, sharex=True, figsize=(10, 8))

# Gain
ax_mag.loglog(freq, mag)
ax_mag.set_ylabel('Magnitude (Gain)')
ax_mag.grid(True, which='both', linestyle='--', alpha=0.5)

# Phase
ax_phase.semilogx(freq, phase)
ax_phase.set_ylabel('Phase [deg]')
ax_phase.set_xlabel('Frequency [Hz]')
ax_phase.set_ylim(-180, 180)
ax_phase.set_yticks(np.arange(-180, 181, 90))
ax_phase.grid(True, which='both', linestyle='--', alpha=0.5)

plt.show()

スペクトログラム (q_transform)

時系列データの色地図（横軸：時間、縦軸：周波数、色：強さ）を描きます。 グリッチノイズや重力波信号を視覚的に探すのに最強のツールです。

# Q-transform (Q値変換) の計算
q_spec = data.q_transform(frange=(20, 500), qrange=(4, 64))

# プロット
plot = q_spec.plot()
ax = plot.gca()
ax.set_yscale('log')  # 縦軸を対数に
plot.colorbar(label="Normalized energy")
plot.show()

属性へのアクセス (Attributes)

データの中身の数値そのものや、時間情報を取り出したい場合です。

例: Numpy配列としての利用

import numpy as np

# 最大値を探す
max_val = np.max(data.value)
print(f"最大値: {max_val}")