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)

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)

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.

Examples

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)

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

get_sn_average(func)

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

get_sn_moment(power)

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)

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

ligo.skymap.bayestar.filter.abscissa(series)

Produce the independent variable for a lal TimeSeries or FrequencySeries.

ligo.skymap.bayestar.filter.autocorrelation(H, out_duration, normalize=True)

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

Parameters:

Hlal.REAL8FrequencySeries

Signal PSD series.

Scallable

Noise power spectral density function.

Returns:

acornumpy.ndarray

The complex-valued autocorrelation sequence.

sample_ratefloat

The sample rate.

ligo.skymap.bayestar.filter.ceil_pow_2(n)

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

Examples

>>> ceil_pow_2(128.0)
128.0
>>> ceil_pow_2(0.125)
0.125
>>> ceil_pow_2(129.0)
256.0
>>> ceil_pow_2(0.126)
0.25
>>> ceil_pow_2(1.0)
1.0

ligo.skymap.bayestar.filter.get_approximant_and_orders_from_string(s)

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)

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

ligo.skymap.bayestar.filter.lal_ndebug()

Temporarily disable lal error messages, except for memory errors.

ligo.skymap.bayestar.filter.truncated_ifft(y, nsamples_out=None)

Truncated inverse FFT.

See http://www.fftw.org/pruned.html 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:

log(nsamples)/log(nsamples_out)

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.

Parameters:

ynumpy.ndarray

Complex input vector.

nsamples_outint, optional

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

Returns:

xnumpy.ndarray

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

Examples

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)
True
>>> np.allclose(x, truncated_ifft(y, 1024), rtol=1e-15)
True
>>> np.allclose(x[:128], truncated_ifft(y, 128), rtol=1e-15)
True
>>> np.allclose(x[:1], truncated_ifft(y, 1), rtol=1e-15)
True
>>> np.allclose(x[:32], truncated_ifft(y, 32), rtol=1e-15)
True
>>> np.allclose(x[:63], truncated_ifft(y, 63), rtol=1e-15)
True
>>> np.allclose(x[:25], truncated_ifft(y, 25), rtol=1e-15)
True
>>> 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)

Unwrap phases while skipping NaN or infinite values.

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

Examples

>>> 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)

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.