科学技術計算の Python から GWexpy へ#
NumPy 配列にサンプリング、開始時刻、単位を持たせ、チャンネルの辞書をコレクションのメソッドでまとめて解析します。 前提知識は NumPy 配列、Python の辞書、基本的な描画です。 GWpy の知識は必要ありません。 先にインストールを済ませてください。 例は標準の依存パッケージと合成データを使います。 学習時間の目安は 10〜15 分です。 小さな例をノート PC で数秒程度で実行することを想定していますが、環境によって異なります。
以下の Python ブロックをノートブックで順番に実行するか、1 つの .py ファイルにまとめて Python で実行します。
図は現在の作業ディレクトリに保存されます。
Before：配列と別々のメタデータ#
NumPy 配列はサンプル値を保存します。 スペクトルの計算にはサンプルレート、軸ラベルや時間選択には単位と時刻原点も必要です。 ここでは、それぞれを別の変数に保存しています。
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import welch
fs = 512
t0 = 1400000000
unit = "V"
time = np.arange(16 * fs) / fs
rng = np.random.default_rng(10)
values = np.sin(2 * np.pi * 40 * time) + rng.normal(0, 0.3, len(time))
frequency, psd = welch(
values, fs=fs, nperseg=2 * fs, noverlap=fs, window="hann",
detrend="constant", scaling="density", average="mean",
)
fig, ax = plt.subplots()
ax.loglog(frequency[1:], np.sqrt(psd[1:]))
ax.set(xlabel="Frequency [Hz]", ylabel=r"ASD [V/$\sqrt{\mathrm{Hz}}$]")
fig.savefig("numpy-asd.png")
plt.close(fig)
welch はパワースペクトル密度を返し、その平方根が図に示す ASD です。
values 自体は単位を持たないため、図のラベルを明示しています。
After：TimeSeries#
同じサンプルを、その解釈に必要なメタデータとともに包みます。
from gwexpy.timeseries import TimeSeries, TimeSeriesDict
series = TimeSeries(
values, sample_rate=fs, t0=t0, unit=unit, name="Sensor A"
)
spectrum = series.asd(fftlength=2, overlap=1, window="hann", method="welch")
plot = spectrum.plot(xlim=(1, 256))
plot.savefig("gwexpy-asd.png")
print(series.sample_rate, series.t0, series.unit)
print(spectrum.df, spectrum.unit)
入力配列は同じです。
series は sample_rate、t0、unit を保持し、返されたスペクトルには周波数軸と ASD の単位が付いています。
表示される周波数間隔は 0.5 Hz、スペクトルの単位は V/√Hz です。
GWexpy の fftlength と overlap は秒単位の長さです。
上の SciPy の nperseg と noverlap は、それに対応するサンプル数を指定します。
t0 は最初のサンプルの GPS 時刻です。
例えば series.crop(t0 + 2, t0 + 6) は、その時刻の 2 秒後から 4 秒間を選択します。
右端の境界にあるサンプルは含みません。
Before：辞書とスペクトル計算のループ#
NumPy で複数チャンネルを扱う場合、配列を辞書に保存し、スペクトル計算を繰り返すことがあります。
rng_b = np.random.default_rng(20)
arrays = {
"Sensor A": values,
"Sensor B": np.sin(2 * np.pi * 40 * time)
+ rng_b.normal(0, 0.8, len(time)),
}
numpy_spectra = {}
for name, data in arrays.items():
f, power = welch(
data, fs=fs, nperseg=2 * fs, noverlap=fs, window="hann",
detrend="constant", scaling="density", average="mean",
)
numpy_spectra[name] = (f, np.sqrt(power))
After：TimeSeriesDict#
各チャンネルを作るときにメタデータを付け、スペクトル設定をコレクションに適用します。
channels = TimeSeriesDict({
name: TimeSeries(data, sample_rate=fs, t0=t0, unit=unit, name=name)
for name, data in arrays.items()
})
spectra = channels.asd(fftlength=2, overlap=1, window="hann", method="welch")
plot = spectra.plot(xlim=(1, 256))
plot.gca().legend()
plot.savefig("channels-asd.png")
first_channel = channels["Sensor A"]
plain_values = first_channel.value
channels["Sensor A"] は TimeSeries、spectra["Sensor A"] はその FrequencySeries を選択します。
ASD の図には共通の 40 Hz の線が現れ、Sensor B のノイズフロアが高くなるはずです。
この例の全チャンネルはサンプルレート、開始時刻、長さ、単位が共通です。
実際の配列を包む際は、各チャンネルの実際のメタデータを指定してください。
.value は NumPy を受け取る関数に渡せるサンプル配列を返します。
外部関数が時刻や単位を必要とする場合は、それらも別途渡してください。
他の科学技術計算オブジェクトへの変換は相互運用ガイドを参照してください。
コレクションから次の解析へ#
作成・メタデータ・時間区間の選択・描画は TimeSeries の基礎、時間軸を揃えたチャンネルの行列表現は TimeSeriesMatrix の基礎へ進んでください。 ファイルを使って設定を保存する解析は、Commissionerのワークフローで実行できます。