Metadata-Version: 2.4
Name: matchedfilter
Version: 0.1.0a1
Summary: Fast single-threaded batched matched filter with peak-only output (x86-64)
License-Expression: MIT
Project-URL: Homepage, https://github.com/ahnitz/matchedfilter
Keywords: matched-filter,correlation,fft,avx512,signal-processing
Classifier: Development Status :: 3 - Alpha
Classifier: Operating System :: POSIX :: Linux
Classifier: Environment :: Console
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: C
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.20
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Provides-Extra: bench
Requires-Dist: numpy>=1.20; extra == "bench"
Dynamic: license-file

# 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](#caveats).

```python
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

```bash
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:

| n | matchedfilter | numpy | |
|---:|---:|---:|---:|
| 1024 | 0.61 µs | 18.98 µs | 31x |
| 4096 | 2.08 µs | 37.81 µs | 18x |
| 16384 | 9.98 µs | 127.04 µs | 13x |
| 65536 | 48.87 µs | 568.56 µs | 12x |

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:

```bash
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.

```python
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](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

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

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

## License

MIT
