matchedfilter

A fast single-threaded matched filter for x86. You give it a batch of data segments and a batch of templates; it correlates every pair and hands back only the peaks.

Not returning the full correlation is the point. Most searches threshold the output and throw the rest away, and once you say so up front the filter can skip work that could not have produced a peak anyway.

Status: work in progress. The API still moves, and there is a known calibration weakness in the hierarchical filter. See Caveats.

import matchedfilter as mf

filt = mf.MatchedFilter(16384, ndata=16, ntemplates=64)
filt.set_data(data_spectra)          # (16, 16384) complex64, already FFT'd
filt.set_templates(template_spectra) # (64, 16384) complex64

peaks = filt.run(binsize=1024, threshold=5.5)
peaks["index"], peaks["value"], peaks["magnitude"]

import matchedfilter as mf is the convention used throughout these docs.

Install

pip install --pre matchedfilter

--pre because this is an alpha release. Wheels are built for CPython 3.9 to 3.13 on manylinux x86-64; anywhere else pip falls back to the source distribution, which needs numpy and a C compiler.

x86-64 only for now. The kernels are AVX2 and AVX-512 intrinsics with no portable fallback, so a build on any other architecture stops with an error rather than producing a slow one. AVX-512 is used when the CPU has it and AVX2 otherwise, decided at runtime.

How it works

Inputs are frequency domain: the unnormalised forward transform of each segment, in natural order. Produce them with whatever you already use (numpy, MKL, FFTW); matchedfilter does not need to own that step.

The filter is built once and reused. Ingest conjugates the templates and stores both sides in the layout the correlation loop walks, which costs a few percent of a run and less as the batch grows.

run returns a structured array of shape (ndata, ntemplates, nbins) with fields index, value and magnitude. Bins whose peak fell below the threshold carry index == -1.

Supported lengths are 1024 and the powers of two from 4096 to 1048576.

Performance

Per (data, template) pair, 8x32 batch, one core of a Zen 5 desktop:

numpy is a floor, not a rival. It is there so the comparison runs anywhere. Against MKL or FFTW the margin is much smaller, and part of what is left comes from computing peaks instead of a full correlation. Measure on your own box:

python -m matchedfilter.benchmark

Hierarchical filtering

HierarchicalFilter adds a cheap pre-pass: correlate against a low-frequency slice of the template, and only run the full-length filter where that slice leaves a peak plausible.

This helps only under an assumption about your templates: that enough of the matched-filter output power sits in the low band that a narrow slice gives a usable bound on the full result. For chirp-like templates whose power is concentrated at low frequency that tends to hold. For templates whose power is spread flat across the band, or concentrated high, the slice bounds nothing useful, and the pre-pass is pure added cost. It is worth checking against your own templates before relying on it.

hf = mf.HierarchicalFilter(16384, ndata=16, ntemplates=64,
                           snr=6.0,   # threshold you intend to use
                           fd=1e-3,   # false-dismissal budget
                           band=2048) # width of the cheap slice

hf.set_reference(expected_output_power)      # real frequency series, the OUTPUT
hf.set_templates(template_spectra)
hf.set_data(data_spectra)
peaks = hf.run(binsize=16384, threshold=6.0)

set_reference takes a real frequency series of length n holding the expected power of the filter output in each bin. Only its shape is used; the overall normalisation is divided out.

This is the output, not the template. The two differ whenever the data is coloured, and passing the template's own power will mis-set the gate: a broadband template reconstructing a narrowband signal is the case where it goes wrong by the largest factor.

The gate is one-sided by construction: peaks it reports are bit-identical to the flat filter's. It can only omit, never invent. fd is the budget for how often it is allowed to omit one.

On pure noise at n=4096 the gate runs about 7x faster than the flat filter. The saving scales with how little survives, so it grows with your threshold and falls toward 1x on data where most pairs trigger.

Caveats

• The false-dismissal budget is not currently met at low thresholds. fd is honoured well at snr 6 and above. At snr 5.0 to 5.5 with a coarse band the gate omits more than it should: 1.4% against a 0.1% budget in a 418-template search. The cause is that the gate's recovery factors are measured from a mean frequency series, which is not a bound on any individual realisation. Tracked by an xfail test in tests/test_api.py and written up in docs/hierarchical.md.
• Single-threaded by design. Parallelism is the caller's to arrange.
• x86-64 Linux only, as above. Other architectures are not implemented rather than merely untested.

Development

pip install -e .[test]
pytest              # includes a C test that builds itself from source
python -m matchedfilter.benchmark

docs/ holds the design notes: how the hierarchical gate is calibrated, and measurements of the approaches that were tried and rejected.

License

MIT