# Matched Filter (Convolution)¶

For source detection, SEP supports using a matched filter, which can
give the optimal detection signal-to-noise for objects with some known
shape. This is controlled using the `filter_kernel`

keyword in
`sep.extract`

. For example:

```
kernel = np.array([[1., 2., 3., 2., 1.],
[2., 3., 5., 3., 2.],
[3., 5., 8., 5., 3.],
[2., 3., 5., 3., 2.],
[1., 2., 3., 2., 1.]])
objects = sep.extract(data, thresh, filter_kernel=kernel)
```

If `filter_kernel`

is not specified, a default 3-by-3 kernel
is used. To disable filtering entirely, specify `filter_kernel=None`

.

What array should be used for `filter_kernel`

? It should be
approximately the shape of the objects you are trying to detect. For
example, to optimize for the detection of point sources,
`filter_kernel`

should be set to shape of the point spread function
(PSF) in the data. For galaxy detection, a larger kernel could be
used. In practice, anything that is roughly the shape of the desired
object works well since the main goal is to negate the effects of
background noise, and a reasonable estimate is good enough.

## Correct treatment in the presence of variable noise¶

In Source Extractor, the matched filter is implemented assuming there
is equal noise across all pixels in the kernel. The matched filter
then simplifies to a convolution of the data with the kernel. In
`sep.extract`

, this is also the behavior when there is constant noise
(when `err`

is not specified).

In the presence of independent noise on each pixel, SEP uses a full matched filter implementation that correctly accounts for the noise in each pixel. This is not available in Source Extractor. Some benefits of this method are that detector sensitivity can be taken into account and edge effects are handled gracefully. For example, suppose we have an image with noise that is higher in one region than another. This can often occur when coadding images:

```
# create a small image with higher noise in the upper left
n = 16
X, Y = np.meshgrid(np.arange(n), np.arange(n))
mask = Y > X
error = np.ones((n, n))
error[mask] = 4.0
data = error * np.random.normal(size=(n, n))
# add source to middle of data
source = 3.0 * np.array([[1., 2., 1.],
[2., 4., 2.],
[1., 2., 1.]])
m = n // 2 - 1
data[m:m+3, m:m+3] += source
plt.imshow(data, interpolation='nearest', origin='lower', cmap='bone')
```

Specifying `filter_type='conv'`

will use simple convolution, matching the
behavior of Source Extractor. The object is not detected:

```
>>> objects = sep.extract(data, 3.0, err=error, filter_type='conv')
>>> len(objects)
0
```

Setting `filter_type='matched'`

(the default)
correctly deweights the noisier pixels around the source and detects
the object:

```
>>> objects = sep.extract(data, 3.0, err=error, filter_type='matched')
>>> len(objects)
1
```

## Derivation of the matched filter formula¶

Assume that we have an image containing a single point source. This produces a signal with PSF \(S_i\) and noise \(N_i\) at each pixel indexed by \(i\). Then the measured image data \(D_i\) (i.e. our pixel values) is given by:

Then we want to apply a linear transformation \(T_i\) which gives an output \(Y\):

We use matrix notation from here on and drop the explicit sums. Our objective is to find the transformation \(T_i\) which maximizes the signal-to-noise ratio \(SNR\).

We can expand the denominator as:

Where \(C_{ik}\) is the covariance of the noise between pixels \(i\) and \(k\). Now using the Cauchy-Schwarz inequality on the numerator:

since \(C^T = C\). The signal-to-noise ratio is therefore bounded by:

Choosing \(T = \alpha C^{-1} S\) where \(\alpha\) is an arbitrary normalization constant, we get equality. Hence this choise of \(T\) is the optimal linear tranformation. We normalize this linear transformation so that if there is no signal and only noise, we get an expected signal-to-noise ratio of 1. With this definition, the output \(SNR\) represents the number of standard deviations above the background. This gives:

Putting everything together, our normalized linear transformation is:

And the optimal signal-to-noise is given in terms of the known variables as: