スペクトル解析 (Spectral)

gwexpy.spectral.bootstrap_spectrogram(spectrogram, n_boot=1000, method='median', average=None, ci=0.68, window='hann', fftlength=None, overlap=None, nfft=None, noverlap=None, block_size=None, rebin_width=None, return_map=False, ignore_nan=True, **kwargs)[source]

Estimate robust ASD from a spectrogram using bootstrap resampling.

Error bars are corrected for correlation between overlapping segments based on the window function’s autocorrelation.

Performance Note: This function utilizes Numba for JIT compilation and parallelization if available, significantly accelerating the resampling process.

Parameters:

  • spectrogram (gwpy.spectrogram.Spectrogram)

  • n_boot (int)

  • method (str) – ‘median’ (default) or ‘mean’.

  • average (str, optional) – Alias for method (‘median’ or ‘mean’) for compatibility.

  • ci (float)

  • window (str or array, optional)

  • fftlength (float or Quantity, optional) – FFT segment length in seconds used to generate the spectrogram. Used for VIF (overlap correlation) correction. Cannot be used with nfft.

  • overlap (float or Quantity, optional) – Overlap between FFT segments in seconds. Used for VIF correction. If None and fftlength is given, defaults to the recommended overlap for the window function. Cannot be used with noverlap.

  • nfft (int, optional) – FFT segment length in samples. Alternative to fftlength for precise sample-based specification. Cannot be used with fftlength.

  • noverlap (int, optional) – Overlap length in samples. Must be used with nfft. Cannot be used with overlap.

  • block_size (float, Quantity, or 'auto', optional) – Duration of blocks for block bootstrap in seconds. Can be specified as float (seconds), Quantity with time units, or ‘auto’. If ‘auto’, estimates size based on overlap ratio (ceil(duration/stride)). If None, perform standard bootstrap with analytical VIF correction.

  • rebin_width (float, optional) – Width to rebin frequencies (in Hz) before bootstrapping.

  • return_map (bool, optional) – If True, return a BifrequencyMap of the covariance.

  • ignore_nan (bool, optional) – If True, ignore NaNs during bootstrap averaging. Default is True.

Return type:

FrequencySeries or (FrequencySeries, BifrequencyMap)

Notes

VIF (Variance Inflation Factor) Application

The VIF correction for overlapping Welch segments is applied only when block_size is not specified (standard bootstrap). When block_size is specified, the VIF correction is disabled (factor=1.0) because the block bootstrap is designed to capture serial correlation implicitly.

This behavior prevents “double correction” where both block structure and analytical VIF would over-inflate variance estimates.

Stationarity Assumption

Moving Block Bootstrap assumes the input data is a stationary process. For non-stationary data (e.g., glitches, transient events, or drifting noise floors), the bootstrap confidence intervals may be biased.

Consider: - Pre-processing to remove non-stationarities - Using shorter analysis segments where stationarity holds - Validating stationarity with statistical tests

Block Size Selection

If block_size=None, standard (iid) bootstrap is performed, which ignores serial correlation from Welch overlap. Recommended to set block_size to at least the stride length (segment spacing) for time-correlated spectrograms.

Examples

>>> from gwexpy.spectrogram import Spectrogram
>>> from gwexpy.spectral import bootstrap_spectrogram
>>> import numpy as np
>>> from astropy import units as u
>>>
>>> # Create synthetic spectrogram
>>> np.random.seed(42)
>>> spec_data = np.random.random((100, 50))
>>> spec = Spectrogram(spec_data, dt=1.0*u.s, f0=10*u.Hz, df=1*u.Hz)
>>>
>>> # Bootstrap with time-based parameters (high-level API)
>>> result = bootstrap_spectrogram(
...     spec,
...     n_boot=100,
...     fftlength=4.0,    # 4 seconds
...     overlap=2.0,      # 2 seconds
...     block_size=2.0,   # 2 seconds block
...     window='hann',
...     method='median'
... )
>>> print(result.value.shape)
(50,)
>>>
>>> # Alternatively, use sample-based parameters (low-level API)
>>> result2 = bootstrap_spectrogram(
...     spec,
...     n_boot=100,
...     nfft=4,        # 4 samples (dt=1s, so 4 seconds)
...     noverlap=2,    # 2 samples overlap
...     window='hann',
...     method='median'
... )
>>> print(result2.value.shape)
(50,)

References

See also

Tutorial: Bootstrap Spectrogram Analysis

gwexpy.spectral.calculate_correlation_factor(window, nperseg, noverlap, n_blocks)[source]

Calculate the variance inflation factor for Welch’s method with overlap.

This computes the correction factor by numerically calculating the normalized squared autocorrelation of the window function. The correction accounts for reduced effective degrees of freedom due to correlated (overlapping) segments.

Formula:

factor = sqrt(1 + 2 * sum_{k=1}^{M-1} (1 - k/M) * abs(rho_window(k * S))**2)

where rho_window is the normalized autocorrelation of the window.

Notes

This VIF (Variance Inflation Factor) formula follows Percival & Walden (1993), Eq.(56). It corrects for variance inflation due to overlap in Welch’s method.

Important: This is NOT the regression VIF (1/(1-R²)) used for multicollinearity diagnosis. The name collision has caused confusion, but the implementation is correct for spectral analysis.

This formula has been validated through unit tests and independent technical review (2026-02-01). The “statistically invalid” critique from some models was based on confusion with regression VIF.

References

  • Percival, D.B. & Walden, A.T., Spectral Analysis for Physical Applications (1993), Ch. 7.3.2, Eq.(56)

  • Bendat, J.S. & Piersol, A.G., Random Data (4th ed., 2010)

  • Ingram, A. (2019), Error formulae for the energy-dependent cross-spectrum

Parameters:

  • window (str, tuple, or array_like) – Window function used for FFT.

  • nperseg (int) – Length of each segment (N_fft).

  • noverlap (int) – Number of overlapping samples.

  • n_blocks (int) – Number of segments available for averaging.

Returns:

Multiplicative correction factor for the standard error.

Return type:

float

gwexpy.spectral.estimate_psd(timeseries, *, fftlength=None, overlap=None, window='hann', method='median', **kwargs)[source]

Estimate the one-sided PSD for a regular TimeSeries.

This wraps TimeSeries.psd without changing the underlying algorithm. It enforces NaN-free input and checks that fftlength does not exceed the data duration so that the frequency axis (df = 1 / fftlength) is well-defined.

Parameters:

  • timeseries (gwexpy.timeseries.TimeSeries) – Input series for PSD estimation.

  • fftlength (float or Quantity, optional) – FFT length in seconds. Defines df = 1 / fftlength.

  • overlap (float or Quantity, optional) – Overlap between segments in seconds.

  • window (str or array_like, optional) – Window function name or samples.

  • method (str, optional) – PSD estimation method name registered in GWpy.

Returns:

One-sided PSD with unit (input unit)**2 / Hz.

Return type:

FrequencySeries