Matched Filters

The inner product of a normalized matched filter vector (h) with a normalized measurement vector (m) is a number from 0 (no target) to 1 (target clearly present). Matched filters are built from a target signal (s) and a covariance matrix of the anticipated noise (R).

One application of matched filters is in spectrometry, where they are used to detect spectral signatures. Matched filters are particularly effective here because they use the correlations between spectral channels (common in sensors) to optimize SNR. Figure 1 shows fake spectral measurements of varying noise levels along with their inner products. The inner products decrease as the noise increases. Detection becomes a matter of choosing the right threshold value, for example, balancing false positives against false negatives using a receiver operating characteristic (ROC) curve.

_images/mf_spectralSignature.png

Figure 1. Spectra with varying noise levels and their matched-filter inner products

The general matched filter formula is derived here. We begin by expressing SNR in terms of the inner product (y) of the matched filter and the measurement, which is sum of the signal and noise, \(m = s+n\):

\[\begin{split}y &= \sum_{k}{h^*[k] m[k]} = h^Hm = h^Hs + h^Hn \\ SNR &= \frac{|h^Hs|^2}{E[|h^Hn|^2]}\end{split}\]

The objective is to find the h that maximizes the SNR expression. If the noise is zero-biased, the denominator can be rewritten as

\[E[(h^Hn)(h^Hn)^H] = h^HE[nn^H]h = h^HRh \qquad\textrm{(0)}\]

where R is the noise covariance matrix. The SNR expression can be simplified by first replacing h with \(R^{1/2}h'\), where the apostrophe indicates a transpose operation:

\[SNR = \frac{|h^Hs|^2}{h^HRh} = \frac{|(R^{1/2}h')^Hs|^2}{(R^{1/2}h')^HR(R^{1/2}h')} = \frac{|(R^{1/2}h')^H(R^{-1/2}s)|^2}{(R^{1/2}h')^H(R^{1/2}h')} \qquad\textrm{(1)}\]

The numerator has the same form as the left side of the Cauchy-Schwarz inequality, \(|a^Hb|^2 \leq (a^Ha)(b^Hb)\), which essentially means that the projection of one vector onto another can’t be greater than the longer of the two vectors. Replacing the numerator of (1) with the Cauchy-Schwarz inequality (\(a=R^{1/2}h'\) and \(b=R^{-1/2}s\)) leads to a tidy upper bound for SNR:

\[\begin{split}SNR \leq \frac{[(R^{1/2}h')^H(R^{1/2}h')][(R^{-1/2}s)^H(R^{-1/2}s)]}{(R^{1/2}h')^H(R^{1/2}h')} &= \frac{(R^{-1/2}s)^H(R^{-1/2}s)}{1} \\ &= s^HR^{-1}s.\end{split}\]

This latest result is equal to the rightmost term in (1) when \(R^{1/2}h'=\alpha\ R^{-1/2}s\):

\[SNR = \frac{|(R^{1/2}h')^H(R^{-1/2}s)|^2}{(R^{1/2}h')^H(R^{1/2}h')} = \frac{|(\alpha\ R^{-1/2}s)^H(R^{-1/2}s)|^2}{(\alpha\ R^{-1/2}s)^H(\alpha\ R^{-1/2}s)} = s^HR^{-1}s.\]

In other words, the SNR is maximized when the \(R^{1/2}h' = \alpha\ R^{-1/2}s\) condition is met. Or, equivalently, when

\[h = \alpha\ R^{-1}s \qquad\textrm{(2)}\]

We know from (0) that the noise is \(h^HRh\). Normalizing this, we can solve for alpha:

\[\begin{split}(\alpha\ R^{-1}s)^HR(\alpha\ R^{-1}s) = 1 \\ \alpha^2s^HR^{-1}s = 1 \\ \alpha = \frac{1}{\sqrt{s^HR^{-1}s}} \qquad\textrm{(3)}\end{split}\]

Combining (2) and (3), we get our general matched filter equation:

\[h=\alpha\ R^{-1}s = \frac{R^{-1}s}{\sqrt{s^HR^{-1}s}}.\]

So, to create a matched filter, we need to know the target that we’d like to find (s), and we need to have an inversion of the estimate of the covariance of the background noise.

 1# Python 3.x
 2import numpy as np
 3import matplotlib.pyplot as plt
 4from scipy.linalg import cholesky, inv, eigh
 5from scipy.stats import norm
 6
 7def fakeSpectralSignature(N):
 8    x = np.arange(N)
 9    s = 50*norm.pdf(x,70,10)
10    s += 300*norm.pdf(x,250,30)
11    s += 75*norm.pdf(x,480,10)
12    return s
13
14def covarianceMatrix(N):
15    R = np.zeros([N,N])
16    for r in range(N):
17        for c in range(r,N):
18            x = 1 if r == c else 0.01*np.random.random(1)
19            R[r,c] = x
20            R[c,r] = x
21#    d,v = eigh(R)
22#    c = np.dot(v, np.diag(np.sqrt(d))).transpose() # CC^T=R
23    c = cholesky(R, lower=True) # ~40x faster   
24#    print(c.transpose()@c, '\n', R)
25    return R, c
26
27def correlatedNoise(scalingFactor, c):
28   n = scalingFactor*np.dot(c, norm.rvs(size=c.shape[0]))
29   return n - np.mean(n) # noise should be zero-biased
30
31def project(s, Rinv, c, nFactor, ax):
32    n = correlatedNoise(nFactor, c)
33#    print(np.cov(n))
34    m = n + s
35    h = (Rinv @ s)/(s.transpose() @ Rinv @ s)**0.5
36    h /= np.sum(h**2)**0.5 # make into a unit vector
37    m /= np.sum(m**2)**0.5
38    projection = np.dot(h.transpose(), m)
39    ax.plot(m, 'k')
40    ax.set_xlabel('Spectral Channel')
41    ax.set_title('$\hat{h}^T•\hat{m}$ = %0.3f' % projection)
42    ax.get_yaxis().set_visible(False)
43    
44Nch = 512
45s = fakeSpectralSignature(Nch)
46R,c = covarianceMatrix(Nch)
47Rinv = inv(R)
48
49fig = plt.figure(1); fig.clf()
50project(s, Rinv, c, 0, fig.add_subplot(321))
51project(s, Rinv, c, 0.1, fig.add_subplot(322))
52project(s, Rinv, c, 0.5, fig.add_subplot(323))
53project(s, Rinv, c, 1, fig.add_subplot(324))
54project(s, Rinv, c, 3, fig.add_subplot(325))
55project(s, Rinv, c, 50, fig.add_subplot(326))
56plt.tight_layout()