Matched filters and detection

A known pulse, a gliding tone, is hidden somewhere in 512 samples of noise about as strong as the pulse. Slide a copy of the pulse along the record, and at each position multiply and add. Watch the output: where the copy lines up with the hidden pulse, it jumps.

Slide the template, find the pulse

You can pick out a friend’s voice in a noisy crowd, because you know exactly how it sounds. A radar receiver has the same advantage. It sent the pulse itself, so it knows the shape of the echo it listens for. What it does not know is when the echo comes back.

Here is the problem in symbols. A known pulse $q[n]$, the template, sits somewhere in a record, scaled by an amplitude $A$ and buried in noise:

The template is $N_q$ samples long and starts at $n_0$, which I do not know. The noise $v[n]$ is white, as in “White noise forgets, coloured noise remembers” of Random processes (24.2), with standard deviation $\sigma_v$.

My template is a chirp, the gliding tone of “Reading a chirp” in Spectrograms & the STFT (15.5). It has $N_q=64$ samples, $q[k]=\cos\big(\pi(0.05k+0.20k^2/63)\big)$ for $k$ = 0 to 63. Its frequency, the slope of its phase, glides from $0.05\pi$ to $0.45\pi$ rad/sample.

Its energy, the sum of its squares from “Adding it all up” in How big is a signal (1.3), is $E_q=32.55$. That is close to half of 64, since the squares of a cosine average about one half.

In the first instrument the record has 512 samples. The chirp has amplitude $A=1$ and starts at $n_0=300$, and the noise has $\sigma_v=1$. It is a seeded draw (seed 261), so it is the same on every visit; this draw measures a standard deviation of 0.985.

To find $n_0$, I do what Correlation (24.3) did in “Slide, multiply, add: no flip”. Slide the template along the record, multiply sample by sample, and add. At lag $\ell$ the template’s first sample sits on $x[\ell]$:

This is 24.3’s $R_{xq}[\ell]=\sum_nx[n]\,q[n-\ell]$, with $n=\ell+k$: the record stays still and the template slides. So the peak should land at $\ell=n_0$, where the record lags the template by $n_0$. The lags run from 0 to 448, the last place where all 64 samples of the template fit inside the record.

From statistics this page assumes only this: the noise is white and Gaussian, its power is known, and so is the pulse’s shape. Only the pulse’s position is unknown.

The picture at the top of the page slides the template along this record, one lag at a time. Notice the end frame: the peak, 31.8 at lag 300, stands 5.8 noise standard deviations high, and it is narrow. One lag to either side, at 299 and 301, the output is already down to 17.6 and 18.8.

Look back at the record panel, too. The chirp’s amplitude, 1, is about the noise’s standard deviation, so nothing there shows where it starts. The correlation finds it anyway.

Why the peak stands so tall

At the right lag, $\ell=n_0$, each product of the pulse with the template is $A\,q[k]^2$, never negative. So they all add up, to $A$ times the energy:

The noise’s products have random signs, so they partly cancel. Their sum $\sum_kv[\ell+k]\,q[k]$ adds 64 separate draws, each scaled by $q[k]$. Scaling a draw by $q[k]$ multiplies its variance by $q[k]^2$, since the variance is an average of squares.

By “Sums of draws” in Random variables for signals (24.1), the variances of separate draws add. So the noise in the output has variance

and standard deviation $\sigma_v\sqrt{E_q}$, at every lag.

The peak stands $A\,E_q/(\sigma_v\sqrt{E_q})=A\sqrt{E_q}/\sigma_v$ noise standard deviations high. Square that, and you have a ratio of powers, the output SNR:

Only the pulse’s energy counts, not its shape. For the chirp, $\sqrt{E_q}=5.71$, so with $A=1$ and $\sigma_v=1$ the peak should stand 5.71 standard deviations high.

In this record the peak is 31.8. The chirp alone would give 32.55, and the noise at that lag adds −0.77. Where the template no longer overlaps the pulse, 64 lags or more away, the output’s standard deviation measures 5.52, against 5.71 in theory. So the peak stands $31.8/5.52=5.8$ of them high.

With this draw’s measured 0.985 in place of $\sigma_v=1$, the theory says 5.62. The 322 lags are not separate draws, since neighbouring lags share 63 of their 64 samples, so 5.52 is a rough estimate.

In decibels, as in “Signal and noise, in decibels” of 1.3: during the pulse, the chirp’s power per sample is $E_q/N_q=0.509$ against the noise’s 1. That is an SNR of −2.9 dB. With this draw’s measured noise RMS, 0.987, it is −2.8 dB.

After the filter the SNR is $10\log_{10}32.55=15.1$ dB. The filter has multiplied the SNR by $N_q=64$, which is 18 dB. That holds for any pulse of $N_q$ samples, because $E_q$ is $N_q$ times the pulse’s power per sample.

No filter does better

Could other weights beat the template? Say a filter weights the 64 samples by some $g[k]$ instead of $q[k]$. The same two steps give the pulse part $A\sum_kg[k]\,q[k]$, and the noise’s standard deviation $\sigma_v\sqrt{\sum_kg[k]^2}$. Their ratio is

The last fraction is the normalised correlation of “Between −1 and 1” in 24.3, of $g$ with $q$ at lag 0. It is at most 1, and it reaches 1 only when $g$ is a positive multiple of $q$. Weights outside the pulse’s 64 samples collect only noise, so they cannot help either.

So in white noise no linear filter lifts the pulse higher above the noise than the template does: $A\sqrt{E_q}/\sigma_v$ is the most it can reach. This is the Cauchy–Schwarz inequality, which linear algebra proves.

The template, played backward

So far I have correlated. To do it with a filter, recall “A convolution with one signal reversed” in 24.3: convolving with a reversed signal is correlating with it. The reversal is “Playing it backward: reversal” of Shifting, reversing and scaling time (2.1). The flip in “Flip and slide: a mechanical way to compute it” of Discrete convolution (5.2) undoes it.

So the filter’s impulse response is the template played backward, moved to start at $n=0$:

For the chirp, $h[n]=q[63-n]$: it starts with the template’s last sample. The move by $N_q-1$ keeps $h$ at $n$ = 0 to 63, so the filter is causal and needs no future samples.

Its output, with $k=N_q-1-m$, is the correlation again:

The output is the correlation, 63 samples late. Its peak comes at $n=n_0+63=363$, the pulse’s last sample. The filter can add up the whole pulse only once the whole pulse has arrived.

This filter, matched to one pulse, is the matched filter. Built as a sliding correlation instead, the same receiver is often called a correlation receiver.

Everything here assumed white noise. If the noise is coloured, first pass the record through a filter that makes the noise white, such as the prediction-error filter of Parametric models and linear prediction (25.3). Then match the template as it looks after that filter.

Where to set the threshold

The peak tells me where the pulse is. A receiver must also decide whether there is a pulse at all. A weak echo may stand only a little above the noise, and noise alone sometimes rises high.

Think of a smoke alarm. Set it sensitive, and it rings for toast. Set it dull, and it stays silent through a small fire. A detector has the same dial.

First I make the output’s scale simple. Divide it by its noise standard deviation, $\sigma_v\sqrt{E_q}$. Without a pulse, the result at any lag has mean 0 and standard deviation 1. It is also Gaussian, because a weighted sum of Gaussian draws is Gaussian; I assume this; statistics proves it.

With the pulse at that lag, the same bell moves to

The noise is the same, so the bell keeps its width; only its centre moves. On this page $d$ is this distance between the two bells, in noise standard deviations.

For this section the chirp is weaker, with amplitude $A=0.35$. Then $d=0.35\times5.7055=1.997$, which the instrument shows as 2.00.

Now pick a threshold $\eta$, in these units, and declare “pulse” whenever the output passes it.

Two things can go wrong. Noise alone can pass $\eta$: that is a false alarm. A pulse can stay under it: that is a miss. The share of pulses that do pass is the detection rate.

The false-alarm rate is the area of the no-pulse bell beyond $\eta$. For whole numbers, 24.1’s areas give it. Within 1, 2 and 3 standard deviations lie 68.27 %, 95.45 % and 99.73 % of the draws. The bell is symmetric, so half of the rest lies beyond $+\eta$: 15.87 %, 2.28 % and 0.13 % for $\eta$ = 1, 2 and 3.

For any $\eta$, this tail area has a name. The complementary error function $\operatorname{erfc}$ is a tabulated function that gives Gaussian tail areas; SciPy has it as scipy.special.erfc. The false-alarm rate is

The detection rate is the same area for the moved bell, $\tfrac12\operatorname{erfc}\big((\eta-d)/\sqrt2\big)$.

Each threshold gives a pair of rates. Plot the detection rate against the false-alarm rate, one point per threshold, and the points trace the ROC curve, short for receiver operating characteristic.

Watch the dot trace the curve as the threshold slides from 3 to 1: the detections rise, and so do the false alarms.

After the clip, move the threshold yourself with the slider.

Notice the middle frame: at $\eta=2$ about half the pulses are caught, 49.9 %, at 2.28 % false alarms.

Reading the curve

Check that frame by hand with $d=2$. The bell of the pulse is centred on $\eta=2$, so half of it lies beyond: 50.0 %. The readout says 49.9 % because the instrument uses $d=1.997$, slightly below 2. In the same way, $d=2$ gives 15.9 % at $\eta=3$, against the readout’s 15.8 %; at $\eta=1$ both give 84.1 %.

Lowering $\eta$ moves the dot along the curve towards the corner where both rates are 1: more pulses caught, and more false alarms. Raising it moves the dot back towards the corner where both rates are 0. One rate improves only when the other gets worse.

Now the diagonal. Imagine a detector that ignores the record and declares “pulse” at random, a share $p$ of the time. It catches a share $p$ of the pulses and raises a share $p$ of false alarms, so it sits on the diagonal. Every point of the matched filter’s curve lies above it.

A stronger pulse moves the second bell further out. Then the curve bends further towards the corner with false-alarm rate 0 and detection rate 1, where every pulse is caught and noise never is.

Measured, not assumed

The bells are a model. To test it, I ran 2000 trials without the pulse and 2000 with it, each with 64 fresh draws through the matched filter (seed 2610). The outputs averaged 0.005 and 1.967, against the model’s 0 and 1.997. Their standard deviations were 0.974 and 1.010, against the model’s 1.

At $\eta$ = 3, 2 and 1 the false alarms were 0.10 %, 1.95 % and 15.05 %. The detections were 15.5 %, 48.8 % and 83.6 %.

How close should they be? Each trial either crosses $\eta$ or not: a draw that is 1 with probability $p$ and 0 otherwise. Its mean is $p$, and its variance is $p-p^2$, because its square is itself.

The share over 2000 separate trials adds 2000 such draws and divides by 2000. The variances add, and the division scales the variance by $1/2000^2$, as before. So the share has spread $\sqrt{p(1-p)/2000}$.

At $\eta=1$ that is 0.82 %, and the measured 15.05 % lies 0.82 below the bell’s 15.87 %, just under one spread. Each of the six measured rates lies less than one spread from the bells’ value.

These rates are for one decision. A search like the first instrument’s makes one decision at every lag, so the noise gets many chances. In that record, of the 322 lags at least 64 away from the pulse, one passed $\eta=3$: lag 384, where the output was 17.14, just above $3\sqrt{E_q}=17.12$.

The chirp of the first instrument, with $A=1$, has $d=5.71$, and $\eta=3$ catches 99.7 % of such pulses. So a search over many positions sets its threshold high, and pays for it with a strong pulse.

In practice you first decide how many false alarms you can afford, then set $\eta$ to give that rate. This is the Neyman–Pearson rule.

Worked example

1. A matched filter by hand. Take the template $q$ = 1, 2, −1 at $n$ = 0 to 2, so $N_q=3$ and $E_q=1+4+1=6$. Played backward, $h[n]=q[2-n]$ is −1, 2, 1. Let the record be the template starting at $n_0=2$, with no noise: $x$ = 0, 0, 1, 2, −1 at $n$ = 0 to 4.

Convolving $x$ with $h$ gives 0, 0, −1, 0, 6, 0, −1 at $n$ = 0 to 6. At $n=4$, for example, the products are $(-1)(-1)+2\cdot2+1\cdot1=6$. The peak, 6, is $E_q$, at $n=n_0+N_q-1=4$, the pulse’s last sample. The correlation itself peaks at lag 2, which is $n_0$.

2. Output SNR. For the chirp, $E_q=32.55$, with $A=1$ and $\sigma_v=1$, gives $\mathrm{SNR}_\text{out}=10\log_{10}32.55=15.1$ dB, against −2.9 dB per sample during the pulse. The peak should stand $\sqrt{E_q}=5.71$ standard deviations high; this record’s stands 5.8.

3. One threshold. At $\eta=2$, with $d=2.00$: false alarms $\tfrac12\operatorname{erfc}(2/\sqrt2)=2.28$ %. Detections are half the pulse bell, 50.0 % by hand with $d=2$, and 49.9 % with the full $d=1.997$.

4. How strong must the pulse be? Keep the false alarms at 0.13 %, so $\eta=3$. To catch 84.1 % of pulses, the pulse bell must sit one standard deviation beyond $\eta$, at $d=4$, as 24.1’s 68.27 % shows. That needs $A=4/\sqrt{E_q}=4/5.7055=0.70$, twice the 0.35 of the instrument.

Where you’ll meet this

Radar and sonar correlate each echo with the pulse they sent. They often send a chirp, as here: it lasts long, so it carries much energy, yet its peak is narrow. The chirp of this page lasts 64 samples, but its peak, alone, is only 3 lags wide at half its height. Squeezing a long pulse into a short peak is called pulse compression, the subject of Radar signal processing (33.1).

A GPS receiver correlates what it hears with each satellite’s known code, at many delays, and declares a satellite found when the peak passes a threshold. Image software slides a small template over a picture to find where it appears. Neuroscientists find nerve-cell spikes in noisy electrode recordings by matching a spike’s known shape.

The matched filter finds a pulse whose shape you know. When the signal itself is unknown and you want its whole waveform back from the noise, you need The Wiener filter (26.2).

For more, see A. V. Oppenheim and G. C. Verghese, Signals, Systems and Inference (2015), the chapters on hypothesis testing and signal detection, and H. L. Van Trees, Detection, Estimation, and Modulation Theory, part I. In SciPy, scipy.signal.correlate(x, q, mode='valid') gives the output at every lag, and scipy.special.erfc the tail areas.

Reference card