= gwpy.timeseries.TimeSeries =
<<TableOfContents>>

== TimeSeries とは？ ==
[[https://gwpy.github.io/docs/stable/api/gwpy.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) から取得する場合 ===

[[https://gwdoc.icrr.u-tokyo.ac.jp/DocDB/0165/T2516508/001/read_gwf.html|フレームファイルの時系列信号を読むチュートリアル]] (要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 / スペクトル解析) ==
時系列データ ([[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries/|TimeSeries]]) から、周波数領域のデータ ([[https://gwpy.github.io/docs/stable/api/gwpy.frequencyseries.FrequencySeries/|FrequencySeries]]) を計算する方法です。
単なる FFT ({{{numpy.fft}}}) ではなく、データセットを分割して平均化する '''Welch法''' がデフォルトで使われるため、安定したスペクトル評価が可能です。

=== ASD と PSD (単一チャンネル解析) ===
データのノイズレベルを評価する基本です。
 * '''ASD (振幅スペクトル密度)''': 単位は 1/√Hz (ルートヘルツ分の1)。感度曲線のプロットには通常こちらを使います。
 * '''PSD (パワースペクトル密度)''': 単位は 1/Hz (ヘルツ分の1)。ASDの2乗です。

[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries/#gwpy.timeseries.TimeSeries.asd|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 (クロススペクトル密度) ====
[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries/#gwpy.timeseries.TimeSeries.csd|gwpy.timeseries.TimeSeries.csd]]
2つの信号の共通成分の強さを表します。

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

==== Coherence (コヒーレンス) ====
[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries/#gwpy.timeseries.TimeSeries.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) ===
[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries/#gwpy.timeseries.TimeSeries.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) ==
データの中身の数値そのものや、時間情報を取り出したい場合です。

|| '''属性名''' || '''型''' || '''説明''' ||
|| {{{data.value}}} || numpy.ndarray || データの中身（数値配列）。Numpyとして計算したい時に使う。 ||
|| {{{data.t0}}} || astropy.time.Time || 開始時刻。{{{data.t0.value}}} でGPSの数値(float)が取れる。 ||
|| {{{data.dt}}} || astropy.units.Quantity || サンプリング間隔 (例: 0.00024 s)。 ||
|| {{{data.sample_rate}}} || astropy.units.Quantity || サンプリング周波数 (例: 4096.0 Hz)。 ||
|| {{{data.times}}} || astropy.units.Quantity || 横軸の時刻配列。{{{data.times.value}}} で配列として取得。 ||

=== 例: Numpy配列としての利用 ===
{{{
import numpy as np

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