ml4gw.spectral

This module provides functions for calculation of spectral densities and for whitening.

Several implementation details are derived from the scipy csd and welch implementations. For more info, see

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html

and

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.csd.html

Functions

`fast_spectral_density`(x, nperseg, nstride, ...)	Compute the power spectral density of a multichannel timeseries or a batch of multichannel timeseries, or the cross power spectral density of two such timeseries.
`median`(x, axis)	Implements a median calculation that matches numpy's behavior for an even number of elements and includes the same bias correction used by scipy's implementation.
`minimum_phase_whitening_filter`(psd[, n_fft, ...])	Construct a causal, minimum-phase whitening filter.
`normalize_by_psd`(X, psd, sample_rate, pad[, ...])
`spectral_density`(x, nperseg, nstride, ...[, ...])	Compute the power spectral density of a multichannel timeseries or a batch of multichannel timeseries.
`truncate_inverse_power_spectrum`(psd, ...[, ...])	Truncate the length of the time domain response of a whitening filter built using the specified `psd` so that it has maximum length `fduration` seconds.
`whiten`(X, psd, fduration, sample_rate[, ...])	Whiten a batch of timeseries using the specified background one-sided power spectral densities (PSDs), modified to have the desired time domain response length `fduration` and possibly to highpass/lowpass filter.

ml4gw.spectral.fast_spectral_density(x, nperseg, nstride, window, scale, average='median', y=None)

Compute the power spectral density of a multichannel timeseries or a batch of multichannel timeseries, or the cross power spectral density of two such timeseries. This implementation is non-exact for the two lowest frequency bins, since it implements centering using the mean for the entire timeseries rather than on a per- window basis. The benefit of this is a faster implementation, which might be beneficial for cases where the lowest frequency components are going to be discarded anyway.

Parameters:

x (Float[Tensor, 'time'] | Float[Tensor, 'channel time'] | Float[Tensor, 'batch channel time']) -- The timeseries tensor whose power spectral density to compute, or for cross spectral density the timeseries whose fft will be conjugated. Can have shape (batch_size, num_channels, length * sample_rate), (num_channels, length * sample_rate), or (length * sample_rate).
nperseg (int) -- Number of samples included in each FFT window
nstride (int) -- Stride between FFT windows
window (Float[Tensor, '{nperseg//2+1}']) -- Window array to multiply by each FFT window before FFT computation. Should have length nperseg // 2 + 1.
scale (float) -- Scale factor to multiply the FFT'd data by, related to desired units for output tensor (e.g. letting this equal 1 / (sample_rate * (window**2).sum()) will give output units of density, :math``text{Hz}^-1``.
average (str) -- How to aggregate the contributions of each FFT window to the spectral density. Allowed options are 'mean' and 'median'.
y (Float[Tensor, 'time'] | Float[Tensor, 'channel time'] | Float[Tensor, 'batch channel time'] | None) -- Timeseries tensor to compute cross spectral density with x. If left as None, x's power spectral density will be returned. Otherwise, if x is 1D, y must also be 1D. If x is 2D, the assumption is that this represents a single multi-channel timeseries, and y must be either 2D or 1D. In the former case, the cross-spectral densities of each channel will be computed individually, so y must have the same shape as x. Otherwise, this will compute the CSD of each of x's channels with y. If x is 3D, this will be assumed to be a batch of multi-channel timeseries. In this case, y can either be 3D, in which case each channel of each batch element will have its CSD calculated or 2D, which has two different options. If y's 0th dimension matches x's 0th dimension, it will be assumed that y represents a batch of 1D timeseries, and for each batch element this timeseries will have its CSD with each channel of the corresponding batch element of x calculated. Otherwise, it sill be assumed that y represents a single multi-channel timeseries, in which case each channel of y will have its CSD calculated with the corresponding channel in x across _all_ of x's batch elements.

Return type:

Float[Tensor, 'frequency'] | Float[Tensor, 'channel frequency'] | Float[Tensor, 'batch channel frequency']

Returns:

Tensor of power spectral densities of x or its cross spectral density with the timeseries in y.

ml4gw.spectral.median(x, axis)

Implements a median calculation that matches numpy's behavior for an even number of elements and includes the same bias correction used by scipy's implementation.

Return type:

']

Parameters:

x (Float[Tensor, '... size'])
axis (int)

ml4gw.spectral.minimum_phase_whitening_filter(psd, n_fft=None, validate=True)

Construct a causal, minimum-phase whitening filter.

The filter is obtained by spectral factorization of the inverse power spectral density using a folded real cepstrum. Its one-sided frequency response has magnitude 1 / sqrt(psd), while its phase minimizes the delay of the corresponding impulse response.

Parameters:

psd (Float[Tensor, 'num_ifos frequency']) -- A strictly positive, finite one-sided power spectral density. Any leading dimensions are treated as batch dimensions. Use double precision for unscaled physical PSDs whose values may fall below the representable range of torch.float32.
n_fft (int | None) -- Length of the desired impulse response. If omitted, an even length of 2 * (psd.size(-1) - 1) is inferred.
validate (bool) -- If True, verify that all PSD values are finite and positive. Disable only when these properties are guaranteed by the caller, since validating a GPU tensor synchronizes with the host.

Return type:

Tensor

Returns:

A real-valued impulse response with the same leading dimensions as: psd and length n_fft.

Raises:

ValueError -- If the PSD shape is inconsistent with n_fft, or, when validate is True, if it contains a non-finite or non-positive value.

Notes

The folded-cepstrum construction follows the zlw MPWhiteningFilter implementation at commit b81871cd.

ml4gw.spectral.normalize_by_psd(X, psd, sample_rate, pad, crop=True)

Parameters:

X (Float[Tensor, 'batch num_ifos time'])
psd (Float[Tensor, 'num_ifos frequency'])
sample_rate (float)
pad (int)
crop (bool)

ml4gw.spectral.spectral_density(x, nperseg, nstride, window, scale, average='median')

Compute the power spectral density of a multichannel timeseries or a batch of multichannel timeseries. This implementation is exact for all frequency bins, but slower than the fast implementation.

Parameters:

x (Float[Tensor, 'time'] | Float[Tensor, 'channel time'] | Float[Tensor, 'batch channel time']) -- The timeseries tensor whose power spectral density to compute, or for cross spectral density the timeseries whose fft will be conjugated. Can have shape (batch_size, num_channels, length * sample_rate), (num_channels, length * sample_rate), or (length * sample_rate).
nperseg (int) -- Number of samples included in each FFT window
nstride (int) -- Stride between FFT windows
window (Float[Tensor, '{nperseg//2+1}']) -- Window array to multiply by each FFT window before FFT computation. Should have length nperseg // 2 + 1.
scale (float) -- Scale factor to multiply the FFT'd data by, related to desired units for output tensor (e.g. letting this equal 1 / (sample_rate * (window**2).sum()) will give output units of density, \(\text{Hz}^-1\).
average (str) -- How to aggregate the contributions of each FFT window to the spectral density. Allowed options are 'mean' and 'median'.

Return type:

Float[Tensor, 'frequency'] | Float[Tensor, 'channel frequency'] | Float[Tensor, 'batch channel frequency']

ml4gw.spectral.truncate_inverse_power_spectrum(psd, fduration, sample_rate, highpass=None, lowpass=None)

Truncate the length of the time domain response of a whitening filter built using the specified psd so that it has maximum length fduration seconds. This is meant to mitigate the impact of sharp features in the background PSD causing time domain responses longer than the segments to which the whitening filter will be applied.

Implementation details adapted from here.

Parameters:

psd (Float[Tensor, 'num_ifos frequency']) -- The one-sided power spectral density used to construct a whitening filter.
fduration (Float[Tensor, 'time'] | float) -- Desired length in seconds of the time domain response of a whitening filter built using this PSD, or a window of this length to taper the edges of the time domain response of the filter. If passed as a float, a Hann window of this length will be used.
sample_rate (float) -- Rate at which the time domain data to which the whitening filter will be applied has been sampled.
highpass (float | None) -- If specified, will zero out the frequency response of all frequencies below this value in Hz. If left as None, no highpass filtering will be applied.
lowpass (float | None) -- If specified, will zero out the frequency response of all frequencies above this value in Hz. If left as None, no lowpass filtering will be applied.

Return type:

Float[Tensor, 'num_ifos frequency']

Returns:

The PSD with its time domain response truncated: to fduration and any filtered frequencies tapered.

ml4gw.spectral.whiten(X, psd, fduration, sample_rate, highpass=None, lowpass=None, crop=True)

Whiten a batch of timeseries using the specified background one-sided power spectral densities (PSDs), modified to have the desired time domain response length fduration and possibly to highpass/lowpass filter.

Parameters:

X (Float[Tensor, 'batch num_ifos time']) -- batch of multichannel timeseries to whiten
psd (Float[Tensor, 'num_ifos frequency']) -- PSDs use to whiten the data. The frequency response of the whitening filter will be roughly the inverse of the square root of this PSD, ensuring that data from the same distribution will have approximately uniform power after whitening. If 2D, each batch element in X will be whitened using the same PSDs. If 3D, each batch element will be whitened by the PSDs contained along the 0th dimenion of psd, and so the first two dimensions of X and psd should match.
fduration (Float[Tensor, 'time'] | float) -- Desired length in seconds of the time domain response of a whitening filter built using this PSD, or a window of this length to taper the edges of the time domain response of the filter. If passed as a float, a Hann window of this length will be used. Moreover, half of this length will be removed from each edge of the whitened timeseries to account for filter settle-in time.
sample_rate (float) -- Rate at which the data in X has been sampled
highpass (float | None) -- The frequency in Hz at which to highpass filter the data, setting the frequency response in the whitening filter to 0. If left as None, no highpass filtering will be applied.
lowpass (float | None) -- The frequency in Hz at which to lowpass filter the data, setting the frequency response in the whitening filter to 0. If left as None, no lowpass filtering will be applied.
crop (bool) -- If True, crop fduration / 2 seconds of data from both sides of the time dimension to remove the corruption from the filter. If False, return the full timeseries.

Return type:

Float[Tensor, 'batch num_ifos time']

Returns:

Batch of whitened multichannel timeseries with fduration / 2: seconds optionally trimmed from each side.