See it work

Every figure below was produced by running the library while this page was built. The code that produced them is at the bottom, taken from the module that ran -- not transcribed, so the two cannot disagree.

7/7

peaks agree with numpy

7/7

behaved as expected

cases within 1.0 of the threshold

What you are looking at

Each panel is one (data, template) pair: white noise of 512 samples, transformed, correlated against a band-limited template. The blue curve is the full correlation computed by numpy -- all 512 lags, the answer the library must agree with. The dashed line is the threshold the caller asked for. The circle is the single sample the library reported.

The library returns only that circle. Computing the blue curve is the work it is allowed to skip, which is the entire reason it is faster, so the curve is here as the check rather than as output.

case

What each case shows

case	injected at	loudest sample	reported lag	reported snr	agrees with numpy
pure noise	-	3.26	nothing	-	yes
pure noise	-	2.69	nothing	-	yes
snr 4.2, well below	137	4.38	nothing	-	yes
snr 4.8, just below	64	3.45	nothing	-	yes
snr 5.3, just above	301	5.80	301	5.80	yes
snr 6.0, clear	200	6.55	200	6.55	yes
snr 9.0, loud	411	9.05	411	9.05	yes

The cases near the threshold are the ones worth reading. A signal injected at snr 4.2 does not arrive measuring 4.2: the noise it lands in moves it, by about a unit either way. So an injection below the threshold can clear it and one above can fail to, and both happen here. What the filter is responsible for is narrower and is what these plots check -- reporting the loudest sample, and reporting it only when it crosses. It agrees with numpy in every case, including the ones where the decision is close.

Concretely: injected at snr 4.2, loudest sample 4.38, below the threshold of 5.0, so silent; injected at snr 4.8, loudest sample 3.45, below the threshold of 5.0, so silent. That is the detection statistics, not a fault in the filter -- the flat filter reports nothing there either, which is exactly what the blue curve shows. It is also why the hierarchical mode is tuned against measured false-dismissal rather than against a model: near the threshold is where a cheap first pass could lose something, so that is where it has to be measured.

The pure-noise cases report nothing at all. That is the filter working: the loudest noise sample never reached the threshold, so there was no peak to report, and the library says so rather than handing back its largest fluctuation.

The code

This is the source of the functions that ran, read from the module at build time.

def make_template(n, rng):
    """A unit-norm, band-limited template spectrum.

    Band-limited because a real search's templates are, and unit-norm so that
    the filter's output reads directly as a signal-to-noise ratio.
    """
    band = n // 4
    amp = np.zeros(n)
    amp[:band] = np.exp(-np.arange(band) / (band / 3.0))
    h = (amp * np.exp(2j * np.pi * rng.random(n))).astype(np.complex64)
    return h / np.linalg.norm(h)

def make_data(n, h, rng, snr=0.0, lag=0):
    """White noise, optionally with one copy of the template buried in it.

    The filter works on spectra, so this returns a spectrum. Multiplying by
    the phase ramp is exactly a circular shift of `lag` samples in time, which
    is how the signal is placed without leaving the frequency domain.
    """
    d = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
    if snr:
        ramp = np.exp(-2j * np.pi * lag * np.arange(n) / n)
        d = d + (snr * h * ramp).astype(np.complex64)
    return d

def correlate(d, h):
    """The full correlation, by numpy, as the reference to check against.

    This is the whole answer: n samples, one per lag. The library computes the
    same thing but reports only the loudest sample per bin, which is what
    makes it fast -- so this is what it must agree with.
    """
    return np.abs(np.fft.ifft(d * np.conj(h)) * len(d))

def run_case(n, rng, snr, lag, threshold):
    """Filter one (data, template) pair and check the peak against numpy."""
    h = make_template(n, rng)
    d = make_data(n, h, rng, snr=snr, lag=lag)

    filt = mf.MatchedFilter(n, ndata=1, ntemplates=1)
    filt.set_data(d[None, :])
    filt.set_templates(h[None, :])
    peak = filt.run(binsize=n, threshold=threshold)[0, 0, 0]

    rho = correlate(d, h)
    return peak, rho

Run it yourself with python -m matchedfilter.demo.

← Using it

Numerical accuracy →