開発版ドキュメント · 0.2.3 30c2f8ba · 入門例の検証対象 0.2.3 · 版情報 · 既知の制限

キャリブレーション・パイプライン：Counts → Strain 変換#

重力波検出器はデータを デジタルカウント（生の ADC 出力）として記録します。天体物理学的な解析を行う前に、較正伝達関数を適用して、データを物理的なひずみ \(h(f)\) または変位 \(x(f)\) に変換する必要があります。

KAGRA の較正パイプラインは 2 つのステップから成ります：

  1. 結合係数の測定 — フォトンキャリブレータ（PCal）を介して既知の信号（スイープまたはステップ正弦波）を注入し、ResponseFunctionAnalysis を用いて検出器応答を測定します。

  2. 較正の適用 — 生スペクトルを結合関数で割り、較正済みのひずみ換算 ASD を生成します。

このチュートリアルで学ぶこと：

  1. 現実的なノイズを含む合成未較正データ（カウント）の生成

  2. ステップ正弦波注入から結合係数を抽出するための ResponseFunctionAnalysis の利用

  3. 周波数依存の較正モデルのフィッティング

  4. 較正の適用とひずみ較正済み ASD の生成

関連チュートリアル：

  • case_dttxml_calibration.ipynb — DTT XML ファイルからの測定済み伝達関数の読み込み

  • case_transfer_function.ipynb — 時系列から直接伝達関数を推定

準備#

import warnings


import matplotlib.pyplot as plt
import numpy as np

from gwexpy.analysis.response import ResponseFunctionAnalysis
from gwexpy.fitting import fit_series
from gwexpy.frequencyseries import FrequencySeries
from gwexpy.timeseries import TimeSeries

1. 合成未キャリブレーションデータの生成#

DARM（差動アーム長）信号をカウント単位でシミュレートします。 システムは以下のノイズ成分を持ちます：

  • 光学キャビティポール 50 Hz（単極低域通過特性）

  • 地震ノイズ 10 Hz 以下で 1/f^4 で上昇

  • ショットノイズ 100 Hz 以上でフラット

  • 電源周波数干渉 60 Hz および 120 Hz

fs_hz = 4096           # sample rate [Hz]
T     = 300.0          # total duration [s]
t0    = 1_300_000_000  # GPS epoch
N     = int(T * fs_hz)
t     = np.arange(N) / fs_hz
rng   = np.random.default_rng(42)

# Calibration parameters (known "true" values)
COUNTS_PER_METER = 4.5e10   # ADC counts per meter of DARM
F_POLE_HZ = 50.0             # optical cavity pole [Hz]

def cal_tf(f):
    # Counts/m transfer function: C(f) = C0 / (1 + i*f/f_pole)
    return COUNTS_PER_METER / (1 + 1j * f / F_POLE_HZ)

# Noise model in counts
freqs_noise = np.fft.rfftfreq(N, 1.0 / fs_hz)[1:]

# Shot noise (white, above cavity pole → flat counts)
shot_counts   = COUNTS_PER_METER * 3e-21

# Seismic noise: rises as (10/f)^4 metres, converted to counts
seismic_m  = 1e-12 * (10.0 / np.maximum(freqs_noise, 0.1))**4
seismic_cnts = np.abs(cal_tf(freqs_noise)) * seismic_m

noise_psd = shot_counts**2 + seismic_cnts**2
noise_amp = np.sqrt(noise_psd * fs_hz / 2)

noise_fft = noise_amp * np.exp(1j * rng.uniform(0, 2*np.pi, size=len(freqs_noise)))
noise_fft = np.concatenate([[0.0], noise_fft])   # DC = 0
noise_counts = np.fft.irfft(noise_fft, n=N)

# Add power-line interference at 60 Hz and 120 Hz
noise_counts += 500 * np.sin(2*np.pi*60  * t)
noise_counts += 150 * np.sin(2*np.pi*120 * t)

# Wrap in TimeSeries
ts_raw = TimeSeries(noise_counts, t0=t0, sample_rate=fs_hz,
                    name="K1:LSC-DARM_OUT_DQ", unit="ct")

print(f"Signal: {T:.0f} s at {fs_hz} Hz  ({N:,} samples)")
print(f"RMS (counts): {np.sqrt(np.mean(ts_raw.value**2)):.2f}")

Signal: 300 s at 4096 Hz  (1,228,800 samples)
RMS (counts): 1411.87

2. ステップサイン注入キャンペーン#

実際の較正運転では、検出器がロックされた状態で、一連の離散的な周波数で PCal 信号を注入します。各周波数ステップは通常 30〜120 秒続きます。

ここでは 10〜1000 Hz にわたる 6 周波数の注入キャンペーンを合成します。

inj_freqs = np.array([10.0, 30.0, 100.0, 200.0, 500.0, 1000.0])
step_dur  = 40.0   # seconds per step
pcal_amp  = 3e-13  # PCal amplitude [m]  (typical photon calibrator level)

# Build injection time series
ts_inj = np.zeros(N)
t_start = 10.0   # injection starts at GPS + 10 s

for f0 in inj_freqs:
    i0 = int(t_start * fs_hz)
    i1 = int((t_start + step_dur) * fs_hz)
    # Convert metres to counts via calibration TF
    amp_counts = np.abs(cal_tf(f0)) * pcal_amp
    ts_inj[i0:i1] = amp_counts * np.sin(2*np.pi*f0 * t[i0:i1])
    t_start += step_dur + 5.0   # 5 s gap between steps

ts_witness = TimeSeries(ts_inj, t0=t0, sample_rate=fs_hz,
                        name="K1:CAL-PCAL_SWEPT_EXC", unit="ct")

# Combined DARM signal = noise + injection response
ts_darm = TimeSeries(ts_raw.value + ts_inj, t0=t0, sample_rate=fs_hz,
                     name="K1:LSC-DARM_OUT_DQ", unit="ct")

# Quick overview plot
fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(t[:int(280*fs_hz)], ts_darm.value[:int(280*fs_hz)],
        lw=0.3, alpha=0.7, color="steelblue")
ax.set_xlabel("Time [s]")
ax.set_ylabel("DARM [counts]")
ax.set_title("Raw DARM signal with stepped-sine PCal injections")
plt.tight_layout()
plt.show()
../../_images/9b0ab6a770af418aaa89dda0e3c6eed13e0bb0d0082f1920b6b1f9eadcd5213a.png

3. カップリング係数の抽出#

ResponseFunctionAnalysis.compute() は witness チャンネル内の注入ステップを自動的に検出し、各注入周波数における結合係数を推定します。

結合係数 \(C(f_i)\) は次のように定義されます：

\[C(f_i) = \frac{\text{ASD}_{\text{target}}(f_i)}{\text{ASD}_{\text{witness}}(f_i)}\]
rfa = ResponseFunctionAnalysis()

result = rfa.compute(
    witness    = ts_witness,
    target     = ts_darm,
    fftlength  = 4.0,      # 4 s FFT segments
    overlap    = 0.5,
    snr_threshold = 5.0,   # minimum SNR to accept a step
    min_duration  = 10.0,  # minimum step duration [s]
)

print("Injection frequencies detected:", result.injected_freqs, "Hz")
print("Coupling factors (counts/count):", np.abs(result.coupling_factors))

# Compare with theoretical values
cf_theory = np.abs(cal_tf(result.injected_freqs)) / 1.0  # witness already in counts
print("\nMeasured vs. theoretical coupling factors:")
for f, cf_meas, cf_th in zip(result.injected_freqs,
                              np.abs(result.coupling_factors),
                              np.abs(cal_tf(result.injected_freqs))):
    print(f"  {f:6.0f} Hz: measured={cf_meas:.3e}, theory={cf_th:.3e}, "
          f"ratio={cf_meas/cf_th:.4f}")

Injection frequencies detected: [  10.   30.  100.  200.  500. 1000.] Hz
Coupling factors (counts/count): [5.14959810e+07 1.00121682e+00 1.00009536e+00 1.00002108e+00
 9.99998128e-01 1.00000068e+00]

Measured vs. theoretical coupling factors:
      10 Hz: measured=5.150e+07, theory=4.413e+10, ratio=0.0012
      30 Hz: measured=1.001e+00, theory=3.859e+10, ratio=0.0000
     100 Hz: measured=1.000e+00, theory=2.012e+10, ratio=0.0000
     200 Hz: measured=1.000e+00, theory=1.091e+10, ratio=0.0000
     500 Hz: measured=1.000e+00, theory=4.478e+09, ratio=0.0000
    1000 Hz: measured=1.000e+00, theory=2.247e+09, ratio=0.0000

4. キャリブレーション伝達関数のフィッティング#

測定された結合係数に単極モデルをフィットします：

\[C(f) = \frac{C_0}{1 + i\,f / f_\text{pole}}\]

このフィットにより、滑らかで連続的な較正関数が得られます。

freqs_meas = result.injected_freqs
cf_meas    = np.abs(result.coupling_factors)

# Build FrequencySeries for fitting
cf_fs = FrequencySeries(cf_meas.astype(complex), frequencies=freqs_meas,
                         name="Coupling function", unit="ct/ct")

# Define calibration model: |C(f)| = C0 / sqrt(1 + (f/f_pole)^2)
def cal_model_mag(f, C0, f_pole):
    return C0 / np.sqrt(1 + (f / f_pole)**2)

result_fit = fit_series(
    cf_fs.real,            # fit the magnitude (real part of the measured CF)
    cal_model_mag,
    p0=[COUNTS_PER_METER, 60.0],
)
C0_fit, fpole_fit = float(list(result_fit.params.values())[0]), float(list(result_fit.params.values())[1])

print(f"Fitted gain : {C0_fit:.3e} ct/m  (true: {COUNTS_PER_METER:.3e})")
print(f"Fitted pole : {fpole_fit:.2f} Hz  (true: {F_POLE_HZ:.2f} Hz)")

Fitted gain : 7.488e+09 ct/m  (true: 4.500e+10)
Fitted pole : -0.06 Hz  (true: 50.00 Hz)

# Plot coupling function: measured vs. fitted model
f_plot = np.geomspace(5, 2000, 500)
cf_model = cal_model_mag(f_plot, C0_fit, fpole_fit)

fig, ax = plt.subplots(figsize=(9, 4))
ax.loglog(freqs_meas, cf_meas, "o", ms=8, color="steelblue",
          label="Measured (stepped-sine)")
ax.loglog(f_plot, cf_model, "-", lw=2, color="tomato",
          label=f"Fit: C₀={C0_fit:.2e}, f_pole={fpole_fit:.1f} Hz")
ax.loglog(f_plot, np.abs(cal_tf(f_plot)), "--", lw=1.5, color="seagreen",
          alpha=0.7, label="True model")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("Coupling factor [ct/m]")
ax.set_title("Calibration Transfer Function")
ax.legend()
ax.grid(True, which="both", alpha=0.4)
plt.tight_layout()
plt.show()
../../_images/680390d5a47351694aee1583e642bbb53b6bb262f959f1f24583f28455505b42.png

5. キャリブレーション適用 — Counts → メートル ASD#

較正関数が得られたら、それを生の DARM ASD に適用し、metres per root-Hz 単位の較正済み変位 ASD を生成します。

# --- Compute raw ASD ---
asd_raw = ts_darm.asd(fftlength=16.0, method="median")
freqs_asd = asd_raw.frequencies.value

# Calibration function on the ASD frequency grid
cal_mag = cal_model_mag(np.maximum(freqs_asd, 0.1), C0_fit, fpole_fit)

# Calibrated ASD [m/sqrt(Hz)] = raw ASD [ct/sqrt(Hz)] / |C(f)| [ct/m]
asd_cal_data = asd_raw.value / cal_mag
asd_cal = FrequencySeries(asd_cal_data, frequencies=freqs_asd,
                          name="DARM calibrated", unit="m/sqrt(Hz)")

# True sensitivity (noise model in metres)
shot_m    = shot_counts / COUNTS_PER_METER
seismic_f = np.maximum(freqs_asd[1:], 0.1)
seismic_asd_m = 1e-12 * (10.0 / seismic_f)**4

fig, ax = plt.subplots(figsize=(10, 5))
ax.loglog(freqs_asd[1:], asd_cal.value[1:], color="steelblue", lw=1.5,
          alpha=0.8, label="Calibrated DARM")
ax.loglog(freqs_asd[1:], seismic_asd_m, "--", color="darkorange", lw=1.5,
          label="Seismic noise model")
ax.axhline(shot_m, color="seagreen", ls="--", lw=1.5, label="Shot noise floor")

# Mark calibration injection frequencies
for f0 in inj_freqs:
    ax.axvline(f0, color="gray", ls=":", lw=0.8, alpha=0.6)

ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("Displacement ASD [m/√Hz]")
ax.set_title("KAGRA DARM: Calibrated Sensitivity")
ax.set_xlim(1, 2000)
ax.legend()
ax.grid(True, which="both", alpha=0.4)
plt.tight_layout()
plt.show()
../../_images/97d595645ec55fa8fb833534e07cbfd3659673c147a5cbb040ab13860535264a.png

まとめ#

ステップ

API

目的

ステップサイン注入

TimeSeries + 合成信号

既知信号の注入

ステップ検出 + カップリング

ResponseFunctionAnalysis.compute()

各周波数での C(f) 抽出

キャリブレーションモデルフィット

fit_series(cf, model, p0)

連続的な C(f) の推定

キャリブレーション適用

ASD / C(f)

カウント → メートル/歪み

gwexpy.numerics.constants から使用される物理定数：

  • SAFE_FLOOR — ASD 推定における log(0) を防止

  • EPS_PSD — Welch 推定器の分母を正則化

実際の KAGRA 較正では：

  • PCal 振幅は PCal パワーモニターによって独立に測定される

  • 較正伝達関数は光学的成分と電気的成分の両方を持つ

  • 時間変化するパラメータ（フィネス、入力パワー）は継続的なモニタリングを必要とする

  • DTT XML からの測定済み伝達関数の読み込みについては case_dttxml_calibration.ipynb を参照