最初の解析#
Python コードで 2 つの信号を生成し、時間波形と振幅スペクトル密度（ASD）の図を保存します。 Python やスペクトル解析の知識は前提にしていません。 先にインストールを済ませてください。 学習時間の目安は 20〜30 分です。 小さな合成データを手元の環境で数秒程度で処理することを想定していますが、環境によって異なります。
検出器のチャンネルやスペクトルを扱っている方は、Python と実験概念の対応から同じ例を実行できます。
Python ファイルを実行する#
Python スクリプトは、名前が .py で終わるテキストファイルです。
Python はその文を上から順に実行します。
quickstart.py をダウンロードして、新しい作業フォルダに保存してください。
そのフォルダでターミナルを開き、インストールガイドの環境を有効にして、次を実行します。
python quickstart.py
このコマンドはターミナルに入力します。
quickstart.py の中にある行は Python コードです。
実行が終わったら、同じフォルダの asd.png を画像ビューアで開いてください。
クイックスタートで期待される図を確認できます。
ノートブックでも Python を実行できます。 スクリプトをコードセルに貼り付けて実行してください。 後のセルは先のセルで作った名前を使うため、順番に実行します。 図はノートブックの作業ディレクトリに保存されます。
名前、関数、チャンネル#
ダウンロードしたスクリプトを、この説明と並べて読んでください。
from gwexpy.noise.wave import gaussian, sine
from gwexpy.timeseries import TimeSeriesDict
settings = dict(duration=16, sample_rate=512, t0=0, unit="V")
tone = sine(frequency=40, **settings)
channels = TimeSeriesDict(
{
"Sensor A": tone + gaussian(std=0.3, seed=10, **settings),
"Sensor B": tone + gaussian(std=0.8, seed=20, **settings),
}
)
spectra = channels.asd(fftlength=2, overlap=1, window="hann", method="welch")
plot = spectra.plot(xlim=(1, 256), ylabel=r"ASD [V/$\sqrt{\mathrm{Hz}}$]")
plot.gca().legend()
plot.savefig("asd.png")
from ... import ... で始まる行は、ライブラリの道具を名前で使えるようにします。
tone = sine(...) のような代入は、関数を呼び出し、返されたオブジェクトを tone という名前で保存します。
括弧の中は関数への入力で、frequency=40 は 40 Hz の正弦波を指定します。
変数には数値、文字列、測定値とメタデータを含むオブジェクトなどを保存できます。
settings は、名前と値を組にした辞書です。
ここでは長さ 16 秒、サンプルレート 512 Hz、開始時刻 0、単位 V を記録します。
**settings は、名前付きの各値を信号生成関数に渡します。
これにより、生成する全チャンネルの時間軸と単位が揃います。
sine(...) は繰り返す振動、gaussian(...) はランダムなノイズを生成します。
seed で乱数列を固定すると、再実行時も同じサンプルを生成できます。
+ 演算子は正弦波とノイズをサンプルごとに足し、模擬センサー信号を作ります。
TimeSeries は 1 チャンネルの値、サンプリング情報、開始時刻、単位を持ちます。
TimeSeriesDict は名前ごとにチャンネルをまとめ、channels["Sensor A"] で 1 つを選択できます。
channels.asd(...) は全チャンネルの ASD を計算し、スペクトルのコレクションを返します。
channels.asd のドットは、そのオブジェクトが持つ操作や属性を選択します。
時間に対する信号を描く#
次のコードを quickstart.py の末尾に追加し、python quickstart.py をもう一度実行します。
first_second = channels.copy().crop(0, 1)
time_plot = first_second.plot(ylabel="Voltage [V]")
time_plot.gca().legend()
time_plot.savefig("timeseries.png")
print(channels["Sensor A"].sample_rate)
print(channels["Sensor A"].unit)
timeseries.png を開いてください。
横軸は時間、縦軸は電圧です。
振動にノイズが重なり、Sensor B の変動が大きくなるはずです。
.crop(0, 1) は 0 秒以上 1 秒未満の区間を選びます。
.crop() は対象のコレクションを更新します。
ここでは先に .copy() で別のコレクションを作るため、channels に全データが残ります。
2 つの print(...) 文は、ターミナルにサンプルレートと単位を表示します。
512 Hz と V が表示されるはずです。
これらは図のラベルだけでなく、チャンネル自体が保持している値です。
ASD の軸を読む#
ASD は、周波数帯域幅の平方根あたりの変動振幅を表します。 横軸は周波数（Hz、1 秒あたりの周期数）、縦軸は V/√Hz です。 40 Hz 付近の高いピークは、両チャンネルにある振動を示します。 周囲の広帯域フロアはノイズによるものです。 Sensor B はノイズの標準偏差が大きいため、フロアが高くなります。
スクリプトのスペクトル軸は対数目盛です。 等しい距離が等しい比率を表します。 例えば 10 Hz から 100 Hz までと、1 Hz から 10 Hz までの距離は同じです。 時間波形の図は変動がいつ起きたかを示し、ASD はその周波数成分を示します。
ASD の推定には長さ 2 秒の区間（fftlength=2）、Hann 窓、Welch 平均、1 秒のオーバーラップを使います。
区間長から、周波数ビンの間隔は 1 / 2 = 0.5 Hz になります。
窓関数は区間境界によるスペクトル漏れを減らし、オーバーラップは隣り合う区間でサンプルを再利用します。
ASD のピークの高さはこれらの設定に依存するため、元の正弦波の V 単位の振幅として読まないでください。
Python と実験概念の対応#
検出器の解析で使う用語と、上の実行可能なスクリプトを対応付けます。 Python が初めての方は、ファイルの実行と名前と関数で必要な構文を確認できます。
実験での概念 |
このレッスンの Python 表現 |
|---|---|
サンプル値を持つ 1 チャンネル |
|
チャンネル一覧 |
|
サンプルレートと区間開始時刻 |
|
工学単位 |
各チャンネルに付ける |
時間選択 |
|
FFT 長とオーバーラップ |
|
ASD トレース |
|
図の保存 |
|
ここでは合成データの便宜的な原点として t0=0 を使います。
記録した GW データでは、通常 t0 は秒単位の GPS 時刻を表します。
切出し範囲はチャンネルの開始時刻と同じ座標を使います。
GPS 時刻 t0 に始まるチャンネルの最初の 1 秒は、.crop(t0, t0 + 1) で選択できます。
条件を 1 つずつ変える#
正弦波の周波数を 40 Hz から 70 Hz に変え、ファイルを再実行して asd.png を確認します。
共通のスペクトルのピークが 70 Hz に移るはずです。
40 Hz に戻し、Sensor B の std だけを 0.8 から 0.3 に変えます。
広帯域フロアは同程度になりますが、シードが独立しているためノイズの線は一致しません。
さらに読む#
Commissionerのワークフロー：保存チャンネルを読み、ASD とコヒーレンスを比較します。
TimeSeries の基礎：時系列の操作を練習します。
はじめに：別の学習経路を選びます。