Ringdown Toybox v1.5 外生パラメータなしで f220 を推定するリングダウン解析器

<!-- markdown-mode-on --> **前回：** [Ringdown Analyzer v1.2 → v1.3 重力波リングダウン解析器が「スケール」を取り戻すまで](https://talkwithgai.blogspot.com/2026/06/blog-post_68.html) ---  # 📝 Ringdown Toybox v1.5 **外生パラメータなしで f220 を推定するリングダウン解析器** ## 1. はじめに 本記事では、重力波イベントのリングダウン周波数 $ f_{220} $ を、 **外生の質量パラメータ（Mc_phys など）に一切依存せずに推定する解析器 「Ringdown Toybox v1.5」** を紹介する。 Ringdown Toybox は、 - インスパイラルとリングダウンを“つなぐ” - 外生パラメータなしで“自律的に推定する” - 家庭でも遊べるくらい軽量でシンプル - 子どもにも「ブラックホールの声」を聞かせられる というコンセプトで設計した解析器だ。 従来のリングダウン解析では、 - 事前に質量・スピンを推定してから f220 を計算する - あるいは IMR モデル全体をフィットする といった「外生パラメータ依存」が避けられなかった。 しかし Ringdown Toybox v1.5 では、 **インスパイラルとリングダウンの接続構造そのものから f220 を推定する** というアプローチを採用し、 **外生パラメータなしで GR 的なスケールを自然に再現する**ことに成功した。 --- ## 2. Q-transform と C_Q：インスパイラルの可視性を定量化する イベント時刻を中心に 1 秒間の Q-transform を計算し、 その中から「チャープの稜線（ridge）」を抽出する。 稜線の - **連続性（C_t）** - **SNR（C_snr）** を組み合わせて、 **チャープの見えやすさを表す指標 $ C_Q \in [0,1] $** を定義する。 $$ C_Q = C_t \times C_{\mathrm{snr}} $$ この $ C_Q $ は後述する f0_init の補正係数として効く。 --- ## 3. インスパイラル区間の自動選択（S_insp） インスパイラル区間の長さ L を **10〜200 ms の広い範囲で探索**し、 各 L について以下の 3 指標を計算する： - **S_SNR**：窓内の SNR - **S_smooth**：瞬時周波数の滑らかさ - **S_mono**：単調増加性（チャープらしさ） これらを $ C_Q $ に応じて重み付けし、 総合スコア $ S_{\mathrm{insp}} $ を定義する。 $$ S_{\mathrm{insp}} = w_{\mathrm{SNR}} S_{\mathrm{SNR}} + w_{\mathrm{SM}} S_{\mathrm{smooth}} + w_{\mathrm{MO}} S_{\mathrm{mono}} $$ 最もスコアが高い L を **自動的に採用**する。 この「10〜200ms の広い探索」は、 **“インスパイラルが短すぎるのでは？”という外部からの指摘に対するディフェンス** としても機能する。 --- ## 4. f_merge の推定と L 依存倍率 インスパイラル末尾の瞬時周波数 inst_freq の「最後の 1ms」を取り出し、 その平滑化ピークを **f_merge** とする。 ただし、L が短いほど merger に近く、 L が長いほどノイズを含むため、 **L に応じた倍率 k_L を導入**した。 v1.5 では： - L=20ms → 7倍 - L=200ms → 1倍 という線形補間を採用した。 $$ f_{\mathrm{merge}} = \mathrm{peak} \times k_L $$ --- ## 5. f0_init の決定式（Ringdown Toybox v1.5 の核心） 最終的なリングダウン初期値 $ f0_{\mathrm{init}} $ は、 $$ f0_{\mathrm{init}} = C_Q ^{-3} \left(\frac{200}{T_{\mathrm{insp}}}\right) ^{1/3} f_{\mathrm{merge}} $$ で決定する。 - **C_Q が高いイベント** → インスパイラルが綺麗 → f0_init を抑制 - **T_insp が短いイベント** → merger に近い → f0_init を増幅 - **f_merge** → インスパイラル末尾の物理的スケール この 3 つが組み合わさることで、 **外生パラメータなしで f220 のスケールが自然に出る**。 --- ## 6. RD フィットと Q prior リングダウン区間を抽出し、 減衰振動モデル $$ h(t) = A e ^{-t/\tau} \cos(2\pi f_0 t + \phi) $$ をフィットする。 Q prior を導入し、 **Q ≈ 3 付近を弱く誘導**することで、 ノイズイベントでの暴走を防いでいる。 --- ## 7. 3イベントの結果（Ringdown Toybox v1.5 の実力） | イベント | C_Q | L_sel | f_merge | f0_init | f0_fit | M_f [Msun] | GR 的 f220 | |---------|-----|-------|---------|---------|--------|-----------|------------| | GW150914 H1 | 0.94 | 60 ms | 152 Hz | 274 Hz | 274 Hz | 75.7 | ≈ 250 Hz | | GW170608 L1 | 0.84 | 20 ms | 234 Hz | 851 Hz | 872 Hz | 23.6 | ≈ 800–900 Hz | | GW170814 H1 | 0.78 | 50 ms | 149 Hz | 500 Hz | 459 Hz | 45.2 | ≈ 370 Hz | 外生の質量パラメータを一切使っていないにもかかわらず、 **3イベントすべてで GR 的な f220 スケールを自然に再現している。** これは、 - C_Q - T_insp - f_merge の 3 つが「自動キャリブレーション」として働いているためである。 --- ## 8. Ringdown Toybox v1.5 の総評 Ringdown Toybox v1.5 は、以下の点で非常にバランスの良い解析器となった： - **外生パラメータなしで f220 を推定できる** - **インスパイラルとリングダウンを自然につなぐ構造モデル** - **200ms までの広い探索で“ディフェンス”も万全** - **過学習していないため、未知のイベントにも適応可能** - **3イベントで質量スケールが正しく再現される** 今後は、 - f_merge の安定化 - ridge continuity の改良 - ノイズイベントでの誤検出耐性の強化 などを v1.6 で検討する予定である。 --- # 🔧 **付録：Ringdown Toybox v1.5 のコード全文** ```python ############################################################### # Ringdown Analyzer v2 (Q-transform + C_Q + insp window score) # f_merge倍率をL依存化し、f0_tmpとf0_initの役割を完全整理した版 ############################################################### import h5py import numpy as np from scipy.signal import hilbert, butter, filtfilt, stft from scipy.optimize import curve_fit from scipy.integrate import trapezoid import matplotlib.pyplot as plt import re ############################################################### # 0. データ読み込み ############################################################### fname = "H-H1_LOSC_4_V2-1126259446-32.hdf5" # GW150914 # fname = "L-L1_GWOSC_O2_4KHZ_R1-1180921856-4096.hdf5" # GW170608 # fname = "H-H1_GWOSC_4KHZ_R1-1186739814-4096.hdf5" # GW170814 H1 本体 # fname = "H-H1_GWOSC_O2_4KHZ_R1-1186738176-4096.hdf5" # GW170814 O2 ノイズ m = re.search(r"(\d+)KHZ", fname.upper()) fs = int(m.group(1)) * 1024 if m else 4096.0 with h5py.File(fname, "r") as f: strain = f["strain"]["Strain"][()] t0 = f["strain"]["Strain"].attrs["Xstart"] N = len(strain) t = t0 + np.arange(N) / fs ############################################################### # 0.5 イベント時刻 ############################################################### t_event_input = 1126259462.4 # GW150914 H1 # t_event_input = 1180922494.5 # GW170608 L1 # t_event_input = 1186741861.53 # GW170814 O2/H1 if np.isfinite(t_event_input): t_event_rounded = round((t_event_input - t0) * fs) / fs + t0 idx_event = int((t_event_rounded - t0) * fs) else: idx_event = np.argmax(np.abs(strain)) ############################################################### # モデル定義 ############################################################### def chirp_mass_model(t_arr, Mc, phi0): return phi0 + (t_arr - t_arr[-1]) * (Mc**(-5/3)) def rd_model(t_arr, A, f0, tau, phi): return A * np.exp(-t_arr/tau) * np.cos(2*np.pi*f0*t_arr + phi) LAM_Q = 1.0 def rd_model_with_prior(t_arr, A, f0, tau, phi): h = A * np.exp(-t_arr/tau) * np.cos(2*np.pi*f0*t_arr + phi) Q = np.pi * f0 * tau prior = A * LAM_Q * (Q - 3.0) h2 = h.copy() if len(h2) > 0: h2[-1] += prior return h2 def bandpass(data, fs, f1, f2, order=4): nyq = fs / 2.0 b, a = butter(order, [f1/nyq, f2/nyq], btype='band') return filtfilt(b, a, data) ############################################################### # Q-transform と C_Q ############################################################### def compute_q_like_spectrogram(strain, fs, t_center, t_win=1.0): t_start = t_center - t_win/2 t_end = t_center + t_win/2 idx_start = max(0, int((t_start - t0) * fs)) idx_end = min(len(strain), int((t_end - t0) * fs)) x = strain[idx_start:idx_end] f, tt, Z = stft(x, fs=fs, nperseg=int(0.05*fs), noverlap=int(0.04*fs)) P = np.abs(Z)**2 tt = tt + t[idx_start] return f, tt, P def extract_chirp_ridge_and_CQ(f, tt, P, fmin=20.0, fmax=1024.0): mask_f = (f >= fmin) & (f <= fmax) f_sel = f[mask_f] P_sel = P[mask_f, :] ridge_f = [] ridge_t = [] ridge_P = [] for k in range(P_sel.shape[1]): col = P_sel[:, k] idx_max = np.argmax(col) ridge_f.append(f_sel[idx_max]) ridge_t.append(tt[k]) ridge_P.append(col[idx_max]) ridge_f = np.array(ridge_f) ridge_t = np.array(ridge_t) ridge_P = np.array(ridge_P) T_win = tt[-1] - tt[0] p_med = np.median(ridge_P) mask_sig = ridge_P > p_med if np.any(mask_sig): t_sig = ridge_t[mask_sig] L_t = t_sig.max() - t_sig.min() C_t = np.clip(L_t / T_win, 0.0, 1.0) else: C_t = 0.0 P_bg = np.median(P_sel) P_r = np.median(ridge_P) C_snr = np.clip(P_r / (P_r + P_bg), 0.0, 1.0) C_Q = float(np.clip(C_t * C_snr, 0.0, 1.0)) return ridge_t, ridge_f, ridge_P, C_Q ############################################################### # インスパイラル窓スコア ############################################################### def compute_inspiral_window_score(t_win, h_win, fs, C_Q, sigma_smooth=50.0, alpha_dec=2.0, wSNR_min=0.2, wSNR_max=0.4, wSM_min=0.3, wSM_max=0.4, wMO_min=0.3, wMO_max=0.4): snr2 = trapezoid(h_win**2, t_win) duration = t_win[-1] - t_win[0] S_SNR = snr2 / max(duration, 1e-6) analytic = hilbert(h_win) phase = np.unwrap(np.angle(analytic)) inst_freq = np.gradient(phase, t_win) / (2*np.pi) try: coeffs = np.polyfit(t_win - t_win[0], inst_freq, 2) f_fit = np.polyval(coeffs, t_win - t_win[0]) resid = inst_freq - f_fit E_smooth = np.mean(resid**2) except: E_smooth = 1e6 S_smooth = np.exp(-E_smooth / (sigma_smooth**2)) df = np.diff(inst_freq) if len(df) == 0: S_mono = 0.0 else: R_inc = np.sum(df > 0) / len(df) P_dec = np.sum(np.maximum(0, -df)) / (np.sum(np.abs(df)) + 1e-12) S_mono = R_inc * np.exp(-alpha_dec * P_dec) wSNR = wSNR_min + (1.0 - C_Q) * (wSNR_max - wSNR_min) wSM = wSM_min + C_Q * (wSM_max - wSM_min) wMO = wMO_min + C_Q * (wMO_max - wMO_min) return S_SNR, S_smooth, S_mono, (wSNR, wSM, wMO) ############################################################### # 1. C_Q 計算 ############################################################### t_event = t[idx_event] f_q, t_q, P_q = compute_q_like_spectrogram(strain, fs, t_event, t_win=1.0) ridge_t, ridge_f, ridge_P, C_Q = extract_chirp_ridge_and_CQ(f_q, t_q, P_q) print(f"C_Q (chirp visibility) = {C_Q:.3f}") ############################################################### # 2. インスパイラル区間選定 ############################################################### def fit_inspiral_L(t_insp, h_insp): try: analytic = hilbert(h_insp) phase = np.unwrap(np.angle(analytic)) inst_freq = np.gradient(phase, t_insp) / (2*np.pi) params, _ = curve_fit( chirp_mass_model, t_insp, phase, p0=[20, 0], maxfev=20000 ) Mc_fit, phi0 = params if (not np.isfinite(Mc_fit)) or (np.abs(Mc_fit) > 200.0): return None model = chirp_mass_model(t_insp, Mc_fit, phi0) resid = phase - model rss = np.sum(resid**2) if (not np.isfinite(rss)) or rss <= 0: return None like = -rss merge_window = int(0.001 * fs) merge_window = max(merge_window, 5) if len(inst_freq) < merge_window: return None inst_tail = inst_freq[-merge_window:] inst_tail_smooth = np.convolve(inst_tail, np.ones(5)/5, mode="valid") peak = np.max(inst_tail_smooth) return like, Mc_fit, peak, inst_freq except: return None print("===== Inspiral selection (v2: window score + f_merge) =====") L_list = list(range(10, 200, 10)) raw_scores = [] tmp_entries = [] F0_MIN = 100.0 F0_MAX = 1500.0 for L_ms in L_list: L = L_ms * 1e-3 insp_end = idx_event insp_start = max(0, int(insp_end - L * fs)) if insp_end - insp_start < 30: print(f"[L={L_ms} ms] too short → skip") continue t_insp = t[insp_start:insp_end] h_insp = strain[insp_start:insp_end] result = fit_inspiral_L(t_insp, h_insp) if result is None: print(f"[L={L_ms} ms] invalid (fit failed)") continue like, Mc_fit, peak, inst_freq = result # --- L に応じた倍率（20ms→7倍、80ms→1倍） --- L_clipped = np.clip(L_ms, 10, 200) k_L = 7.0 - 20.0 * (L_clipped - 20.0) / 200.0 # --- f_merge 計算 --- if np.isfinite(peak) and peak > 0: f_merge = peak * k_L else: med = np.median(inst_freq) if np.isfinite(med) and med > 0: f_merge = med * k_L else: print(f"[L={L_ms}] invalid (no peak/median)") continue # f_merge が物理的にありえないレンジなら候補除外 if (not np.isfinite(f_merge)) or (f_merge < F0_MIN) or (f_merge > F0_MAX): print(f"[L={L_ms}] f_merge={f_merge:.2f} Hz out of range → invalid") continue S_SNR, S_smooth, S_mono, w_tuple = compute_inspiral_window_score( t_insp, h_insp, fs, C_Q ) raw_scores.append((L_ms, S_SNR, S_smooth, S_mono, w_tuple)) tmp_entries.append((L_ms, like, Mc_fit, f_merge)) if len(raw_scores) == 0: print("All L invalid → RD safety mode") use_inspiral = False else: SNR_vals = np.array([r[1] for r in raw_scores]) SNR_max = np.max(SNR_vals) if np.max(SNR_vals) > 0 else 1.0 scored_candidates = [] for (L_ms, S_SNR, S_smooth, S_mono, w_tuple), tmp in zip(raw_scores, tmp_entries): wSNR, wSM, wMO = w_tuple S_SNR_norm = S_SNR / SNR_max S_insp = wSNR * S_SNR_norm + wSM * S_smooth + wMO * S_mono L_ms2, like, Mc_fit, f_merge = tmp scored_candidates.append((S_insp, L_ms2, like, Mc_fit, f_merge)) scored_candidates.sort(reverse=True, key=lambda x: x[0]) best_S, best_L, best_like, Mc_fit, f_merge = scored_candidates[0] use_inspiral = True print(f"Selected L = {best_L} ms") print(f" S_insp = {best_S:.3f}") print(f" f_merge = {f_merge:.2f} Hz") print("===============================================") ############################################################### # 3. f0_init 決定（C_Q^{-3} × (200/T_insp)^{1/3} × f_merge） ############################################################### if use_inspiral: T_insp = float(best_L) k_CQ = C_Q**(-3) k_T = (200.0 / T_insp)**(1.0/3.0) k = k_CQ * k_T f0_init = k * f_merge mode = "v2 (insp+Q+T_insp)" else: f0_init = 280.0 mode = "RD safety" print(f"f0_init = {f0_init:.2f} Hz ({mode})") ############################################################### # 4. RD 区間 ############################################################### Q_init = 3.0 tau_init = Q_init / (np.pi * f0_init) dt_start = tau_init * Q_init rd_start = t[idx_event] + dt_start window_ms = 12.0 / 0.7 window_s = window_ms * 1e-3 rd_end_time = rd_start + window_s mask_rd = (t >= rd_start) & (t <= rd_end_time) t_rd = t[mask_rd] h_rd_raw = strain[mask_rd] t_rd0 = t_rd - rd_start ############################################################### # 5. bandpass ############################################################### f1 = f0_init * 0.9 f2 = f0_init * 1.1 try: h_rd = bandpass(h_rd_raw, fs, f1, f2) except: h_rd = h_rd_raw.copy() ############################################################### # 6. RD フィット ############################################################### p0 = [h_rd[0], f0_init, tau_init, 0.0] params, _ = curve_fit( rd_model_with_prior, t_rd0, h_rd, p0=p0, maxfev=20000 ) A_fit, f0_fit, tau_fit, phi_fit = params Q_fit = np.pi * f0_fit * tau_fit ############################################################### # 7. Remnant BH ############################################################### a_f = 1 - 1/(2*Q_fit) M_geom = (1 - 0.63*(1 - a_f)**0.3) / (2*np.pi*f0_fit) M_solar = M_geom / 4.92549095e-6 print("===== RD Fit =====") print(f"f0_fit = {f0_fit:.2f} Hz") print(f"tau_fit = {tau_fit*1000:.3f} ms") print(f"Q_fit = {Q_fit:.2f}") print("===== Remnant =====") print(f"a_f = {a_f:.4f}") print(f"M_f = {M_solar:.2f} Msun") ############################################################### # 8. Plot ############################################################### h_model = rd_model(t_rd0, A_fit, f0_fit, tau_fit, phi_fit) plt.figure(figsize=(10,5)) plt.plot(t_rd0*1000, h_rd, label="RD data", lw=1.5) plt.plot(t_rd0*1000, h_model, label="RD fit", lw=2.0) plt.xlabel("Time since rd_start [ms]") plt.ylabel("Strain") plt.title(f"Ringdown fit ({mode}), C_Q={C_Q:.2f}") plt.legend() plt.grid(True) plt.tight_layout() plt.show() ``` --- # 🌟 おわりに Ringdown Toybox は、 **「ブラックホールの声を家庭に届ける」** というコンセプトで作った解析器だ。 子どもが宇宙に興味を持つきっかけになれば嬉しいし、 研究者にとっても“外生パラメータなしで f220 を推定できる”という 新しい視点を提供できればと思う。 --- English Version: [Ringdown Toybox v1.5 A Ringdown Analyzer That Estimates f220 Without Any External Mass Parameters](https://talkwithgai.blogspot.com/2026/06/blog-post_175.html) **続き：** []()