Note
このページは Jupyter Notebook から生成されました。 ノートブックをダウンロード (.ipynb)
[1]:
# Skipped in CI: Colab/bootstrap dependency install cell.
信号抽出: 色付き雑音からの微弱信号回収
このケーススタディは legacy の weak-signal notebook を公開 gallery 向けに整理したものです。物理的に見たいのは、生波形では見えない狭帯域信号をどうやって回収するかです。
ここでは弱い単色線を注入し、ASD で候補周波数を見つけ、band-pass で取り出してから簡単な正弦波モデルを当てます。
関連 API: Time Series API, Signal API
関連理論: 検証済みアルゴリズム, 数値安定性
[2]:
import matplotlib.pyplot as plt
import numpy as np
from gwexpy import TimeSeries
duration = 32
sample_rate = 4096
signal_freq = 123.4
signal_amp = 0.5
noise_std = 5.0
t = np.linspace(0, duration, int(duration * sample_rate), endpoint=False)
# The injected tone is far below the broadband noise floor, which is why it disappears in the raw waveform.
noise = np.random.normal(0, noise_std, size=len(t))
clean_signal = signal_amp * np.sin(2 * np.pi * signal_freq * t)
ts = TimeSeries(noise + clean_signal, t0=0, sample_rate=sample_rate, name="Noisy Data", unit="V")
1. 候補帯域を見つける
ASD の平均化は推定分散を減らすので、時間波形では見えない狭帯域線が統計的に見えるようになります。
[3]:
plot = ts.plot()
plot.gca().set_xlim(0, 0.1)
plot.gca().set_title("Raw Time Series (Zoomed)")
plt.show()
asd = ts.asd(fftlength=4, method="welch")
plot = asd.plot()
plot.gca().set_xlim(10, 1000)
plot.gca().set_yscale("log")
plot.gca().set_title("ASD")
plt.show()
2. フィルタしてフィットする
band-pass は注入線を残しつつ不要な広帯域成分を落とすためのものです。帯域を狭くしすぎると波形を歪め、広くしすぎると雑音を戻してしまいます。
関連 API: Time Series API, Fitting API
関連理論: 検証済みアルゴリズム
[4]:
# Keep the band around the detected line so the fit sees the signal-dominated portion of the data.
filtered_ts = ts.bandpass(110, 130).crop(1, 1.1)
plot = filtered_ts.plot()
plot.gca().set_title("Band-passed TimeSeries")
plt.show()
def sine_model(t, amp, freq, phase, amp2, freq2, phase2):
return amp * np.sin(2 * np.pi * freq * t + phase) * (1 + amp2 * np.sin(2 * np.pi * freq2 * t + phase2))
p0 = {"amp": 0.5, "freq": 123, "phase": 0, "amp2": 0.05, "freq2": 10, "phase2": 0}
limits = {"freq": (110, 130), "amp": (0.1, 5), "amp2": (0, 0.2), "freq2": (0, 50)}
result = filtered_ts.fit(sine_model, p0=p0, limits=limits, sigma=0.01)
print(result)
result.plot()
plt.show()
┌─────────────────────────────────────────────────────────────────────────┐
│ Migrad │
├──────────────────────────────────┬──────────────────────────────────────┤
│ FCN = 3.241e+05 (χ²/ndof = 802.2)│ Nfcn = 420 │
│ EDM = 3.68e-07 (Goal: 0.0002) │ │
├──────────────────────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Below EDM threshold (goal x 10) │
├──────────────────────────────────┼──────────────────────────────────────┤
│ SOME parameters at limit │ Below call limit │
├──────────────────────────────────┼──────────────────────────────────────┤
│ Hesse ok │ Covariance accurate │
└──────────────────────────────────┴──────────────────────────────────────┘
┌───┬────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │
├───┼────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ amp │ 224.7e-3 │ 0.7e-3 │ │ │ 0.1 │ 5 │ │
│ 1 │ freq │ 120.966 │ 0.017 │ │ │ 110 │ 130 │ │
│ 2 │ phase │ 10.83 │ 0.11 │ │ │ │ │ │
│ 3 │ amp2 │200.000e-3 │ 0.004e-3 │ │ │ 0 │ 0.2 │ │
│ 4 │ freq2 │ 9.35 │ 0.04 │ │ │ 0 │ 50 │ │
│ 5 │ phase2 │ 5.19 │ 0.24 │ │ │ │ │ │
└───┴────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
┌────────┬─────────────────────────────────────────────────────────────────────────────────────┐
│ │ amp freq phase amp2 freq2 phase2 │
├────────┼─────────────────────────────────────────────────────────────────────────────────────┤
│ amp │ 5.04e-07 -0.2e-6 1.4e-6 -161.87e-18 -5.1e-6 33.8e-6 │
│ freq │ -0.2e-6 0.000292 -1.86e-3 -41.40e-18 -0.02e-3 0.14e-3 │
│ phase │ 1.4e-6 -1.86e-3 0.0119 274.33e-18 0.0001 -0.001 │
│ amp2 │ -161.87e-18 -41.40e-18 274.33e-18 2.93e-19 1.74137e-15 -11.43654e-15 │
│ freq2 │ -5.1e-6 -0.02e-3 0.0001 1.74137e-15 0.00132 -0.0086 │
│ phase2 │ 33.8e-6 0.14e-3 -0.001 -11.43654e-15 -0.0086 0.0566 │
└────────┴─────────────────────────────────────────────────────────────────────────────────────┘
