Using matchedfilter

This guide covers inputs, full correlations, output bins, device selection and hierarchical calibration. The examples execute during the documentation build. Run them locally with python -m matchedfilter.tutorial.

Install

pip install matchedfilter

The current release is alpha. Pin its version for reproducible work. Wheels target CPython 3.10–3.14 on Linux x86-64 and macOS arm64. Source installation requires NumPy and a C compiler.

Inputs and normalization

set_data() and set_templates() accept frequency-domain arrays in natural order, with shapes (ndata, n) and (ntemplates, n). Use complex64 spectra from an unnormalized forward FFT, such as numpy.fft.fft. The filter computes the unnormalized inverse transform of data * conj(template).

Declare the dimensions when constructing the filter. Plans reuse storage; first-use GPU dispatches and adaptive CPU layouts can still allocate.

A complete example

import numpy as np
import matchedfilter as mf

n, ndata, ntemplates = 4096, 8, 32
rng = np.random.default_rng(0)

# Inputs are SPECTRA: the unnormalised forward transform of each segment.
# With unit-norm templates and unit-variance noise, the reported
# magnitude reads directly as a signal-to-noise ratio.
templates = (rng.standard_normal((ntemplates, n))
             + 1j * rng.standard_normal((ntemplates, n))).astype(np.complex64)
templates /= np.linalg.norm(templates, axis=1, keepdims=True)
data = (rng.standard_normal((ndata, n))
        + 1j * rng.standard_normal((ndata, n))).astype(np.complex64)

# Bury one copy of template 3 in segment 5, at lag 900.
ramp = np.exp(-2j * np.pi * 900 * np.arange(n) / n)
data[5] += (9.0 * templates[3] * ramp).astype(np.complex64)

filt = mf.MatchedFilter(n, ndata=ndata, ntemplates=ntemplates)
filt.set_data(data)
filt.set_templates(templates)
peaks = filt.run(binsize=n, threshold=6.0)

print("peaks.shape", peaks.shape, " fields", peaks.dtype.names)
found = np.argwhere(peaks["index"] >= 0)
for d, t, b in found:
    print("segment %d x template %d: lag %d, snr %.2f"
          % (d, t, peaks["index"][d, t, b], np.abs(peaks["value"])[d, t, b]))
print("everything else is below threshold and reports index -1")

output

peaks.shape (8, 32, 1)  fields ('index', 'value')
segment 5 x template 3: lag 900, snr 9.07
everything else is below threshold and reports index -1

For unit-norm template spectra and independent noise with unit variance in each real and imaginary component, the output magnitude has the normalization used in these SNR examples. Other input normalizations change the scale of threshold. Check the normalization of your data and template bank.

What comes back

import numpy as np
import matchedfilter as mf

n = 1024
rng = np.random.default_rng(4)
h = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
h /= np.linalg.norm(h)
d = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
d += (8.0 * h * np.exp(-2j * np.pi * 300 * np.arange(n) / n)).astype(np.complex64)

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

v, m = complex(pk["value"]), float(np.abs(pk["value"]))
print("index      %d" % int(pk["index"]))
print("value      %+.4f%+.4fj      the complex sample, so phase is available"
      % (v.real, v.imag))
print("magnitude  %.6f" % m)
print("abs(value) %.6f   differs by %.1e" % (abs(v), abs(abs(v) - m)))
print()
print("magnitude is the number the threshold was compared against, so it")
print("is the one to re-test against; recomputing abs(value) can land a")
print("few ULPs the other side of a threshold.")

output

index      300
value      +7.1331+1.0169j      the complex sample, so phase is available
magnitude  7.205201
abs(value) 7.205201   differs by 2.9e-07

magnitude is the number the threshold was compared against, so it
is the one to re-test against; recomputing abs(value) can land a
few ULPs the other side of a threshold.

Output bins and thresholds

run() returns an array with shape (ndata, ntemplates, nbins):

Field	Type	Meaning
`index`	int64	lag of the strongest sample in the bin; −1 if dismissed
`value`	complex64	complex correlation at that lag; zero if dismissed

binsize is the number of lags in each output bin. A full window with binsize=n returns one peak per pair. A partial final bin is included.

One peak per window

import numpy as np
import matchedfilter as mf

n = 4096
rng = np.random.default_rng(1)
h = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
h /= np.linalg.norm(h)
d = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
d += (7.0 * h * np.exp(-2j * np.pi * 2500 * np.arange(n) / n)).astype(np.complex64)

filt = mf.MatchedFilter(n, 1, 1)
filt.set_data(d[None, :])
filt.set_templates(h[None, :])

for bs in (4096, 1024, 256):
    pk = filt.run(binsize=bs, threshold=0.0)
    lags = [int(i) for i in pk["index"][0, 0]]
    signal = next(j for j, i in enumerate(lags) if i == 2500)
    print("binsize %5d -> %3d peaks, lags %-34s signal in bin %d"
          % (bs, pk.shape[2], str(lags[:4]) + (" ..." if len(lags) > 4 else ""),
             signal))
print()
print("Every bin reports its own loudest lag, so binsize is the time")
print("resolution of the output. A search that clusters candidates at a")
print("fixed resolution sets binsize to it and does no clustering of its")
print("own. binsize=n is the other extreme: one peak for the whole pair.")

output

binsize  4096 ->   1 peaks, lags [2500]                             signal in bin 0
binsize  1024 ->   4 peaks, lags [425, 1236, 2500, 4007]            signal in bin 2
binsize   256 ->  16 peaks, lags [230, 425, 680, 772] ...           signal in bin 9

Every bin reports its own loudest lag, so binsize is the time
resolution of the output. A search that clusters candidates at a
fixed resolution sets binsize to it and does no clustering of its
own. binsize=n is the other extreme: one peak for the whole pair.

threshold applies to the peak magnitude. Dismissed bins keep their place in the output, so peaks[d, t, j] always refers to bin j.

Keeping every lag

CorrelationFilter accepts the same spectra and bank dimensions as MatchedFilter. Its run() returns a complex64 array shaped (ndata, ntemplates, n) in natural lag order. It computes the unnormalised inverse of data * conj(template), with no threshold or peak search:

import numpy as np
import matchedfilter as mf

n = 1024
template = np.zeros(n, np.float32)
template[0] = 1
data = np.roll(template, 37)
full = mf.CorrelationFilter(n)
full.set_data(np.fft.fft(data).astype(np.complex64)[None, :])
full.set_templates(np.fft.fft(template).astype(np.complex64)[None, :])
values = full.run()                  # shape (1, 1, 1024)
print(np.argmax(np.abs(values[0, 0])))  # 37

For larger banks, run(data=(start, count), templates=(start, count)) selects a rectangular subrange.

For overlap-save output, construct the filter with valid=(lo, hi), the half-open lag interval valid in each FFT block. run_series(series) then uses starts 0, hi-lo, 2*(hi-lo), ..., zero-pads the final input block, and returns one continuous complex64 array shaped (selected_templates, len(series)). The first lo samples start at zero; subsequent calls overwrite every valid sample and leave that invalid prefix untouched. The filter owns and reuses the result, so copy it if it must outlive the next call. CPU and GPU write valid lags into this result during correlation, without assembling a block-output cube. The GPU result uses host-cached shared storage for direct GPU writes and practical NumPy access.

Use run_blocks(series, starts, templates=None, out=None) for explicit block layouts. It gathers and zero-pads the named blocks and returns all lags as (blocks, selected_templates, n). For that form, a writable, C-contiguous complex64 out can reuse caller storage. Both forms consume the data slots, so a later run() needs another set_data() call. The older run_series(series, starts, ...) spelling remains supported.

Full results can be large. A single 2^22 correlation is 32 MiB; a 128×512 bank at that length is 2 TiB. Without out, the class raises before allocating more than 512 MiB. Select a bank subrange or provide storage to process larger results in batches.

Thresholding

import numpy as np
import matchedfilter as mf

n = 4096
rng = np.random.default_rng(2)
h = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
h /= np.linalg.norm(h)
d = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
d += (6.5 * h * np.exp(-2j * np.pi * 1000 * np.arange(n) / n)).astype(np.complex64)

filt = mf.MatchedFilter(n, 1, 1)
filt.set_data(d[None, :])
filt.set_templates(h[None, :])

for thr in (0.0, 4.0, 5.0, 5.5, 6.0):
    pk = filt.run(binsize=256, threshold=thr)
    kept = int((pk["index"][0, 0] >= 0).sum())
    loud = float(np.abs(pk["value"])[0, 0].max())
    print("threshold %.1f -> %2d of %d bins reported%s"
          % (thr, kept, pk.shape[2],
             (", loudest snr %.2f" % loud) if kept else ""))
print()
print("The injection was at snr 6.5 and measures 5.60: noise subtracts")
print("from a peak as often as it adds. A threshold set at the value you")
print("injected will lose about half of them.")
print("Bins below the threshold still occupy their slot, with index -1,")
print("so peaks[d, t, j] is always bin j and needs no searching.")

output

threshold 0.0 -> 16 of 16 bins reported, loudest snr 5.60
threshold 4.0 ->  2 of 16 bins reported, loudest snr 5.60
threshold 5.0 ->  1 of 16 bins reported, loudest snr 5.60
threshold 5.5 ->  1 of 16 bins reported, loudest snr 5.60
threshold 6.0 ->  0 of 16 bins reported

The injection was at snr 6.5 and measures 5.60: noise subtracts
from a peak as often as it adds. A threshold set at the value you
injected will lose about half of them.
Bins below the threshold still occupy their slot, with index -1,
so peaks[d, t, j] is always bin j and needs no searching.

window=(start, end) restricts the searched lags to [start, end). Use it to exclude the invalid wrap-around region when processing overlapping blocks.

Bounding the lags searched

import numpy as np
import matchedfilter as mf

n = 4096
rng = np.random.default_rng(3)
h = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
h /= np.linalg.norm(h)
d = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
# two signals: one inside the window we will ask for, one outside it
for lag, snr in ((500, 12.0), (2000, 8.0)):
    d += (snr * h * np.exp(-2j * np.pi * lag * np.arange(n) / n)).astype(np.complex64)

filt = mf.MatchedFilter(n, 1, 1)
filt.set_data(d[None, :])
filt.set_templates(h[None, :])

# .copy() is not optional here: run() hands back a buffer it reuses, so
# without it the second call would silently rewrite the first answer.
whole = filt.run(binsize=n, threshold=5.0).copy()
part = filt.run(binsize=n, threshold=5.0, window=(1024, 3072))
print("all %d lags     -> lag %d, snr %.1f"
      % (n, int(whole["index"][0, 0, 0]), np.abs(whole["value"])[0, 0, 0]))
print("lags 1024..3072 -> lag %d, snr %.1f"
      % (int(part["index"][0, 0, 0]), np.abs(part["value"])[0, 0, 0]))
print()
print("The louder signal at lag 500 is outside the window and is not")
print("reported. An overlap-save caller passes its valid span here so the")
print("wrap-around region is never searched -- and for the hierarchical")
print("filter that matters twice, because every extra lag is another")
print("chance for noise to force a full correlation.")

output

all 4096 lags     -> lag 500, snr 11.9
lags 1024..3072 -> lag 2000, snr 7.1

The louder signal at lag 500 is outside the window and is not
reported. An overlap-save caller passes its valid span here so the
wrap-around region is never searched -- and for the hierarchical
filter that matters twice, because every extra lag is another
chance for noise to force a full correlation.

Output buffer lifetime

Results can reuse internal storage. Call .copy() to retain them after the next filtering call. raw=True returns separate index and value arrays.

The output buffer is reused

import numpy as np
import matchedfilter as mf

n = 1024
rng = np.random.default_rng(6)
h = (rng.standard_normal(n) + 1j * rng.standard_normal(n)).astype(np.complex64)
h /= np.linalg.norm(h)
d = (rng.standard_normal((2, n)) + 1j * rng.standard_normal((2, n))).astype(np.complex64)
d[0] += (9.0 * h * np.exp(-2j * np.pi * 100 * np.arange(n) / n)).astype(np.complex64)
d[1] += (9.0 * h * np.exp(-2j * np.pi * 700 * np.arange(n) / n)).astype(np.complex64)

filt = mf.MatchedFilter(n, 2, 1)
filt.set_templates(h[None, :])

filt.set_data(d)
first = filt.run(binsize=n, threshold=5.0)
kept = first.copy()
filt.set_data(d[::-1].copy())
second = filt.run(binsize=n, threshold=5.0)

print("same object?      ", first is second, "  shares memory?",
      np.shares_memory(first, second))
print("first  (now)      ", [int(i) for i in first["index"][:, 0, 0]])
print("copy taken earlier", [int(i) for i in kept["index"][:, 0, 0]])
print()
print("Six allocations per call are nothing beside a 2^20 transform, but")
print("a caller driving small batches pays them every time -- at 37")
print("templates they were 15 of the 21 microseconds a call took. So the")
print("buffer is kept and copying is the caller's job. Anything that has")
print("to outlive the next run needs .copy().")

output

same object?       True   shares memory? True
first  (now)       [700, 100]
copy taken earlier [100, 700]

Six allocations per call are nothing beside a 2^20 transform, but
a caller driving small batches pays them every time -- at 37
templates they were 15 of the 21 microseconds a call took. So the
buffer is kept and copying is the caller's job. Anything that has
to outlive the next run needs .copy().

Filtering a time series

Flat and hierarchical peak filters use the same automatic layout: MatchedFilter(n, valid=(lo, hi)) or HierarchicalFilter(n, valid=(lo, hi)), then run_series(series, ...). Each block searches only its valid lag window; the final window clips at the end of the series. The result remains block-major peaks, with unused bins in the clipped final block dismissed (index=-1). Peak-only output has no continuous sample series to assemble. Automatic peak indices are absolute positions in the supplied series; explicit-block calls retain their block-local lag indices.

Use run_blocks(series, starts, win_start, win_end, ...) for irregular block layouts and per-block windows. starts gives each block's input position; win_start and win_end give its lag window. Both forms gather and pad blocks, compute forward FFTs, and filter on the selected device. The older run_series(series, starts, win_start, win_end, ...) spelling remains supported with the same block-local indices.

Templates must be set first. A later run() requires another set_data() call because series execution reuses the data slots. GPU execution is synchronous. See the source documentation for GPU FFTs and shared storage for memory and interoperability details.

Hierarchical filtering

HierarchicalFilter screens pairs using a coarse frequency band and refines candidates with the full filter. Its output has the same format as run(); screening can omit detections. snr and fd specify the signal strength and requested false-dismissal budget used for calibration.

set_reference(power) supplies the expected output-power spectrum, one value per frequency bin. Include the effect of the data's noise spectrum; a template's power alone is generally insufficient for colored noise.

The hierarchical mode

import numpy as np
import matchedfilter as mf

n, ndata, ntemplates = 4096, 8, 32
rng = np.random.default_rng(5)

# The expected power of the filter OUTPUT, bin by bin -- not the
# template's own power. The two differ whenever the data is coloured.
k = np.arange(1, n // 2)
reference = np.zeros(n, np.float32)
reference[1:n // 2] = k ** (-7.0 / 3.0) / ((0.015 * n / k) ** 4 + 1.0)
reference /= reference.sum()

amp = np.sqrt(reference)
templates = (amp * np.exp(2j * np.pi * rng.random((ntemplates, n)))).astype(np.complex64)
templates /= np.linalg.norm(templates, axis=1, keepdims=True)
data = (rng.standard_normal((ndata, n))
        + 1j * rng.standard_normal((ndata, n))).astype(np.complex64)
ramp = np.exp(-2j * np.pi * 1700 * np.arange(n) / n)
data[2] += (11.0 * templates[7] * ramp).astype(np.complex64)

hf = mf.HierarchicalFilter(n, ndata=ndata, ntemplates=ntemplates,
                           snr=6.0, fd=1e-3)
hf.set_reference(reference)          # drives the whole configuration
hf.set_templates(templates)
hf.set_data(data)
peaks = hf.run(binsize=n, threshold=6.0)

band, taps = hf.config
print("it chose band %d, taps %d" % (band, taps))
print("escalated to the full correlation: %.1f%% of pairs"
      % (100 * hf.refine_rate))
for d_, t_, b_ in np.argwhere(peaks["index"] >= 0):
    print("found segment %d x template %d: lag %d, snr %.2f"
          % (d_, t_, peaks["index"][d_, t_, b_],
             np.abs(peaks["value"])[d_, t_, b_]))
print()
print("Peaks it reports are identical to the flat filter's. It can omit,")
print("never invent, and fd is the budget for how often it may omit one.")

output

it chose band 512, taps 8
escalated to the full correlation: 0.0% of pairs
found segment 2 x template 7: lag 1700, snr 11.65

Peaks it reports are identical to the flat filter's. It can omit,
never invent, and fd is the budget for how often it may omit one.

Automatic configuration requires a covering measured calibration file. Missing coverage raises an error. Alternatively, supply band when constructing the filter and call set_coarse_threshold(value) before execution. Both explicit parameters are required; no model or default threshold substitutes for them. An explicit threshold carries no measured false-dismissal guarantee.

When it refuses

import numpy as np
import matchedfilter as mf

n = 4096
reference = np.zeros(n, np.float32)
reference[1:n // 2] = 1.0
hf = mf.HierarchicalFilter(n, 1, 2, snr=5.0, fd=1e-6)
hf.set_reference(reference)
try:
    hf.run(binsize=n, threshold=5.0)
except ValueError as e:
    # Not every ValueError from run() carries the " -- " the tuning
    # refusal uses. Indexing [1] blindly turned a perfectly clear
    # error into an IndexError raised by the example itself, which
    # is a worse failure than the one being demonstrated.
    msg = str(e)
    head = msg.split(" -- ")[1] if " -- " in msg else msg
    print("ValueError:", head.split(". ")[0])
print()
print("An fd of 1e-6 is below what the shipped tables resolve, so it")
print("declines. Passing band, oversample and taps yourself always works.")

output

ValueError: no data: call set_data() before run()

An fd of 1e-6 is below what the shipped tables resolve, so it
declines. Passing band, oversample and taps yourself always works.

Validate the supplied calibration against your template population. A coarse pass is most useful when it dismisses many pairs; dense survivors can make it slower than flat filtering.

Devices and supported sizes

import matchedfilter as mf

print(mf.devices())
filt = mf.MatchedFilter(4096, ndata=16, ntemplates=64, device="gpu")

The default is CPU unless MF_DEVICE is set. device="cpu", "gpu", "gpu:1" and "auto" are supported. Linux and Windows use Vulkan; macOS uses Metal. GPU execution requires a compatible driver.

Capability	CPU	GPU
Peak-only transform sizes	powers of two, 64–1,048,576	powers of two, 64–65,536; device limits apply
Full-output transform sizes	powers of two, 1,024–4,194,304	powers of two, 1,024–4,194,304; device limits apply
Flat / hierarchical filtering	yes	yes
`run_series()`	yes	yes
Input spectra / output values	complex64	complex64
Flat and refinement arithmetic	float32	float32
Coarse screening	float32	may use reduced precision
Output bins	arbitrary count	split internally above 2048 bins
Execution	one CPU thread	selected device

GPU workgroup and shared-memory limits can reject individual sizes. Larger bin counts may repeat GPU transforms because output is processed in chunks. Unsupported requests raise an error; they do not switch devices.

CPU and GPU results can differ slightly because arithmetic ordering and coarse screening differ. Near-equal maxima can select different lags. Transform support is independent of hierarchical calibration coverage.

Shared arrays and input ownership

filter.empty_shared(shape) allocates NumPy arrays backed by GPU-accessible storage for GPU filters. Suitable contiguous complex64 banks can bind without an extra input copy. Call the setter again after changing a shared bank so cached coarse templates are refreshed.

For a full-correlation result that NumPy will read or scale, allocate out = filter.empty_shared((ndata, ntemplates, n), readback=True) and pass it to CorrelationFilter.run(out=out). On Vulkan, the default shared allocation favors GPU writes and can be slow for CPU reads. readback=True favors cached CPU access; ordinary NumPy out is another choice when a separate copy is acceptable. Metal uses shared storage for both settings. Automatic continuous run_series(series) already owns a readback-friendly result.

Host DLPack arrays are supported. Arbitrary CUDA or ROCm device allocations are not imported through this interface. Finish producer writes before filtering; cross-library GPU stream synchronization is not provided. Returned peaks use ordinary NumPy storage.

The CPU/GPU contract describes buffer ownership, cache lifetime and interoperability limits.

Platform checks

Linux x86-64, Linux arm64 and macOS arm64 are exercised in CI. That coverage does not establish performance or compatibility on every physical GPU. Metal was also tested on an Apple M2; historical measurements are in the GPU design notes. macOS x86-64 is not currently tested.

matchedfilter.targets() lists the available CPU targets, matchedfilter.backend() reports the selected target, and set_target() or MF_ISA selects a specific target for testing.

← Overview

See it work →