Analysis

Stability: Stable

Coupling Function Analysis Module for gwexpy.

Estimates the coupling function (CF) with flexible threshold strategies. Threshold strategies are defined in gwexpy.analysis.threshold. The result container is defined in gwexpy.analysis.coupling_result.

class gwexpy.analysis.coupling.CouplingFunctionAnalysis[source]

Bases: object

Estimate Coupling Functions (CFs).

classmethod from_time_windows(data: TimeSeriesDict, bkg_window: tuple[float, float], inj_window: tuple[float, float], witness: str | None = None, fftlength: float = 2.0, overlap: float = 0, threshold_strategy: ThresholdStrategy | float = 3.0, frange: tuple[float, float] | None = None, percentile_factor: float = 2.6, bkg_stride: float | None = None, memory_limit: int = 2147483648, n_jobs: int | None = None, **kwargs: Any) → CouplingResult | dict[str, CouplingResult][source]

Compute Coupling Function by specifying time windows explicitly.

Parameters:

  • data (TimeSeriesDict) – Input data containing the full time range (both background and injection).

  • bkg_window (tuple of float) – GPS time window (t_start, t_end) for the background period.

  • inj_window (tuple of float) – GPS time window (t_start, t_end) for the injection period.

  • witness (str, optional) – The name of the witness channel. If None, the first channel is used.

  • fftlength (float) – FFT length in seconds.

  • overlap (float) – Overlap in seconds (default 0).

  • threshold_strategy (ThresholdStrategy or float) – Threshold strategy used for both Witness and Target. If a float is provided, it is interpreted as a RatioThreshold.

  • frange (tuple of float, optional) – Frequency range (fmin, fmax) to evaluate.

  • percentile_factor (float) – Correction factor for PercentileThreshold (Appendix B.1).

  • bkg_stride (float, optional) – Stride in seconds for building the background SegmentTable.

  • memory_limit (int) – Maximum memory in bytes for the background SegmentTable (default 2 GB).

  • n_jobs (int, optional) – Number of parallel jobs. None means 1; -1 uses all cores.

  • **kwargs – Additional keyword arguments forwarded to PSD computation.

Returns:

A CouplingResult for a single target, or a dict for multiple targets.

Return type:

CouplingResult or dict of CouplingResult

Examples

>>> from gwexpy.analysis.coupling import CouplingFunctionAnalysis
>>> from gwexpy import TimeSeriesDict, TimeSeries
>>> import numpy as np
>>> t = np.linspace(1000000000, 1000000400, 40000)
>>> data = TimeSeriesDict()
>>> data['V1:ENV_WIT'] = TimeSeries(np.random.randn(40000), t0=1000000000, sample_rate=100)
>>> data['V1:ITM_SIGNAL'] = TimeSeries(np.random.randn(40000), t0=1000000000, sample_rate=100)
>>> result = CouplingFunctionAnalysis.from_time_windows(
...     data,
...     bkg_window=(1000000000, 1000000100),
...     inj_window=(1000000200, 1000000300),
...     witness="V1:ENV_WIT",
...     fftlength=4.0,
... )
classmethod from_time_windows_batch(data: TimeSeriesDict, bkg_window: tuple[float, float], inj_windows: list[tuple[float, float]], witness: str | None = None, fftlength: float = 2.0, overlap: float = 0, threshold_strategy: ThresholdStrategy | float = 3.0, frange: tuple[float, float] | None = None, percentile_factor: float = 2.6, bkg_stride: float | None = None, memory_limit: int = 2147483648, n_jobs: int | None = None, **kwargs: Any) → list[CouplingResult | dict[str, CouplingResult]][source]

Compute Coupling Functions for multiple injection windows in a batch.

Parameters:

  • data (TimeSeriesDict) – Input data containing the full time range.

  • bkg_window (tuple of float) – GPS time window (t_start, t_end) for the background period. Commonly used across all injection windows.

  • inj_windows (list of tuple of float) – List of injection windows, each as (t_start, t_end).

  • witness (str, optional) – The name of the witness channel. If None, the first channel is used.

  • fftlength (float) – FFT length in seconds.

  • overlap (float) – Overlap in seconds (default 0).

  • threshold_strategy (ThresholdStrategy or float) – Threshold strategy used for both Witness and Target. If a float is provided, it is interpreted as a RatioThreshold.

  • frange (tuple of float, optional) – Frequency range (fmin, fmax) to evaluate.

  • percentile_factor (float) – Correction factor for PercentileThreshold (Appendix B.1).

  • bkg_stride (float, optional) – Stride in seconds for building the background SegmentTable.

  • memory_limit (int) – Maximum memory in bytes for the background SegmentTable (default 2 GB).

  • n_jobs (int, optional) – Number of parallel jobs. None means 1; -1 uses all cores.

  • **kwargs – Additional keyword arguments forwarded to PSD computation.

Returns:

List of results corresponding to each injection window.

Return type:

list of CouplingResult or dict of CouplingResult

auto_calibrate_percentile_factor(raw_bkg: TimeSeries, target: TimeSeries, fftlength: float, overlap: float = 0, percentile: float = 99.7, factor_range: tuple[float, float] = (0.1, 10.0), **kwargs: Any) → dict[str, Any][source]

Automatically calibrate the percentile correction factor (Appendix B).

Finds the factor c such that the background noise distribution best matches the expected statistical floor, aiming for reduced chi2 β‰ˆ 1.

compute(data_inj: TimeSeriesDict, data_bkg: TimeSeriesDict, fftlength: float, witness: str | None = None, frange: tuple[float, float] | None=None, overlap: float = 0, threshold_witness: ThresholdStrategy = <gwexpy.analysis.threshold.RatioThreshold object>, threshold_target: ThresholdStrategy = <gwexpy.analysis.threshold.RatioThreshold object>, percentile_factor: float = 2.6, bkg_stride: float | None = None, memory_limit: int = 2147483648, n_jobs: int | None = None, **kwargs: object) → CouplingResult | dict[str, CouplingResult][source]

Compute Coupling Function(s) from TimeSeriesDicts.

Parameters:

  • data_inj (TimeSeriesDict) – Injection data (Witness + Targets).

  • data_bkg (TimeSeriesDict) – Background data (Witness + Targets).

  • fftlength (float) – FFT length in seconds.

  • witness (str, optional) – The name (key) of the witness channel. If None, the FIRST channel in data_inj is used.

  • frange (tuple of float, optional) – Frequency range (fmin, fmax) to evaluate CF and CF upper limit. Values outside the range are set to NaN.

  • overlap (float, optional) – Overlap in seconds (default 0).

  • threshold_witness (ThresholdStrategy) – Strategy to determine if Witness is excited.

  • threshold_target (ThresholdStrategy) – Strategy to determine if Target is excited.

  • percentile_factor (float, optional) – Correction factor used by PercentileThreshold.

  • bkg_stride (float, optional) – Stride in seconds for background SegmentTable construction.

  • memory_limit (int, optional) – Maximum background SegmentTable memory budget in bytes.

  • n_jobs (int, optional) – Number of jobs for parallel processing. None means 1 unless in a joblib.parallel_config context. -1 means using all processors.

  • **kwargs – Additional keyword arguments forwarded to PSD computation.

class gwexpy.analysis.coupling.CouplingResult(cf: FrequencySeries, psd_witness_inj: FrequencySeries, psd_witness_bkg: FrequencySeries, psd_target_inj: FrequencySeries, psd_target_bkg: FrequencySeries, valid_mask: np.ndarray, witness_name: str, target_name: str, cf_ul: FrequencySeries | None = None, ts_witness_bkg: TimeSeries | None = None, ts_target_bkg: TimeSeries | None = None, ts_witness_inj: TimeSeries | None = None, ts_target_inj: TimeSeries | None = None, fftlength: float | None = None, overlap: float | None = None)[source]

Bases: object

Result object for a SINGLE Witness -> Target pair.

This object stores all Amplitude Spectral Densities (ASDs), the resulting Coupling Function (CF), and metadata needed to reproduce and visualize the injection analysis for a given witness-target pair.

cf

The measured Coupling Function (CF) magnitude spectrum.

Type:

FrequencySeries

cf_ul

The upper limit of the Coupling Function magnitude, if applicable.

Type:

FrequencySeries or None

psd_witness_inj

The Power Spectral Density (PSD) of the witness channel during injection.

Type:

FrequencySeries

psd_witness_bkg

The background PSD of the witness channel.

Type:

FrequencySeries

psd_target_inj

The PSD of the target channel during injection.

Type:

FrequencySeries

psd_target_bkg

The background PSD of the target channel.

Type:

FrequencySeries

valid_mask

Boolean array indicating where the CF measurement is considered valid.

Type:

numpy.ndarray

witness_name

Name of the witness channel.

Type:

str

target_name

Name of the target channel.

Type:

str

ts_witness_bkg

The original background TimeSeries for the witness channel.

Type:

TimeSeries or None

ts_target_bkg

The original background TimeSeries for the target channel.

Type:

TimeSeries or None

ts_witness_inj

The original injection TimeSeries for the witness channel.

Type:

TimeSeries or None

ts_target_inj

The original injection TimeSeries for the target channel.

Type:

TimeSeries or None

fftlength

The FFT length (in seconds) used to compute the spectra.

Type:

float or None

overlap

The overlap duration (in seconds) used to compute the spectra.

Type:

float or None

property frequencies: IndexLike

Return the frequency axis of the coupling function result.

plot_cf(figsize: tuple[float, float] | None = None, xlim: tuple[float, float] | None = None, **kwargs: object) → Plot[source]

Plot the Coupling Function and its Upper Limit.

plot(figsize: tuple[float, float] = (10, 12), xlim: tuple[float, float] | None = None) → Plot[source]

Create a diagnostic plot showing ASDs and the resulting CF.

spectral_stats() → SpectralStats[source]

Return summary statistical information for background ASD.

Currently, it constructs the statistics from the background/injection PSDs stored in the CouplingResult, using background ASD as mean and the absolute difference from injection as sigma. Sigma has an epsilon lower bound to avoid division by zero.

to_csv(filepath: str | PathLike) → None[source]

Save coupling factors to a CSV file.

Columns: frequency, cf, cf_ul, significance, inj_asd, bkg_asd

classmethod from_csv(filepath: str | PathLike) → CouplingResult[source]

Restore a CouplingResult from a CSV file (for round-tripping).

Only supports files exported by to_csv(). Restores cf, cf_ul, and ASDs; other fields are filled with minimal dummy values.

to_txt(filepath: str | PathLike) → None[source]

Save coupling factors in NInjA.py-compatible text format.

Format:

# frequency(Hz) coupling_factor uncertainty significance
1.0 0.123 0.045 2.73
...

The uncertainty is calculated as cf_ul - cf (NaN if cf_ul is missing).

plot_significance(threshold: float = 3.0, freq_min: float | None = None, freq_max: float | None = None, figsize: tuple[float, float] = (12, 6)) → Any[source]

Plot the significance spectrum ((ASD_inj - ASD_bkg) / ASD_bkg vs. frequency).

Parameters:

  • threshold (float) – Value for the horizontal threshold line (default: 3.0). If <= 0, no line is drawn.

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • figsize (tuple of float, optional) – Figure size.

Returns:

The resulting figure instance.

Return type:

matplotlib.figure.Figure

plot_asdgram(freq_min: float | None = None, freq_max: float | None = None, vmin: float | None = None, vmax: float | None = None, figsize: tuple[float, float] = (14, 6)) → Any[source]

Plot ASD spectrogram with percentile overlays (2-column layout).

The layout shows the injection ASD spectrogram in the left panel and the background ASD spectrogram in the right panel. Both are overlaid with 50th, 90th, and 99th percentile lines.

Parameters:

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • vmin (float, optional) – Colorbar scale limits. Automatic if None.

  • vmax (float, optional) – Colorbar scale limits. Automatic if None.

  • figsize (tuple of float, optional) – Figure size.

Raises:

ValueError – If ts_witness_inj or ts_witness_bkg is None.

plot_rms(fmin: float | None = None, fmax: float | None = None, fftlength: float | None = None, overlap: float | None = None, channels: str = 'both', show_windows: bool = True, figsize: tuple[float, float] = (12, 4), **kwargs: Any) → Any[source]

Plot rolling RMS time-series (band-limited).

Displays the rolling RMS for witness and target channels over time. If show_windows=True, background (grey) and injection (red) intervals are shaded using axvspan.

RMS definition (root frequency-integral of PSD):

RMS(t) = sqrt( trapz(PSD(t, f), f) )  for f in [fmin, fmax]
Parameters:

  • fmin (float, optional) – Lower frequency limit for RMS calculation [Hz]. Defaults to 0 Hz.

  • fmax (float, optional) – Upper frequency limit for RMS calculation [Hz]. Defaults to Nyquist.

  • fftlength (float, optional) – FFT length for the spectrogram [s]. Uses self.fftlength if None.

  • overlap (float, optional) – Overlap for the spectrogram [s]. Uses self.overlap if None.

  • channels ({"witness", "target", "both"}) – Which channels to plot. β€œboth” creates a 2-panel layout.

  • show_windows (bool) – If True, shades background and injection intervals using axvspan (default: True).

  • figsize (tuple of float, optional) – Figure size.

  • **kwargs – Additional keyword arguments for matplotlib’s plot().

Returns:

The created figure.

Return type:

matplotlib.figure.Figure

Raises:

ValueError – If the required TimeSeries for the selected channels is None, or if the β€˜channels’ argument is invalid.

plot_snrgram(freq_min: float | None = None, freq_max: float | None = None, snrmax: float = 100.0, figsize: tuple[float, float] = (12, 6)) → Any[source]

Plot SNR spectrogram ((ASD_inj - median_bkg) / median_bkg).

Parameters:

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • snrmax (float) – Colorbar maximum limit (clamp value). Defaults to 100.0.

  • figsize (tuple of float, optional) – Figure size.

Raises:

ValueError – If ts_witness_inj or ts_witness_bkg is None.

classmethod from_txt(filepath: str | PathLike) → CouplingResult[source]

Restore a CouplingResult from a plain-text NInjA compatible file.

Only supports files exported by to_txt(). Restores cf and cf_ul; other fields are filled with minimal dummy values.

class gwexpy.analysis.coupling.PercentileThreshold(percentile: float = 99.7, factor: float = 2.6, freq_align: str = 'clip')[source]

Bases: ThresholdStrategy

Threshold strategy based on empirical percentile of background distribution.

This strategy follows Appendix B of the PEM injection paper, using the 99.7th percentile of background segments and a correction factor to account for finite-averaging and χ² distribution scaling.

Parameters:

  • percentile (float, default=99.7) – The percentile of the background distribution (0-100). 99.7% equivalent to 3-sigma for Gaussian noise.

  • factor (float, default=2.6) – Correction factor (multiplier) for the percentile value. The value 2.6 is recommended in Appendix B.1 to set reduced χ² β‰ˆ 1.

  • freq_align ({'clip', 'interpolate'}, default='clip') – Frequency alignment strategy for background segments.

check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the mask where injection PSD exceeds the percentile threshold.

threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: Any) → np.ndarray[source]

Return the PSD threshold from the empirical background distribution.

class gwexpy.analysis.coupling.RatioThreshold(ratio: float = 2.0)[source]

Bases: ThresholdStrategy

Checks if P_inj > ratio * P_bkg_mean.

Statistical Assumptions:

  • No specific statistical distribution is assumed.

  • Tests if injection power exceeds the background level by a fixed factor.

Usage:

  • Best for simple, physical excess screening where precise statistical significance is less critical.

  • Extremely fast as it requires no variance estimation.

check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the mask where injection PSD exceeds the ratio threshold.

threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the PSD threshold implied by the fixed ratio.

class gwexpy.analysis.coupling.SigmaThreshold(sigma: float = 3.0)[source]

Bases: ThresholdStrategy

Checks if P_inj > P_bkg + sigma * std_error.

Statistical Assumptions

  • Background Power Spectral Density (PSD) at each bin approximately follows a Gaussian distribution (valid when n_avg is sufficiently large).

  • The parameter n_avg represents the number of independent averages (e.g., in Welch’s method).

  • Assumes standard deviation of the noise reduces as 1 / sqrt(n_avg).

Meaning of Threshold

  • threshold = mean + sigma * (mean / sqrt(n_avg))

  • This is a statistical significance test, NOT a physical upper limit.

  • It identifies frequencies where the injection is statistically distinguishable from background variance.

Gaussian Approximation Validity

Welch PSD estimates follow a χ² distribution with 2K degrees of freedom (K = n_avg). The Gaussian approximation is valid when K β‰₯ 10 (approximately).

For K < 10, consider: - Using PercentileThreshold (empirical distribution, no Gaussian assumption) - Increasing FFT averaging by using longer data or shorter fftlength

References

  • Welch, P.D. (1967): PSD estimation via overlapped segment averaging

  • Bendat & Piersol, Random Data (4th ed., 2010), Ch. 11

Warning

This method relies heavily on the Gaussian and stationary assumptions. It may be unreliable if: - The background contains significant non-Gaussian features (glitches) - n_avg is small (< ~10), where the central limit theorem has not converged - There are strong spectral lines (non-stationary or deterministic signals)

In such cases, PercentileThreshold is recommended as it uses the empirical distribution.

check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the mask where injection PSD exceeds the sigma threshold.

threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the PSD threshold implied by the Gaussian approximation.

class gwexpy.analysis.coupling.ThresholdStrategy[source]

Bases: ABC

Abstract base class for excess detection strategies.

abstractmethod check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return a boolean mask where P_inj is considered β€˜excess’.

Parameters:

  • psd_inj (FrequencySeries) – PSD of the injection data.

  • psd_bkg (FrequencySeries) – PSD of the background data (usually mean or median).

  • raw_bkg (TimeSeries, optional) – Raw background time series data. Required for PercentileThreshold to calculate the distribution of PSDs across segments.

  • **kwargs – Strategy-specific options forwarded by the caller.

abstractmethod threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the PSD threshold values used by this strategy.

gwexpy.analysis.coupling.estimate_bkg_mem_bytes(duration: float, fftlength: float, stride: float, fs: float) → tuple[int, int][source]

Estimate memory required for a background SegmentTable.

Returns:

  • bytes_est (int) – Estimated memory in bytes (including 20% overhead).

  • n_rows (int) – Number of rows.

gwexpy.analysis.coupling.estimate_coupling(data_inj: TimeSeriesDict, data_bkg: TimeSeriesDict | None = None, fftlength: float = 2.0, witness: str | None = None, frange: tuple[float, float] | None = None, overlap: float = 0, threshold_witness: ThresholdStrategy | float = 25.0, threshold_target: ThresholdStrategy | float = 4.0, percentile_factor: float = 2.6, bkg_stride: float | None = None, memory_limit: int = 2147483648, n_jobs: int | None = None, bkg_window: tuple[float, float] | None = None, inj_window: tuple[float, float] | None = None, **kwargs: Any) → CouplingResult | dict[str, CouplingResult][source]

Estimate CF from injection and background data.

Parameters:

  • data_inj (TimeSeriesDict) – Injection data (Witness + Targets). If bkg_window / inj_window are specified, please pass data covering the full time range. In this case, data_bkg is not required.

  • data_bkg (TimeSeriesDict, optional) – Background data (Witness + Targets). Can be omitted if bkg_window is specified.

  • fftlength (float) – FFT length in seconds.

  • witness (str, optional) – The name (key) of the witness channel. If None, the FIRST channel in data_inj is used.

  • frange (tuple of float, optional) – Frequency range (fmin, fmax) to evaluate CF and CF upper limit. Values outside the range are set to NaN.

  • overlap (float, optional) – Overlap in seconds (default 0).

  • threshold_witness (ThresholdStrategy or float) – Strategy to determine if Witness is excited. If a float is given, it is interpreted as a ratio for RatioThreshold.

  • threshold_target (ThresholdStrategy or float) – Strategy to determine if Target is excited. If a float is given, it is interpreted as a ratio for RatioThreshold.

  • percentile_factor (float, optional) – Correction factor for PercentileThreshold (Appendix B.1). Default is 2.6.

  • bkg_stride (float, optional) – Stride in seconds for background SegmentTable construction. Defaults to fftlength.

  • memory_limit (int, optional) – Maximum memory in bytes for the background SegmentTable. Default is 2 GB.

  • n_jobs (int, optional) – Number of jobs for parallel processing. None means 1; -1 means all processors.

  • bkg_window (tuple of float, optional) – GPS time window (t_start, t_end) for the background period. If specified, background data is cropped from data_inj, so data_bkg is not required. Recommended to be used with inj_window.

  • inj_window (tuple of float, optional) – GPS time window (t_start, t_end) for the injection period. Works as part of the time-window API when used with bkg_window.

  • **kwargs – Additional keyword arguments forwarded to PSD computation.

CouplingResult: result container for a single Witness -> Target pair.

class gwexpy.analysis.coupling_result.CouplingResult(cf: FrequencySeries, psd_witness_inj: FrequencySeries, psd_witness_bkg: FrequencySeries, psd_target_inj: FrequencySeries, psd_target_bkg: FrequencySeries, valid_mask: np.ndarray, witness_name: str, target_name: str, cf_ul: FrequencySeries | None = None, ts_witness_bkg: TimeSeries | None = None, ts_target_bkg: TimeSeries | None = None, ts_witness_inj: TimeSeries | None = None, ts_target_inj: TimeSeries | None = None, fftlength: float | None = None, overlap: float | None = None)[source]

Bases: object

Result object for a SINGLE Witness -> Target pair.

This object stores all Amplitude Spectral Densities (ASDs), the resulting Coupling Function (CF), and metadata needed to reproduce and visualize the injection analysis for a given witness-target pair.

cf

The measured Coupling Function (CF) magnitude spectrum.

Type:

FrequencySeries

cf_ul

The upper limit of the Coupling Function magnitude, if applicable.

Type:

FrequencySeries or None

psd_witness_inj

The Power Spectral Density (PSD) of the witness channel during injection.

Type:

FrequencySeries

psd_witness_bkg

The background PSD of the witness channel.

Type:

FrequencySeries

psd_target_inj

The PSD of the target channel during injection.

Type:

FrequencySeries

psd_target_bkg

The background PSD of the target channel.

Type:

FrequencySeries

valid_mask

Boolean array indicating where the CF measurement is considered valid.

Type:

numpy.ndarray

witness_name

Name of the witness channel.

Type:

str

target_name

Name of the target channel.

Type:

str

ts_witness_bkg

The original background TimeSeries for the witness channel.

Type:

TimeSeries or None

ts_target_bkg

The original background TimeSeries for the target channel.

Type:

TimeSeries or None

ts_witness_inj

The original injection TimeSeries for the witness channel.

Type:

TimeSeries or None

ts_target_inj

The original injection TimeSeries for the target channel.

Type:

TimeSeries or None

fftlength

The FFT length (in seconds) used to compute the spectra.

Type:

float or None

overlap

The overlap duration (in seconds) used to compute the spectra.

Type:

float or None

property frequencies: IndexLike

Return the frequency axis of the coupling function result.

plot_cf(figsize: tuple[float, float] | None = None, xlim: tuple[float, float] | None = None, **kwargs: object) → Plot[source]

Plot the Coupling Function and its Upper Limit.

plot(figsize: tuple[float, float] = (10, 12), xlim: tuple[float, float] | None = None) → Plot[source]

Create a diagnostic plot showing ASDs and the resulting CF.

spectral_stats() → SpectralStats[source]

Return summary statistical information for background ASD.

Currently, it constructs the statistics from the background/injection PSDs stored in the CouplingResult, using background ASD as mean and the absolute difference from injection as sigma. Sigma has an epsilon lower bound to avoid division by zero.

to_csv(filepath: str | PathLike) → None[source]

Save coupling factors to a CSV file.

Columns: frequency, cf, cf_ul, significance, inj_asd, bkg_asd

classmethod from_csv(filepath: str | PathLike) → CouplingResult[source]

Restore a CouplingResult from a CSV file (for round-tripping).

Only supports files exported by to_csv(). Restores cf, cf_ul, and ASDs; other fields are filled with minimal dummy values.

to_txt(filepath: str | PathLike) → None[source]

Save coupling factors in NInjA.py-compatible text format.

Format:

# frequency(Hz) coupling_factor uncertainty significance
1.0 0.123 0.045 2.73
...

The uncertainty is calculated as cf_ul - cf (NaN if cf_ul is missing).

plot_significance(threshold: float = 3.0, freq_min: float | None = None, freq_max: float | None = None, figsize: tuple[float, float] = (12, 6)) → Any[source]

Plot the significance spectrum ((ASD_inj - ASD_bkg) / ASD_bkg vs. frequency).

Parameters:

  • threshold (float) – Value for the horizontal threshold line (default: 3.0). If <= 0, no line is drawn.

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • figsize (tuple of float, optional) – Figure size.

Returns:

The resulting figure instance.

Return type:

matplotlib.figure.Figure

plot_asdgram(freq_min: float | None = None, freq_max: float | None = None, vmin: float | None = None, vmax: float | None = None, figsize: tuple[float, float] = (14, 6)) → Any[source]

Plot ASD spectrogram with percentile overlays (2-column layout).

The layout shows the injection ASD spectrogram in the left panel and the background ASD spectrogram in the right panel. Both are overlaid with 50th, 90th, and 99th percentile lines.

Parameters:

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • vmin (float, optional) – Colorbar scale limits. Automatic if None.

  • vmax (float, optional) – Colorbar scale limits. Automatic if None.

  • figsize (tuple of float, optional) – Figure size.

Raises:

ValueError – If ts_witness_inj or ts_witness_bkg is None.

plot_rms(fmin: float | None = None, fmax: float | None = None, fftlength: float | None = None, overlap: float | None = None, channels: str = 'both', show_windows: bool = True, figsize: tuple[float, float] = (12, 4), **kwargs: Any) → Any[source]

Plot rolling RMS time-series (band-limited).

Displays the rolling RMS for witness and target channels over time. If show_windows=True, background (grey) and injection (red) intervals are shaded using axvspan.

RMS definition (root frequency-integral of PSD):

RMS(t) = sqrt( trapz(PSD(t, f), f) )  for f in [fmin, fmax]
Parameters:

  • fmin (float, optional) – Lower frequency limit for RMS calculation [Hz]. Defaults to 0 Hz.

  • fmax (float, optional) – Upper frequency limit for RMS calculation [Hz]. Defaults to Nyquist.

  • fftlength (float, optional) – FFT length for the spectrogram [s]. Uses self.fftlength if None.

  • overlap (float, optional) – Overlap for the spectrogram [s]. Uses self.overlap if None.

  • channels ({"witness", "target", "both"}) – Which channels to plot. β€œboth” creates a 2-panel layout.

  • show_windows (bool) – If True, shades background and injection intervals using axvspan (default: True).

  • figsize (tuple of float, optional) – Figure size.

  • **kwargs – Additional keyword arguments for matplotlib’s plot().

Returns:

The created figure.

Return type:

matplotlib.figure.Figure

Raises:

ValueError – If the required TimeSeries for the selected channels is None, or if the β€˜channels’ argument is invalid.

plot_snrgram(freq_min: float | None = None, freq_max: float | None = None, snrmax: float = 100.0, figsize: tuple[float, float] = (12, 6)) → Any[source]

Plot SNR spectrogram ((ASD_inj - median_bkg) / median_bkg).

Parameters:

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • snrmax (float) – Colorbar maximum limit (clamp value). Defaults to 100.0.

  • figsize (tuple of float, optional) – Figure size.

Raises:

ValueError – If ts_witness_inj or ts_witness_bkg is None.

classmethod from_txt(filepath: str | PathLike) → CouplingResult[source]

Restore a CouplingResult from a plain-text NInjA compatible file.

Only supports files exported by to_txt(). Restores cf and cf_ul; other fields are filled with minimal dummy values.

class gwexpy.analysis.coupling_result.CouplingResultCollection(mapping: dict[str, Any] | None = None)[source]

Bases: dict

Container for managing multiple CouplingResult objects.

Examples

>>> from gwexpy.analysis.coupling_result import CouplingResultCollection
>>> results = CouplingResultCollection()
>>> isinstance(results, dict)
True
>>> # results['WIT-TGT'] = coupling_result_1
>>> # results.to_summary_csv("summary.csv")
to_summary_csv(filepath: str | PathLike) → None[source]

Export all results to a single summary CSV file.

Columns: channel_pair, frequency, cf, cf_ul, significance, inj_asd, bkg_asd

plot_comparison(freq_min: float | None = None, freq_max: float | None = None, threshold: float = 3.0, figsize: tuple[float, float] = (12, 8)) → Any[source]

Plot a comparison of multiple coupling factors.

Parameters:

  • freq_min (float, optional) – Frequency range to display [Hz].

  • freq_max (float, optional) – Frequency range to display [Hz].

  • threshold (float) – Value for the horizontal threshold line (default: 3.0).

  • figsize (tuple of float, optional) – Figure size.

Returns:

The created figure instance.

Return type:

matplotlib.figure.Figure

Response Function Analysis Module for gwexpy.

This module implements the Response Function Model (RFM) based on Stepped Sine (Discrete) Injections. It prioritizes statistical significance by calculating averaged ASDs for each stable frequency step.

class gwexpy.analysis.response.ResponseFunctionResult(spectrogram_inj: Spectrogram, spectrogram_bkg: Spectrogram, injected_freqs: ndarray, step_times: ndarray, coupling_factors: ndarray, witness_name: str, target_name: str, table: SegmentTable | None = None)[source]

Bases: object

Result object for Stepped Sine Response Function Analysis.

Stores the full spectral data for all injection steps efficiently.

spectrogram_inj

The measured target-channel ASD spectrogram during the injection steps. Its frequency axis is the target ASD axis and can include bins above the witness Nyquist frequency; those high-frequency bins are raw target ASD values, not coupling-factor estimates.

Type:

Spectrogram

spectrogram_bkg

The target-channel background ASD spectrogram for the injection steps, on the same target frequency axis as spectrogram_inj.

Type:

Spectrogram

injected_freqs

Array of injected frequencies [Hz] corresponding to each step.

Type:

numpy.ndarray

step_times

Array of GPS start times corresponding to each step.

Type:

numpy.ndarray

coupling_factors

Representative coupling factors evaluated at each injection frequency.

Type:

numpy.ndarray

witness_name

Name of the witness channel.

Type:

str

target_name

Name of the target channel.

Type:

str

table

Underlying segment table containing detailed analysis information.

Type:

SegmentTable, optional

spectrogram_inj: Spectrogram
spectrogram_bkg: Spectrogram
injected_freqs: ndarray
step_times: ndarray
coupling_factors: ndarray
witness_name: str
target_name: str
table: SegmentTable | None = None
plot(ax: Axes | None = None, **kwargs: Any) → Axes[source]

Plot the Coupling Factor vs Injected Frequency (The Transfer Function).

plot_map(ax: Axes | None = None, **kwargs: Any) → Axes[source]

Plot the 2D Response Map (Injected Freq vs Target Spectrum).

Useful to check for non-linear couplings.

plot_projection_summary(freq_min: float | None = None, freq_max: float | None = None, figsize: tuple[float, float] = (14, 6)) → Any[source]

Plot the ASD spectra across injection steps.

Create a single-panel figure that overlays the ASD (or PSD) spectra for all injection steps. This visualization helps inspect step-to-step consistency and detect non-linear or time-dependent response.

Parameters:

  • freq_min (float, optional) – Frequency range in Hz to include in the plot. If None, the full spectrogram frequency range is used.

  • freq_max (float, optional) – Frequency range in Hz to include in the plot. If None, the full spectrogram frequency range is used.

  • figsize (tuple, optional) – Matplotlib figure size as (width, height). Defaults to (14, 6).

Returns:

The created figure instance.

Return type:

matplotlib.figure.Figure

plot_response_matrix(freq_min: float | None = None, freq_max: float | None = None, figsize: tuple[float, float] = (14, 10)) → Any[source]

Plot the 2D response matrix and diagnostic cross-sections (3-panel layout).

The figure consists of:

  • Main panel: time (injection step) vs response frequency, color = ASD amplitude.

  • Right-side panel: temporal evolution (time series) at the central frequency bin.

  • Top panel: frequency profile at the central time step.

Parameters:

  • freq_min (float, optional) – Frequency range in Hz to visualize.

  • freq_max (float, optional) – Frequency range in Hz to visualize.

  • figsize (tuple, optional) – Figure size.

Returns:

The created figure.

Return type:

matplotlib.figure.Figure

plot_snapshot(freq: float | None = None, step_index: int | None = None, ax: Axes | None = None) → Axes[source]

Plot ASDs and Upper Limits for a SPECIFIC injection step.

class gwexpy.analysis.response.ResponseFunctionAnalysis[source]

Bases: object

Analysis engine for Stepped Sine Injections.

compute(witness: TimeSeries, target: TimeSeries, segments: list[tuple[float, float, float]] | None = None, fftlength: float = 4.0, overlap: float = 0.0, auto_detect: bool = True, snr_threshold: float = 10.0, min_duration: float = 5.0, trim_edge: float = 1.0, witness_bkg: TimeSeries | None = None, target_bkg: TimeSeries | None = None, bkg_window: tuple[float, float] | None = None, n_jobs: int | None = None, memory_limit: int = 2147483648, **kwargs: object) → ResponseFunctionResult[source]

Perform the analysis.

Parameters:

  • witness (TimeSeries) – Witness channel time series.

  • target (TimeSeries) – Target channel time series.

  • segments (list of tuple, optional) – Manually specified injection steps as (t_start, t_end, freq) tuples. If None and auto_detect=True, steps are detected automatically.

  • fftlength (float) – FFT length [seconds].

  • overlap (float) – Overlap [seconds] (default 0).

  • auto_detect (bool) – If True and segments is None, auto-detect injection steps.

  • snr_threshold (float) – SNR threshold for step detection.

  • min_duration (float) – Minimum step duration [seconds] for detection.

  • trim_edge (float) – Edge trim [seconds] for each detected step.

  • witness_bkg (TimeSeries, optional) – Explicit background for the witness channel. Takes precedence over bkg_window.

  • target_bkg (TimeSeries, optional) – Explicit background for the target channel. Takes precedence over bkg_window.

  • bkg_window (tuple of float, optional) – A tuple of (t_start, t_end) GPS times for the background interval. Used if witness_bkg / target_bkg are not specified. The specified time range will be cropped from witness / target and used as background data. Has precedence over auto-detection.

  • n_jobs (int, optional) – Number of parallel jobs.

  • memory_limit (int) – Memory limit for batch processing [bytes].

  • **kwargs (dict) – Additional keyword arguments passed to the spectrum calculation.

gwexpy.analysis.response.detect_step_segments(witness: TimeSeries, fftlength: float = 1.0, snr_threshold: float = 10.0, min_duration: float = 5.0, trim_edge: float = 1.0, freq_tolerance: float = 2.0) → list[tuple[float, float, float]][source]

Automatically detect time segments where the injection frequency is constant.

Parameters:

  • witness (TimeSeries) – Time series data containing the injection signal.

  • fftlength (float, optional) – FFT length for tracking spectrogram (default: 1.0s, giving Ξ”f=1 Hz).

  • snr_threshold (float, optional) – Minimum SNR for injection detection (default: 10.0).

  • min_duration (float, optional) – Minimum duration for a valid step (default: 5.0s).

  • trim_edge (float, optional) – Time to trim from step edges (default: 1.0s).

  • freq_tolerance (float, optional) –

    Maximum frequency change within a step [Hz] (default: 2.0).

    Note: For reliable detection, freq_tolerance should be at least 2 * Ξ”f where Ξ”f = 1/fftlength. With the default fftlength=1.0s (Ξ”f=1 Hz), a strict tolerance of Ξ”f risks false step boundaries due to spectral bin fluctuations. The default 2.0 Hz provides adequate margin.

    Reference: Harris, F.J. β€œOn the use of windows for harmonic analysis with the discrete Fourier transform”, Proc. IEEE 66(1), 1978.

Returns:

List of (start_time, end_time, frequency) tuples for each detected step.

Return type:

list of tuple

gwexpy.analysis.response.estimate_response_function(witness: TimeSeries, target: TimeSeries, segments: list[tuple[float, float, float]] | None = None, fftlength: float = 4.0, overlap: float = 0.0, auto_detect: bool = True, snr_threshold: float = 10.0, min_duration: float = 5.0, trim_edge: float = 1.0, witness_bkg: TimeSeries | None = None, target_bkg: TimeSeries | None = None, bkg_window: tuple[float, float] | None = None, n_jobs: int | None = None, memory_limit: int = 2147483648, **kwargs: Any) → ResponseFunctionResult[source]

Backward-compatible wrapper for ResponseFunctionAnalysis.

This function preserves the historical module-level API used by external callers and examples. New code should prefer ResponseFunctionAnalysis().compute(...).

Threshold strategies for Coupling Function Analysis.

This module provides excess-detection strategies used by CouplingFunctionAnalysis: - ThresholdStrategy: Abstract base class - RatioThreshold: Simple power-ratio test - SigmaThreshold: Gaussian significance test - PercentileThreshold: Empirical percentile (Appendix B)

class gwexpy.analysis.threshold.ThresholdStrategy[source]

Bases: ABC

Abstract base class for excess detection strategies.

abstractmethod check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return a boolean mask where P_inj is considered β€˜excess’.

Parameters:

  • psd_inj (FrequencySeries) – PSD of the injection data.

  • psd_bkg (FrequencySeries) – PSD of the background data (usually mean or median).

  • raw_bkg (TimeSeries, optional) – Raw background time series data. Required for PercentileThreshold to calculate the distribution of PSDs across segments.

  • **kwargs – Strategy-specific options forwarded by the caller.

abstractmethod threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the PSD threshold values used by this strategy.

class gwexpy.analysis.threshold.RatioThreshold(ratio: float = 2.0)[source]

Bases: ThresholdStrategy

Checks if P_inj > ratio * P_bkg_mean.

Statistical Assumptions:

  • No specific statistical distribution is assumed.

  • Tests if injection power exceeds the background level by a fixed factor.

Usage:

  • Best for simple, physical excess screening where precise statistical significance is less critical.

  • Extremely fast as it requires no variance estimation.

check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the mask where injection PSD exceeds the ratio threshold.

threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the PSD threshold implied by the fixed ratio.

class gwexpy.analysis.threshold.SigmaThreshold(sigma: float = 3.0)[source]

Bases: ThresholdStrategy

Checks if P_inj > P_bkg + sigma * std_error.

Statistical Assumptions

  • Background Power Spectral Density (PSD) at each bin approximately follows a Gaussian distribution (valid when n_avg is sufficiently large).

  • The parameter n_avg represents the number of independent averages (e.g., in Welch’s method).

  • Assumes standard deviation of the noise reduces as 1 / sqrt(n_avg).

Meaning of Threshold

  • threshold = mean + sigma * (mean / sqrt(n_avg))

  • This is a statistical significance test, NOT a physical upper limit.

  • It identifies frequencies where the injection is statistically distinguishable from background variance.

Gaussian Approximation Validity

Welch PSD estimates follow a χ² distribution with 2K degrees of freedom (K = n_avg). The Gaussian approximation is valid when K β‰₯ 10 (approximately).

For K < 10, consider: - Using PercentileThreshold (empirical distribution, no Gaussian assumption) - Increasing FFT averaging by using longer data or shorter fftlength

References

  • Welch, P.D. (1967): PSD estimation via overlapped segment averaging

  • Bendat & Piersol, Random Data (4th ed., 2010), Ch. 11

Warning

This method relies heavily on the Gaussian and stationary assumptions. It may be unreliable if: - The background contains significant non-Gaussian features (glitches) - n_avg is small (< ~10), where the central limit theorem has not converged - There are strong spectral lines (non-stationary or deterministic signals)

In such cases, PercentileThreshold is recommended as it uses the empirical distribution.

check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the mask where injection PSD exceeds the sigma threshold.

threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the PSD threshold implied by the Gaussian approximation.

class gwexpy.analysis.threshold.PercentileThreshold(percentile: float = 99.7, factor: float = 2.6, freq_align: str = 'clip')[source]

Bases: ThresholdStrategy

Threshold strategy based on empirical percentile of background distribution.

This strategy follows Appendix B of the PEM injection paper, using the 99.7th percentile of background segments and a correction factor to account for finite-averaging and χ² distribution scaling.

Parameters:

  • percentile (float, default=99.7) – The percentile of the background distribution (0-100). 99.7% equivalent to 3-sigma for Gaussian noise.

  • factor (float, default=2.6) – Correction factor (multiplier) for the percentile value. The value 2.6 is recommended in Appendix B.1 to set reduced χ² β‰ˆ 1.

  • freq_align ({'clip', 'interpolate'}, default='clip') – Frequency alignment strategy for background segments.

check(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: object) → np.ndarray[source]

Return the mask where injection PSD exceeds the percentile threshold.

threshold(psd_inj: FrequencySeries, psd_bkg: FrequencySeries, raw_bkg: TimeSeries | None = None, **kwargs: Any) → np.ndarray[source]

Return the PSD threshold from the empirical background distribution.

SpectralStats: Spectral statistics container (mean, sigma, n_avg).

class gwexpy.analysis.stats.SpectralStats(mean: FrequencySeries, sigma: FrequencySeries, n_avg: int)[source]

Bases: object

Spectral statistics container (mean, sigma, n_avg).

mean

Mean ASD or PSD.

Type:

FrequencySeries

sigma

Standard deviation of the ASD or PSD.

Type:

FrequencySeries

n_avg

Number of samples used for averaging.

Type:

int

mean: FrequencySeries
sigma: FrequencySeries
n_avg: int
significance(mu_inj: FrequencySeries) → FrequencySeries[source]

Calculate statistical significance: (ΞΌ_inj - ΞΌ_bkg) / Οƒ_bkg.

Parameters:

mu_inj (FrequencySeries) – The ASD or PSD during the injection period.

Returns:

Frequency spectrum of the significance: (ΞΌ_inj - ΞΌ_bkg) / Οƒ_bkg.

Return type:

FrequencySeries

to_dict() → dict[str, object][source]

Return the statistical information as a dictionary.

Stat-info utilities: association edge extraction and lightweight graph builder.

gwexpy.analysis.stat_info.association_edges(target: Any, matrix: Any, *, method: str = 'pearson', parallel: int | None = None, threshold: float | None = None, threshold_mode: str = 'abs', topk: int | None = None, return_dataframe: bool = True) → Any[source]

Compute association edges between a target TimeSeries and a TimeSeriesMatrix.

Returns a DataFrame (default) with columns ["source", "target", "score", "row", "col", "channel"].

gwexpy.analysis.stat_info.build_graph(edges: Any, *, backend: str = 'networkx', directed: bool = False, weight: str = 'score') → Any[source]

Build a graph object from association edges.

If backend="none", return edges unchanged.