カップリング関数解析チュートリアル

このチュートリアルでは gwexpy.analysis.CouplingFunctionAnalysis を使ったカップリング関数（CF）解析のワークフローを解説します。

カップリング関数とは

カップリング関数は、ウィットネス（環境センサー）のノイズがターゲット（感度チャンネル）にどれだけ結合するかを定量化します:

\[\mathrm{CF}(f) = \sqrt{\frac{P_{\rm target,inj}(f) - P_{\rm target,bkg}(f)}{P_{\rm witness,inj}(f) - P_{\rm witness,bkg}(f)}}\]

ここで \(P_{\rm inj}\) と \(P_{\rm bkg}\) はそれぞれ注入時・バックグラウンド時のパワースペクトル密度です。

典型的な用途: KAGRA/LIGO/Virgo の PEM（Physical Environment Monitor）インジェクション解析

レガシーコード移行メモ:
scipy で PSD を手動計算してカスタム閾値処理していた既存コードは、このワークフローに置き換えられます。
API_MAPPING.md の移行ガイドも参照してください。

[ ]:

import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u

from gwexpy.timeseries import TimeSeries, TimeSeriesDict
from gwexpy.analysis.coupling import (
    CouplingFunctionAnalysis,
    RatioThreshold,
    SigmaThreshold,
    PercentileThreshold,
)

1. 合成インジェクション／バックグラウンドデータの準備

シミュレーション設定:

ウィットネスチャンネル（環境センサー）: 地震計・音響センサー
ターゲットチャンネル 2 本: 結合強度の異なる感度計

インジェクションではウィットネスにトーン信号を加えます。ターゲットへの応答は結合強度に比例します。

[ ]:

rng = np.random.default_rng(0)

fs   = 512
T    = 32.0
n    = int(fs * T)
t    = np.arange(n) / fs
dt   = (1 / fs) * u.s

f_inj1 = 12.5   # Hz
f_inj2 = 25.0   # Hz
noise_level = 0.1

# バックグラウンドデータ（インジェクションなし）
wit_bkg  = rng.normal(0, noise_level, n)
tgt1_bkg = rng.normal(0, noise_level, n)
tgt2_bkg = rng.normal(0, noise_level, n)

# インジェクションデータ（ウィットネスにトーンを加える）
inj_signal = 5.0 * np.sin(2 * np.pi * f_inj1 * t) + 3.0 * np.sin(2 * np.pi * f_inj2 * t)
wit_inj  = inj_signal + rng.normal(0, noise_level, n)

# ターゲット 1: 強結合（CF ≈ 0.4）
cf_true1 = 0.4
tgt1_inj = cf_true1 * inj_signal + rng.normal(0, noise_level, n)

# ターゲット 2: 弱結合（CF ≈ 0.05）
cf_true2 = 0.05
tgt2_inj = cf_true2 * inj_signal + rng.normal(0, noise_level, n)

def make_ts(data, name, unit=u.ct):
    return TimeSeries(data, dt=dt, t0=0*u.s, unit=unit, name=name)

data_bkg = TimeSeriesDict({
    "ACC:FLOOR_X": make_ts(wit_bkg,  "ACC:FLOOR_X", u.m / u.s**2),
    "GW:DARM":     make_ts(tgt1_bkg, "GW:DARM",     u.m),
    "GW:AUX":      make_ts(tgt2_bkg, "GW:AUX",      u.m),
})

data_inj = TimeSeriesDict({
    "ACC:FLOOR_X": make_ts(wit_inj,  "ACC:FLOOR_X", u.m / u.s**2),
    "GW:DARM":     make_ts(tgt1_inj, "GW:DARM",     u.m),
    "GW:AUX":      make_ts(tgt2_inj, "GW:AUX",      u.m),
})

print("バックグラウンドチャンネル:", list(data_bkg.keys()))

2. カップリング関数解析の実行

CouplingFunctionAnalysis.compute() にインジェクション・バックグラウンドの TimeSeriesDict とウィットネスチャンネル名を渡します。

[ ]:

cfa = CouplingFunctionAnalysis()

results = cfa.compute(
    data_inj=data_inj,
    data_bkg=data_bkg,
    fftlength=4.0,
    overlap=0.5,
    witness="ACC:FLOOR_X",
    threshold_witness=RatioThreshold(ratio=10.0),
    threshold_target=RatioThreshold(ratio=2.0),
)

print("解析済みターゲット:", list(results.keys()))

3. 結果の確認

CouplingResult オブジェクトにカップリング関数・上限・診断 PSD が格納されています。

[ ]:

res_darm = results["GW:DARM"]

freqs = res_darm.cf.frequencies.value
idx1 = np.argmin(np.abs(freqs - f_inj1))
idx2 = np.argmin(np.abs(freqs - f_inj2))

print(f"有効 CF ビン数: {res_darm.valid_mask.sum()}")
print(f"f_inj1={f_inj1} Hz での CF: {res_darm.cf.value[idx1]:.3f}  (期待値 ~{cf_true1})")
print(f"f_inj2={f_inj2} Hz での CF: {res_darm.cf.value[idx2]:.3f}  (期待値 ~{cf_true1})")

res_darm.plot_cf(xlim=(5, 100));

[ ]:

# 3 パネル診断プロット: ウィットネス ASD | ターゲット ASD | CF
res_darm.plot(xlim=(5, 100));

[ ]:

# 弱結合チャンネル
res_aux = results["GW:AUX"]
print("GW:AUX の有効 CF ビン数:", res_aux.valid_mask.sum())
res_aux.plot_cf(xlim=(5, 100));

4. 閾値戦略の比較

gwexpy には 3 種類の組み込み閾値戦略があります:

戦略	説明	適したデータ
`RatioThreshold(ratio)`	\(P_{\rm inj} > r \cdot P_{\rm bkg}\)	高速スクリーニング・物理的な根拠がある場合
`SigmaThreshold(sigma)`	ガウス統計的有意性検定	定常・ガウス背景
`PercentileThreshold(pct, factor)`	経験的分布	非ガウス・グリッチが多い背景

[ ]:

strategies = {
    "Ratio×10/2": (RatioThreshold(10.0), RatioThreshold(2.0)),
    "Sigma 3σ":   (SigmaThreshold(3.0),  SigmaThreshold(2.0)),
}

fig, ax = plt.subplots(figsize=(8, 4))

for label, (thr_w, thr_t) in strategies.items():
    res = cfa.compute(
        data_inj=data_inj, data_bkg=data_bkg,
        fftlength=4.0, overlap=0.5,
        witness="ACC:FLOOR_X",
        threshold_witness=thr_w, threshold_target=thr_t,
    )
    cf = res["GW:DARM"].cf
    valid = ~np.isnan(cf.value)
    ax.plot(cf.frequencies.value[valid], cf.value[valid], marker=".", label=label)

ax.axhline(cf_true1, color="gray", linestyle="--", label=f"真の CF={cf_true1}")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(5, 100)
ax.set_xlabel("周波数 [Hz]")
ax.set_ylabel("カップリング関数")
ax.set_title("閾値戦略の比較")
ax.legend()
plt.tight_layout()

5. 上限値（Upper Limit）

ウィットネスが励起されているがターゲットに有意な過剰が見られない周波数帯では、
カップリング関数の上限値が自動計算されます。
ノイズフロア以下でも結合強度の制約を与えるために有用です。

[ ]:

cf_ul = res_darm.cf_ul
ul_valid = ~np.isnan(cf_ul.value)
print(f"上限値ビン数: {ul_valid.sum()}")

# plot_cf には CF と上限値が両方表示される
res_darm.plot_cf(xlim=(5, 100));

まとめ

from gwexpy.analysis.coupling import CouplingFunctionAnalysis, RatioThreshold

cfa = CouplingFunctionAnalysis()
results = cfa.compute(
    data_inj=data_inj,        # インジェクション TimeSeriesDict
    data_bkg=data_bkg,        # バックグラウンド TimeSeriesDict
    fftlength=4.0,
    witness="ACC:FLOOR_X",
    threshold_witness=RatioThreshold(25.0),
    threshold_target=RatioThreshold(4.0),
)

res = results["GW:DARM"]
res.plot()      # 3 パネル診断プロット
res.cf          # FrequencySeries: CF 値
res.cf_ul       # FrequencySeries: 上限値

関連チュートリアル: