matchedfilter

Batched correlation with a choice of full or peak-only output, on CPU or GPU.

matchedfilter correlates data segments against a template bank and returns the full complex correlation or the strongest sample in each output bin. It accepts complex64 spectra or time-series blocks through run_series().

Documentation · Usage guide · Benchmarks

Alpha software: the API may change. Pin a version for reproducible work.

Quick start

pip install matchedfilter

This example finds a template shifted by 37 samples:

import numpy as np
import matchedfilter as mf

n = 1024
rng = np.random.default_rng(1)
template = rng.standard_normal(n)
data = np.roll(template, 37)

filt = mf.MatchedFilter(n, ndata=1, ntemplates=1)
filt.set_templates(np.fft.fft(template).astype(np.complex64)[None, :])
filt.set_data(np.fft.fft(data).astype(np.complex64)[None, :])
peaks = filt.run(binsize=n)
print(peaks["index"][0, 0, 0])  # 37

The result has shape (ndata, ntemplates, nbins), with index and value fields. value is complex; use abs(value) for its magnitude. Bins below the detection threshold have index == -1 and value == 0. Results may reuse storage: copy any result you need to retain across calls.

See the usage guide for normalization, thresholds, search windows and time-series input.

Use CorrelationFilter when downstream code needs every lag:

full = mf.CorrelationFilter(n, ndata=1, ntemplates=1)
full.set_templates(np.fft.fft(template).astype(np.complex64)[None, :])
full.set_data(np.fft.fft(data).astype(np.complex64)[None, :])
correlation = full.run()  # complex64, shape (1, 1, n); peak at lag 37

It uses the same spectral inputs, bank dimensions, selectors and device choice as MatchedFilter. The full-output mode keeps the fused spectral product and batched transforms, but writes all lags. Peak-only filtering avoids those stores; hierarchical filtering also avoids full transforms for pairs dismissed by its coarse gate.

Supported capabilities

Flat filtering and refinement use single precision. GPU hierarchical screening can use reduced-precision kernels. There is no float64 filtering API.

The default device is CPU unless MF_DEVICE is set. Select a GPU explicitly:

print(mf.devices())
gpu_filter = mf.MatchedFilter(4096, device="gpu")

Linux and Windows use Vulkan; macOS uses Metal. A compatible driver is required. Unsupported GPU requests raise an error. See the usage guide for platform and device limits.

Hierarchical filtering

HierarchicalFilter screens each pair using a low-frequency band, then runs the full filter on candidates. It is useful when most pairs can be dismissed and enough signal power lies in that band. Otherwise the screening stage can add work without a useful saving.

Automatic selection uses your expected output-power spectrum and measured calibration files. snr specifies the signal strength and fd the requested false-dismissal budget for that calibration. Validate the calibration against your template population before relying on it.

For example, with spectra and an output-power reference prepared as described in the usage guide:

hf = mf.HierarchicalFilter(4096, ndata=16, ntemplates=64, snr=5.5, fd=1e-2)
hf.set_reference(expected_output_power)
hf.set_templates(template_spectra)
hf.set_data(data_spectra)
peaks = hf.run(binsize=4096, threshold=5.5)

The coarse gate is computed from the full reference profile and requested false-dismissal budget; measured cost files select the band. Alternatively, set an explicit band and coarse threshold. See the gate model for profile assumptions, sampling precision and validation.

Performance

The implementation fuses the frequency-domain product into the inverse transform. Full-output mode writes every lag; peak-only mode scans inside the transform and avoids those writes. Hierarchical mode screens pairs before the full transform. Performance depends on length, bank shape, device and the fraction of pairs requiring refinement.

Measured on a Ryzen AI MAX+ 395 / Radeon 8060S on 2026-10-13 at revision 2b3cac13: 512 data segments × 512 templates (262,144 pairs), 4,096 points, with the bank drawn from the captured PyCBC reference profile. The bars move from a general inverse FFT to fused full output, peak-only output and hierarchical screening. FFTW and rocFFT time only the inverse transform, so they do less work than the filter bars. Full output reuses a caller-supplied result array; GPU timings include synchronization. Hierarchical bars show requested FDR budgets 1e-2, 1e-3 and 1e-4 at SNR threshold 5.5, with initial autotuning passes settled to empirical best configurations; these noise-only timings do not measure FDR. The two panels use separate scales; compare their printed values. This is a workload example, not a speed guarantee.

For CPU and GPU comparisons across eight machines at 16 × 512 and 512 × 512 (including NVIDIA L40S, A100 and A40, AMD Radeon 8060S and Apple M2), with selectable thresholds and output modes, see the hardware comparison.

The flat and hierarchical benchmark pages show results across available transform sizes.

Run the benchmarks

python -m pip install pyfftw  # optional FFTW reference
python -m matchedfilter.benchmark --reps 7 --json bench.json

The default sweep covers all 15 CPU sizes, with GPU timings where supported. Batch sizes shrink at large lengths to bound memory use. --n selects a subset. Missing hierarchical calibration is reported explicitly.

To measure warm overlap-save calls, including automatic peak output at 2,048–8,192 points and continuous full output at 32,768 points, run python tools/bench_series_workloads.py. It checks the output before timing and compares automatic layout with explicit blocks. Use --quick for a shorter continuous series; CI records that version on each benchmark host.

Install mkl-fft and mkl for an optional MKL reference where supported. NumPy provides the correctness reference. Without FFTW or MKL, library timings and correctness checks still run. Reference timing covers the inverse FFT; matchedfilter timing also includes the product and peak scan.

Documentation

• Documentation index: current guides, development tools, and research history.
• Usage guide: inputs, normalization, devices and API behavior.
• Examples: noise and injected signals.
• Numerical accuracy: comparison with a float64 reference.
• Design notes: implementation details and historical experiments.

Status

This is an alpha release. Pin a version when reproducibility matters. Hierarchical calibration covers a subset of input conditions; unsupported requests require additional measurements or explicit coarse parameters.

Wheels target CPython 3.10–3.14 on Linux x86-64 and macOS arm64. Other platforms build from source and need NumPy and a C compiler. CPU filtering uses one thread; callers control parallelism across independent filters.

Contributing

Issues and pull requests are welcome. Useful contributions include device-specific benchmarks, calibration measurements and testing on additional hardware.

git clone https://github.com/ahnitz/matchedfilter
cd matchedfilter
git submodule update --init third_party/highway
pip install -e .
pytest

GPU kernels are written in Slang. Rebuild the shipped SPIR-V and Metal sources with python tools/build_spirv.py --slangc /path/to/slangc. The tooling map identifies benchmark, calibration and historical research scripts. The performance checks include opt-in timing gates and comparison with a saved baseline.

Regression Benchmark Suite & Invariant Verification

To safeguard against throughput drops, threshold contract violations, and autotuning regressions across key interfaces:

# Run via pytest (asserts invariants, noise contracts, and injection recovery fidelity in < 2s):
pytest tests/test_regression_suite.py

# Or run directly for detailed latency, throughput, and status reporting:
python tests/test_regression_suite.py

Development Procedure: Whenever an issue (fidelity bug, threshold contract violation, autotuning failure, or performance regression) is identified and resolved, add a dedicated regression test case and metric assertion to tests/test_regression_suite.py to permanently lock in the improvement.

License

MIT