Matched Filtering (ligo.skymap.bayestar.filter)

Utility functions for BAYESTAR that are related to matched filtering.

class ligo.skymap.bayestar.filter.InterpolatedPSD(f, S, f_high_truncate=1.0, fill_value=inf)[source] [edit on github]

Bases: interp1d

Create a (linear in log-log) interpolating function for a discretely sampled power spectrum S(f).

Initialize a 1-D linear interpolation class.

class ligo.skymap.bayestar.filter.SignalModel(h)[source] [edit on github]

Bases: object

Class to speed up computation of signal/noise-weighted integrals and Barankin and Cramér-Rao lower bounds on time and phase estimation.

Note that the autocorrelation series and the moments are related, as shown below.


Create signal model:

>>> from . import filter
>>> sngl = lambda: None
>>> H = filter.sngl_inspiral_psd(
...     'TaylorF2threePointFivePN', mass1=1.4, mass2=1.4)
>>> S = vectorize_swig_psd_func('SimNoisePSDaLIGOZeroDetHighPower')
>>> W = filter.signal_psd_series(H, S)
>>> sm = SignalModel(W)

Compute one-sided autocorrelation function:

>>> out_duration = 0.1
>>> a, sample_rate = filter.autocorrelation(W, out_duration)

Restore negative time lags using symmetry:

>>> a = np.concatenate((a[:0:-1].conj(), a))

Compute the first 2 frequency moments by taking derivatives of the autocorrelation sequence using centered finite differences. The nth frequency moment should be given by (-1j)^n a^(n)(t).

>>> acor_moments = []
>>> for i in range(2):
...     acor_moments.append(a[len(a) // 2])
...     a = -0.5j * sample_rate * (a[2:] - a[:-2])
>>> assert np.all(np.isreal(acor_moments))
>>> acor_moments = np.real(acor_moments)

Compute the first 2 frequency moments using this class.

>>> quad_moments = [sm.get_sn_moment(i) for i in range(2)]

Compare them.

>>> for i, (am, qm) in enumerate(zip(acor_moments, quad_moments)):
...     assert np.allclose(am, qm, rtol=0.05)

Create a TaylorF2 signal model with the given masses, PSD function S(f), PN amplitude order, and low-frequency cutoff.

get_crb(snr)[source] [edit on github]

Get the Cramér-Rao bound, or inverse Fisher information matrix, describing the phase and time estimation covariance.

get_sn_average(func)[source] [edit on github]

Get the average of a function of angular frequency, weighted by the signal to noise per unit angular frequency.

get_sn_moment(power)[source] [edit on github]

Get the average of angular frequency to the given power, weighted by the signal to noise per unit frequency.

ligo.skymap.bayestar.filter.abs2(y)[source] [edit on github]

Return the absolute value squared, \(|z|^2\) ,for a complex number \(z\), without performing a square root.

ligo.skymap.bayestar.filter.abscissa(series)[source] [edit on github]

Produce the independent variable for a lal TimeSeries or FrequencySeries.

ligo.skymap.bayestar.filter.autocorrelation(H, out_duration, normalize=True)[source] [edit on github]

Calculate the complex autocorrelation sequence a(t), for t >= 0, of an inspiral signal.


Signal PSD series.


Noise power spectral density function.


The complex-valued autocorrelation sequence.


The sample rate.

ligo.skymap.bayestar.filter.ceil_pow_2(n)[source] [edit on github]

Return the least integer power of 2 that is greater than or equal to n.


>>> ceil_pow_2(128.0)
>>> ceil_pow_2(0.125)
>>> ceil_pow_2(129.0)
>>> ceil_pow_2(0.126)
>>> ceil_pow_2(1.0)
ligo.skymap.bayestar.filter.get_approximant_and_orders_from_string(s)[source] [edit on github]

Determine the approximant, amplitude order, and phase order for a string of the form “TaylorT4threePointFivePN”. In this example, the waveform is “TaylorT4” and the phase order is 7 (twice 3.5). If the input contains the substring “restricted” or “Restricted”, then the amplitude order is taken to be 0. Otherwise, the amplitude order is the same as the phase order.

ligo.skymap.bayestar.filter.get_f_lso(mass1, mass2)[source] [edit on github]

Calculate the GW frequency during the last stable orbit of a compact binary.

ligo.skymap.bayestar.filter.lal_ndebug()[source] [edit on github]

Temporarily disable lal error messages, except for memory errors.

ligo.skymap.bayestar.filter.truncated_ifft(y, nsamples_out=None)[source] [edit on github]

Truncated inverse FFT.

See for a discussion of related algorithms.

Perform inverse FFT to obtain truncated autocorrelation time series. This makes use of a folded DFT for a speedup of:


over directly computing the inverse FFT and truncating. Here is how it works. Say we have a frequency-domain signal X[k], for 0 ≤ k ≤ N - 1. We want to compute its DFT x[n], for 0 ≤ n ≤ M, where N is divisible by M: N = cM, for some integer c. The DFT is:

       N - 1
       \           2 π i k n
x[n] =  \     exp[-----------] Y[k]
       /               N
       k = 0

       c - 1   M - 1
       ______  ______
       \       \           2 π i n (m c + j)
     =  \       \     exp[------------------] Y[m c + j]
       /       /                 c M
      /------ /------
       j = 0   m = 0

       c - 1                     M - 1
       ______                    ______
       \           2 π i n j     \           2 π i n m
     =  \     exp[-----------]    \     exp[-----------] Y[m c + j]
       /               N         /               M
      /------                   /------
       j = 0                     m = 0

So: we split the frequency series into c deinterlaced sub-signals, each of length M, compute the DFT of each sub-signal, and add them back together with complex weights.


Complex input vector.

nsamples_outint, optional

Length of output vector. By default, same as length of input vector.


The first nsamples_out samples of the IFFT of x, zero-padded if


First generate the IFFT of a random signal:

>>> nsamples_out = 1024
>>> y = np.random.randn(nsamples_out) + np.random.randn(nsamples_out) * 1j
>>> x = fft.ifft(y)

Now check that the truncated IFFT agrees:

>>> np.allclose(x, truncated_ifft(y), rtol=1e-15)
>>> np.allclose(x, truncated_ifft(y, 1024), rtol=1e-15)
>>> np.allclose(x[:128], truncated_ifft(y, 128), rtol=1e-15)
>>> np.allclose(x[:1], truncated_ifft(y, 1), rtol=1e-15)
>>> np.allclose(x[:32], truncated_ifft(y, 32), rtol=1e-15)
>>> np.allclose(x[:63], truncated_ifft(y, 63), rtol=1e-15)
>>> np.allclose(x[:25], truncated_ifft(y, 25), rtol=1e-15)
>>> truncated_ifft(y, 1025)
Traceback (most recent call last):
ValueError: Input is too short: you gave me an input of length 1024, but you asked for an IFFT of length 1025.
ligo.skymap.bayestar.filter.unwrap(y, *args, **kwargs)[source] [edit on github]

Unwrap phases while skipping NaN or infinite values.

This is a simple wrapper around numpy.unwrap() that can handle invalid values.


>>> t = np.arange(0, 2 * np.pi, 0.5)
>>> y = np.exp(1j * t)
>>> unwrap(np.angle(y))
array([0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. ])
>>> y[3] = y[4] = y[7] = np.nan
>>> unwrap(np.angle(y))
array([0. , 0.5, 1. , nan, nan, 2.5, 3. , nan, 4. , 4.5, 5. , 5.5, 6. ])
class ligo.skymap.bayestar.filter.vectorize_swig_psd_func(str)[source] [edit on github]

Bases: object

Create a vectorized Numpy function from a SWIG-wrapped PSD function. SWIG does not provide enough information for Numpy to determine the number of input arguments, so we can’t just use np.vectorize.