gwadama.fat#

fat.py

Frequency analysis toolkit.

gwadama.fat.bandpass_filter(signal: ndarray, *, f_low: int | float, f_high: int | float, f_order: int | float, fs: int) → ndarray[source]#

Apply a forward-backward digital bandpass filter.

Apply a forward-backward digital bandpass filter to ‘signal’ between frequencies ‘f_low’ and ‘f_high’ with an order of ‘f_order’.

gwadama.fat.correct_phase(phase, time, jump_start, jump_end, correction_factor=1.0)[source]#

Manually correct a phase jump.

Fine-tunes the manual phase correction by adjusting the phase after a phase jump. The phase after the jump is scaled by the correction factor.

Parameters:
phasendarray

The unwrapped phase of the signal.

timendarray

The time array corresponding to the signal.

jump_startfloat

The time where the phase jump starts.

jump_endfloat

The time where the phase jump ends.

correction_factorfloat, optional

The factor to scale the phase correction for fine-tuning. Default is 1.0.

Returns:
corrected_phasendarray

The phase with fine-tuned manual correction applied after the specified jump.

gwadama.fat.find_power_excess(strain, fs, nperseg=16, noverlap=15, return_time=False)[source]#

Estimate the index of maximum broadband power excess in a 1D signal.

This function provides a fast and simple estimate of the time region where the total broadband power in the input signal is maximised. A spectrogram is computed via Welch’s method, and the time slice with the highest total power (summed over all frequencies) is identified.

Formally, at each time step t, the total power is estimated as:

P(t) = ∑_f S(f, t)

where S(f, t) is the estimated power spectral density at frequency f and time t.

The accuracy of this estimate is limited by the spectrogram parameters chosen by the user (nperseg, noverlap), which determine the temporal resolution and frequency leakage.

Parameters:
strainnp.ndarray

Input signal (1D array).

fsint

sampling frequency in Hz.

npersegint, optional

Length of each FFT segment (in samples).

noverlapint, optional

Number of overlapping samples between segments.

return_timebool, optional

If True, also return the estimated time (in seconds) corresponding to the peak index.

Returns:
peak_indexint

Sample index (in the input array) corresponding to the centre of the most power-concentrated region.

peak_timefloat, optional

Time in seconds of the estimated peak, assuming the first sample occurs at t=0. Only returned if return_time=True.

Notes

  • In real detector data, whitening or equalisation is recommended to avoid biasing the result towards dominant background frequencies.

  • The returned index refers to the centre of the peak time bin in the spectrogram, mapped back to the corresponding sample index of the original signal.

  • Spectrogram time resolution is limited by nperseg and noverlap; high accuracy is not guaranteed.

  • This method assumes only one dominant transient is present in the input. It is not suitable for detecting or resolving multiple events.

gwadama.fat.highpass_filter(signal: ndarray, *, f_cut: int | float, f_order: int | float, fs: int) → ndarray[source]#

Apply a forward-backward digital highpass filter.

Apply a forward-backward digital highpass filter to ‘signal’ at frequency ‘f_cut’ with an order of ‘f_order’.

gwadama.fat.instant_frequency(signal, *, fs, phase_corrections=None)[source]#

Computes the instantaneous frequency of a time-domain signal.

Computes the instantaneous frequency of a time-domain signal using the central difference method, with optional phase corrections. If negative frequencies are detected, a warning is raised.

Parameters:
signalndarray

The input time-domain signal.

fsfloat

The sampling frequency of the signal (in Hz).

phase_correctionslist of tuples, optional

A list of phase corrections. Each tuple contains: (jump_start, jump_end, correction_factor). If None, no phase correction is applied.

Returns:
inst_freqndarray

The instantaneous frequency of the signal (in Hz). If negative frequencies are detected, a RuntimeWarning is issued.

gwadama.fat.snr(strain, *, at, psd=None, window=('tukey', 0.5), fbounds=None)[source]#

Compute the signal-to-noise ratio (SNR) of a signal segment.

Compute the SNR of a signal with respect to a background spectrum (psd provided), or assuming the signal is already whitened (psd=None).

Parameters:
strainarray_like

Time-domain signal segment.

atfloat

Sampling interval (seconds per sample).

psdarray_like (freq_array, psd_array), optional

One-sided power spectral density of the background noise. If provided, the SNR is computed using this PSD. If None, the function assumes the signal is already whitened.

windowstr, tuple, or array_like, optional

Window to apply before the FFT. Passed to scipy.signal.get_window if a string/tuple. If an array, it must match the length of strain.

fboundsarray_like (f_min, f_max), optional

Frequency range (Hz) over which to compute the SNR. If None: * with psd, uses the min and max of the PSD frequencies, * otherwise uses (0, Nyquist).

Returns:
snrfloat

The computed signal-to-noise ratio.

Notes

  • The function uses a one-sided rFFT and sums over positive frequencies.

  • In whitened space, PSD is assumed to be identity and the formula reduces to: sqrt(2 * Δt * Δf * Σ |h_w(f)|²).

  • In non-whitened space, the formula used is: sqrt(4 * Δt² * Δf * Σ |h(f)|² / S_n(f)).