import numpy as np
from scipy.signal import butter, lfilter
[docs]
def whitened_snr_scaling(glitch, snr, srate=4096):
"""Scale a glitch signal to the target SNR in the whitened frame."""
glitch = np.asarray(glitch)
if snr is not None:
df = srate / glitch.shape[-1]
glitch_FD = np.fft.rfft(glitch, axis=-1) / srate
true_sigma_sq = 4.0 * df * np.sum(np.multiply(np.conj(glitch_FD), glitch_FD), axis=-1).real
glitch = (glitch.T * snr / np.sqrt(true_sigma_sq)).T
return glitch
[docs]
def quality_factor_conversion(Q, f_0):
"""Convert quality factor Q and central frequency f_0 to decay time tau."""
tau = Q / (np.sqrt(2) * np.pi * f_0)
return tau
[docs]
def rescale(x):
"""Rescale each row of x to the range [-1, 1]."""
abs_max = np.max(np.abs(x), axis=1, keepdims=True)
return x / abs_max
[docs]
def butter_lowpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype="low", analog=False)
return b, a
[docs]
def butter_highpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
d, c = butter(order, normal_cutoff, btype="high", analog=False)
return d, c
[docs]
def butter_filter(data, fs, order=5):
"""Apply a bandpass (20–1024 Hz) Butterworth filter to data."""
low_cutoff = 1024
high_cutoff = 20
b, a = butter_lowpass(low_cutoff, fs, order=order)
low_filtered = lfilter(b, a, data)
d, c = butter_highpass(high_cutoff, fs, order=order)
return lfilter(d, c, low_filtered)
[docs]
def custom_whiten(self, psd, low_frequency_cutoff=None, return_psd=False, **kwds):
"""
Return a whitened PyCBC TimeSeries.
This function is designed to be used with a PyCBC TimeSeries instance
(as a monkey-patched method). Pass ``self`` as the TimeSeries object.
Parameters
----------
psd : FrequencySeries
The power spectral density used for whitening.
low_frequency_cutoff : float, optional
Low frequency cutoff for the inverse spectrum truncation.
return_psd : bool, optional
If True, return the PSD alongside the whitened data.
Returns
-------
white : TimeSeries
The whitened time series.
psd : FrequencySeries, optional
The PSD used (only returned if ``return_psd=True``).
"""
white = (self.to_frequencyseries() / psd**0.5).to_timeseries()
if return_psd:
return white, psd
return white
[docs]
def generate_gaussian_noise(mean, std_dev, num_samples, sample_shape):
"""Generate Gaussian noise samples as a numpy array."""
return np.random.normal(loc=mean, scale=std_dev, size=(num_samples, *sample_shape))