Wiener Filtering: Coherent Noise Subtraction#
A multichannel Wiener filter is useful when several witness sensors contain overlapping views of the same physical noise sources. Treating those witnesses independently can double-count the noise and over-subtract the target.
This compact public case study promotes the legacy MIMO notebook and focuses on the core matrix workflow.
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from gwexpy import FrequencySeriesDict, TimeSeries, TimeSeriesDict
from gwexpy.noise.asd import from_pygwinc
from gwexpy.noise.wave import from_asd
fs = 2048.0
duration = 64.0
t = np.arange(0, duration, 1 / fs)
n = len(t)
np.random.seed(123)
rng = np.random.default_rng(123)
# Two independent physical sources are mixed into two witnesses with different weights, creating crosstalk between sensors.
b_low, a_low = signal.butter(2, 5.0, fs=fs, btype="low")
s1 = signal.lfilter(b_low, a_low, np.random.normal(0, 10.0, n))
b_band, a_band = signal.iirpeak(120.0, 30.0, fs=fs)
s2 = signal.lfilter(b_band, a_band, np.random.normal(0, 2.0, n)) + np.random.normal(0, 0.1, n)
aux1_val = 1.0 * s1 + 0.4 * s2 + np.random.normal(0, 0.05, n)
aux2_val = 0.6 * s1 + 1.0 * s2 + np.random.normal(0, 0.05, n)
asd_main = from_pygwinc("aLIGO", fmin=5.0, fmax=fs / 2, df=1.0 / duration, quantity="strain")
tsd = TimeSeriesDict()
tsd["MAIN"] = from_asd(asd_main, duration, fs, t0=0, rng=rng).highpass(5.0)
tsd["AUX1"] = TimeSeries(aux1_val, sample_rate=fs, unit="V")
tsd["AUX2"] = TimeSeries(aux2_val, sample_rate=fs, unit="V")
# Add witness-coupled noise into MAIN so subtraction has a coherent target to remove.
tsd["MAIN"] += tsd["AUX1"] / tsd["AUX1"].unit * tsd["MAIN"].unit * 2e-22
tsd["MAIN"] += tsd["AUX2"] / tsd["AUX2"].unit * tsd["MAIN"].unit * 5e-22
1. Estimate the multichannel correlation matrices#
The key point is that the witnesses are not independent. The Cxx matrix captures witness-to-witness correlation, and Cyx captures how the target correlates with those witnesses.
This example uses the convention \(C_{ij}=\langle X_i^*X_j\rangle\), so the row filter mapping witnesses to MAIN is \(H=(C_{yx}C_{xx}^{-1})^*\). We use identical Hann windows, overlap, and mean averaging for every matrix entry. The numerical transfer estimates are assigned MAIN/witness units explicitly because the inherited spectral metadata does not infer these mixed-channel units.
The FFT and inverse FFT use their existing normalization; no extra sample-rate or factor-of-two correction is applied. This synthetic example has a nonsingular witness matrix. Check conditioning and choose a justified regularization before applying it to redundant measured witnesses.
aux_names = ["AUX1", "AUX2"]
aux_tsd = TimeSeriesDict({k: tsd[k] for k in aux_names})
# Use the same mean averaging for diagonal PSDs and off-diagonal CSDs.
spectral_options = dict(fftlength=8.0, overlap=4.0, window="hann", average="mean")
cxx = aux_tsd.csd_matrix(**spectral_options)
cyx = TimeSeriesDict({"MAIN": tsd["MAIN"]}).csd_matrix(
other=aux_tsd, **spectral_options
)
# Cij = <conj(X_i) X_j>; conjugate to map witnesses to MAIN.
H_lowres = (cyx @ cxx.inv()).conj()
# Attach the physical output/input units to the numerical transfer estimates.
for j, name in enumerate(aux_names):
H_lowres.meta[0, j].unit = tsd["MAIN"].unit / tsd[name].unit
assert np.all(np.isfinite(H_lowres[0, j].value))
H_lowres.abs().plot(xscale="log", yscale="log").suptitle("Estimated MIMO Coupling (H)")
plt.show()
2. Project and subtract the coherent noise#
Only the coherent, witness-predictable part should be removed. Irreducible target noise should remain.
tsd_fft = tsd.fft()
H = H_lowres.interpolate(tsd_fft["MAIN"].frequencies)
X_mat = FrequencySeriesDict({k: tsd_fft[k] for k in aux_names}).to_matrix()
# Projection estimates the part of MAIN that can be reconstructed from the witnesses.
Y_proj = (H @ X_mat)[0, 0]
projected_ts = Y_proj.ifft()
assert projected_ts.unit == tsd["MAIN"].unit
assert projected_ts.dt == tsd["MAIN"].dt
assert projected_ts.t0 == tsd["MAIN"].t0
assert projected_ts.size == tsd["MAIN"].size
cleaned_ts = tsd["MAIN"] - projected_ts.bandpass(100, 130).real
asd_raw = tsd["MAIN"].asd(fftlength=4.0)
asd_cleaned = cleaned_ts.asd(fftlength=4.0)
plt.figure(figsize=(10, 6))
plt.loglog(asd_raw, label="Original Main")
plt.loglog(asd_cleaned, label="Cleaned")
plt.xlim(10, 1000)
plt.legend()
plt.grid(True, which="both")
plt.show()
