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

ケーススタディ: Bruco と ICA によるノイズ削減#

このノートブックでは、実際の干渉計コミッショニングで使われる エンドツーエンドのノイズ削減ワークフロー を示します：

  1. Bruco — ブルートフォースコヒーレンス スキャンで最相関の補助チャンネルを特定

  2. ICA — 独立成分分析でノイズ源を分離

  3. ノイズ差し引き — DARM からノイズ寄与を除去し、ASD を比較

これは O4b コミッショニング（DARM 116 Hz ライン調査など）で使われたワークフローを再現しています。 物理的な狙いは、「どのセンサに線が見えるか」から一歩進んで、「どの潜在ノイズ源なら DARM から差し引いても他の構造を壊さないか」を判断することです。

前提条件：以下に慣れていること

準備#

import warnings

import warnings
with warnings.catch_warnings():

    # ruff: noqa: I001
    import matplotlib.pyplot as plt
    import numpy as np
    from astropy import units as u

    from gwexpy.analysis import Bruco
    from gwexpy.timeseries import TimeSeries, TimeSeriesDict, TimeSeriesMatrix

1. モックデータ生成#

DARM チャンネルに以下が混入するシナリオをシミュレートします：

  • 広帯域ガウスノイズ（バックグラウンド雑音フロア）

  • 環境センサーから漏れた 116 Hz ラインノイズ

4 つの PEM（Physical Environment Monitor）補助チャンネルは、それぞれ異なる割合で同じ 116 Hz ラインを含みます。 これは O4b コミッショニングの実際のケースを模倣しています。

rng = np.random.default_rng(42)

fs = 512        # sample rate [Hz]
T  = 64.0       # duration [s]
n  = int(fs * T)
t  = np.arange(n) / fs

FREQ_LINE = 116.0  # Hz — the line noise we want to remove

# ----- Independent sources -----
src_line  = np.sin(2 * np.pi * FREQ_LINE * t)   # coherent line noise
src_broad = rng.normal(0, 1, n)                  # broadband floor noise

# ----- DARM: broadband + line leakage -----
DARM_LINE_COUPLING = 0.5
darm_noise = src_broad + DARM_LINE_COUPLING * src_line

target = TimeSeries(
    darm_noise, dt=1 / fs, unit=u.dimensionless_unscaled,
    name="K1:CAL-CS_PROC_DARM_STRAIN_DBL_DQ", t0=0,
)

# ----- 4 auxiliary channels (different line content) -----
aux_configs = {
    "K1:PEM-ACC_PSL_TABLE_PSL1_Y": (0.9, 0.1),  # 90% line, 10% noise
    "K1:PEM-ACC_PSL_TABLE_PSL2_X": (0.7, 0.3),
    "K1:PEM-MIC_PSL_TABLE_PSL1_Z": (0.5, 0.5),
    "K1:PEM-MIC_PSL_TABLE_PSL2_Z": (0.2, 0.8),  # mostly noise
}

aux_dict = TimeSeriesDict({
    name: TimeSeries(
        a_line * src_line + a_noise * rng.normal(0, 1, n),
        dt=1 / fs, unit=u.dimensionless_unscaled, name=name, t0=0,
    )
    for name, (a_line, a_noise) in aux_configs.items()
        })

print(f"DARM   sample rate: {target.sample_rate}")
print(f"Aux channels: {list(aux_dict.keys())}")

DARM   sample rate: 512.0 Hz
Aux channels: ['K1:PEM-ACC_PSL_TABLE_PSL1_Y', 'K1:PEM-ACC_PSL_TABLE_PSL2_X', 'K1:PEM-MIC_PSL_TABLE_PSL1_Z', 'K1:PEM-MIC_PSL_TABLE_PSL2_Z']

# Visualize the line noise in the ASD before cleaning
fig, ax = plt.subplots(figsize=(10, 4))
asd = target.asd(fftlength=4, overlap=2)
ax.loglog(asd.frequencies.value, asd.value, label="DARM (original)", color="steelblue")

for name, ts in aux_dict.items():
    a = ts.asd(fftlength=4, overlap=2)
    ax.loglog(a.frequencies.value, a.value, alpha=0.5, lw=0.8, label=name.split(":")[-1])

ax.axvline(FREQ_LINE, color="red", ls="--", label=f"{FREQ_LINE} Hz line")
ax.set_xlim(1, fs / 2)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("ASD [1/√Hz]")
ax.set_title("ASD before cleaning — line noise visible at 116 Hz")
ax.legend(fontsize=7, ncol=2)
plt.tight_layout()
plt.show()
../../_images/03e4066bcc7f8ad99cc7dec328ec0164dc84a2a1d097e479f7e018c0c3f1b4ef.png

2. Bruco: 相関補助チャンネルの特定#

Bruco.compute() は全補助チャンネルをスキャンし、各周波数ビンで上位 N 件の最相関チャンネルを返します。 これがパイプラインの 第 1 ステージ です：DARM と同じノイズを観測しているセンサーを特定します。

よくある誤り：利用可能なすべての witness を ICA に投入すること。ここでは Bruco を物理的なフィルタとして用い、116 Hz の汚染をもっともらしく共有するチャンネルに分解を集中させます。

bruco = Bruco(target_channel=target.name, aux_channels=[])

result = bruco.compute(
    fftlength=4.0,
    overlap=2.0,
    target_data=target,
    aux_data=aux_dict,
    top_n=4,         # keep top-4 channels per frequency bin
)

print("Bruco scan complete.")
print(f"Result type: {type(result)}")

/home/runner/micromamba/envs/gwexpy/lib/python3.11/site-packages/gwpy/timeseries/core.py:263: UserWarning: xindex was given to TimeSeries(), x0 will be ignored
  new = super().__new__(

/home/runner/micromamba/envs/gwexpy/lib/python3.11/site-packages/gwpy/timeseries/core.py:263: UserWarning: xindex was given to TimeSeries(), x0 will be ignored
  new = super().__new__(

/home/runner/micromamba/envs/gwexpy/lib/python3.11/site-packages/gwpy/timeseries/core.py:263: UserWarning: xindex was given to TimeSeries(), x0 will be ignored
  new = super().__new__(

/home/runner/micromamba/envs/gwexpy/lib/python3.11/site-packages/gwpy/timeseries/core.py:263: UserWarning: xindex was given to TimeSeries(), x0 will be ignored
  new = super().__new__(

Bruco scan complete.
Result type: <class 'gwexpy.analysis.bruco.BrucoResult'>

# Show top channels near the 116 Hz line
df = result.to_dataframe(ranks=[0])
df_line = (
    df[df["frequency"].between(114, 118)]
    .sort_values("coherence", ascending=False)
    .dropna(subset=["channel"])
    .head(10)
    .reset_index(drop=True)
)
print("Top correlated channels near 116 Hz:")
df_line

Top correlated channels near 116 Hz:
frequency rank channel coherence projection
0 116.00 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.995049 0.696766
1 115.75 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.984059 0.334223
2 116.25 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.979161 0.333768
3 114.75 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.459202 0.029435
4 116.75 1 K1:PEM-ACC_PSL_TABLE_PSL2_X 0.403842 0.024020
5 115.00 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.330589 0.020492
6 117.00 1 K1:PEM-ACC_PSL_TABLE_PSL2_X 0.293278 0.020808
7 117.50 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.271092 0.016356
8 114.50 1 K1:PEM-ACC_PSL_TABLE_PSL1_Y 0.269835 0.020907
9 115.50 1 K1:PEM-MIC_PSL_TABLE_PSL1_Z 0.246493 0.014666 
# Coherence projection plot — shows which channels dominate at each frequency
result.plot_projection(coherence_threshold=0.3)
plt.xlim(1, fs / 2)
plt.axvline(FREQ_LINE, color="red", ls="--", lw=1)
plt.title("Bruco noise projection")
plt.show()
../../_images/598fee58393ce2d9360c26607d17b9796ac59f93a7d4a50927af1a888cc52824.png

3. ICA: ノイズ源の分離#

Bruco スキャンから、DARM と 116 Hz 付近で高いコヒーレンスを持つチャンネルが分かりました。 次にこれらのチャンネルを TimeSeriesMatrix に積み上げ、ICA を適用して 独立したノイズ源を分離します。

ICA モデルは混合行列 A を求めます：

X = S · A^T   (X: observed channels, S: independent sources)

この行列を使って DARM からノイズ成分の寄与を差し引けます。

注意すべき失敗モード：ICA 成分は物理的重要度で自動的に並べられるわけではなく、各成分の符号／スケールは任意です。どの成分が差し引きたい線を表すかを判断する前に、必ずスペクトルとチャンネル負荷を確認してください。

注 — ICA の収束：FastICA は一部のデータセットでは既定の反復回数内で収束しないことがあります。ConvergenceWarning が表示された場合は、ICA パラメータの max_iter または tol を増やすことを検討してください。完全には収束していなくても、結果が利用可能な場合があります。

二つ目のよくある失敗モードはランク落ちです。ほぼ同一の witness を ICA に投入すると、分解が不安定になり、差し引きの重みが信頼しにくくなります。情報量はあるが完全には冗長でない少数のチャンネルを選ぶようにしてください。

# Pick the top-2 channels from Bruco result
TOP_CHANNELS = df_line["channel"].dropna().unique()[:2].tolist()
print("Selected channels for ICA:", TOP_CHANNELS)

# Stack DARM + top channels into a TimeSeriesMatrix  (shape: n_ch × 1 × n_samples)
channels = [target] + [aux_dict[ch] for ch in TOP_CHANNELS]
data_3d = np.stack([ch.value for ch in channels], axis=0)[:, np.newaxis, :]

tsm = TimeSeriesMatrix(
    data_3d, dt=1 / fs, unit=u.dimensionless_unscaled, t0=0,
)
print(f"TimeSeriesMatrix shape: {tsm.shape}  (n_channels × 1 × n_samples)")

Selected channels for ICA: ['K1:PEM-ACC_PSL_TABLE_PSL1_Y', 'K1:PEM-ACC_PSL_TABLE_PSL2_X']
TimeSeriesMatrix shape: (3, 1, 32768)  (n_channels × 1 × n_samples)

import warnings
with warnings.catch_warnings():

    # Run ICA
    n_components = len(channels)
    ica_sources, ica_model = tsm.ica(n_components=n_components, return_model=True)

    sk = ica_model.sklearn_model
    print(f"ICA stopped after {sk.n_iter_} iterations")
    print(f"Mixing matrix A shape: {sk.mixing_.shape}")  # (n_channels, n_components)

    # Visualize ICA sources in frequency domain
    fig, axes = plt.subplots(1, n_components, figsize=(4 * n_components, 4), sharey=True)
    for k in range(n_components):
        src_ts = TimeSeries(
            ica_sources.value[k, 0, :], dt=1 / fs,
            unit=u.dimensionless_unscaled, t0=0,
        )
        asd_k = src_ts.asd(fftlength=4, overlap=2)
        axes[k].semilogy(asd_k.frequencies.value, asd_k.value)
        axes[k].axvline(FREQ_LINE, color="red", ls="--", lw=0.8)
        axes[k].set_title(f"ICA source #{k}")
        axes[k].set_xlim(1, fs / 2)
        axes[k].set_xlabel("Frequency [Hz]")

    axes[0].set_ylabel("ASD [arb/√Hz]")
    plt.suptitle("ICA independent components — one should peak at 116 Hz")
    plt.tight_layout()
    plt.show()

/home/runner/micromamba/envs/gwexpy/lib/python3.11/site-packages/sklearn/decomposition/_fastica.py:132: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations.
  warnings.warn(

ICA stopped after 200 iterations
Mixing matrix A shape: (3, 3)
../../_images/182a3ec3bbde5fdda7eadc0579b37e8414faf21fbbfca5a8366decb36e094590.png

4. ノイズ差し引きと ASD 比較#

116 Hz ラインを含む ICA 成分を特定し（ASD を確認して選択）、 混合行列を使って DARM からその寄与を差し引きます。

DARM_clean = DARM - Σ_k  A[0, k] · S_k(t)     (summed over noise components k)

よくある誤り：非ガウス的に見える ICA 成分をすべて差し引くこと。スペクトル内容と witness 負荷が対象とする結合メカニズムに一致する成分だけを差し引いてください。さもないと、線とともに無関係な検出器構造まで除去してしまう恐れがあります。

# Identify the noise component: the one with the largest ASD at FREQ_LINE
fftlength = 4.0
overlap   = 2.0
freqs = np.fft.rfftfreq(int(fftlength * fs), d=1 / fs)
line_bin = np.argmin(np.abs(freqs - FREQ_LINE))

A = sk.mixing_           # (n_channels, n_components)
X_raw = data_3d[:, 0, :].T   # (n_samples, n_channels)
S = sk.transform(X_raw)       # (n_samples, n_components)

# ASD of each component at the line frequency
asd_at_line = []
for k in range(n_components):
    s_ts = TimeSeries(S[:, k], dt=1 / fs, unit=u.dimensionless_unscaled, t0=0)
    asd_k = s_ts.asd(fftlength=fftlength, overlap=overlap)
    val = float(np.interp(FREQ_LINE, asd_k.frequencies.value, asd_k.value))
    asd_at_line.append(val)

noise_component = int(np.argmax(asd_at_line))
print(f"Dominant noise component: #{noise_component}  (ASD at {FREQ_LINE} Hz = {asd_at_line[noise_component]:.4f})")

# Subtract noise component from DARM (channel index 0)
darm_clean = X_raw[:, 0] - A[0, noise_component] * S[:, noise_component]
target_clean = TimeSeries(
    darm_clean, dt=1 / fs, unit=u.dimensionless_unscaled,
    name="DARM_cleaned", t0=0,
)

Dominant noise component: #0  (ASD at 116.0 Hz = 1.5730)

# Compare original vs cleaned ASD
asd_orig  = target.asd(fftlength=fftlength, overlap=overlap)
asd_clean = target_clean.asd(fftlength=fftlength, overlap=overlap)

fig, ax = plt.subplots(figsize=(10, 5))
ax.loglog(asd_orig.frequencies.value,  asd_orig.value,  label="DARM original",  color="steelblue", lw=1.2)
ax.loglog(asd_clean.frequencies.value, asd_clean.value, label="DARM cleaned",   color="orangered",  lw=1.2)
ax.axvline(FREQ_LINE, color="gray", ls="--", lw=0.8, label=f"{FREQ_LINE} Hz")

# Suppression factor at line frequency
ratio = float(np.interp(FREQ_LINE, asd_orig.frequencies.value, asd_orig.value)) /         float(np.interp(FREQ_LINE, asd_clean.frequencies.value, asd_clean.value))
ax.set_xlim(50, 200)
ax.set_ylim(1e-4, 10)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("ASD [1/√Hz]")
ax.set_title(f"Noise reduction at {FREQ_LINE} Hz: ×{ratio:.1f} suppression")
ax.legend()
plt.tight_layout()
plt.show()

print(f"Suppression factor at {FREQ_LINE} Hz: {ratio:.2f}×")
../../_images/af75c714e4f886e354b09e8adedc1d28ae34935717d67157402e8104d32330cb.png 
Suppression factor at 116.0 Hz: 2.74×

まとめ#

ステップ

ツール

返り値

1. コヒーレンス スキャン

Bruco.compute()

各周波数での最相関チャンネル

2. ノイズ源分離

TimeSeriesMatrix.ica()

独立成分 + 混合行列

3. ノイズ差し引き

混合行列の代数演算

クリーン DARM チャンネル

重要なポイント#

  • Bruco は数千のチャンネルを、重要な少数へと効率的に絞り込みます。

  • ICA はコヒーレンスを超えます。複数のセンサーが共通のノイズを共有している場合でも、独立なソースの寄与を分離します。

  • 混合行列 A は、時間領域フィルタリングなしで直接的な差し引きの式を与えます。

  • 安定した差し引きは、ICA ソルバーだけでなく、witness の選択と成分の同定に依存します。

次のステップ#

  • モックデータを実データ（TimeSeries.read() や NDS2）に置き換える。

  • より長いセグメントで解析し、Bruco の投影推定と照合する。

  • Spectrogram.normalize()（SNR スペクトログラム）と組み合わせて、ラインの時間変化を追跡する。