2 A Matched Filter on a single exposure
2.1 The special case of an isolated point source in given location
For an image of a single point source, the total number of counts recorded in each pixel of an image is:
Note that h can be a point spread function (scaled by the photometric solution), but could apply to more arbitrary mappings.
The final image will also contain noise:
When we are concerned with objects near the limits of detection, the photon noise from the sky will dominate; the noise will be uniform over the image.
If we wanted to measure the flux from the point source, we might simply estimate it from the brightest pixel:
We are not limited to using just that single pixel; we can get such an estimate for any pixel with flux:
but the uncertaitly will not be uniform:
We can do a weighted mean of these different estimates to get an estimate with minimal uncertainty:
Note that this is the value of the convolution of the image by its point spread function, evaluated at the position of the point source, and scaled by constant that depends on the point spread function and photometric solution:
Because the noise variance is independent of position, the uncertainty in the estimate of the flux is the the standard deviation of the sky noise in the image, after convolution by the point spread function. Because we know this particulary weighting minimizes the uncertainty in the flux relative to the flux, we can conclude that this weighting maximizes the value of the flux relative to the sky noise. This is precisely what we want for object detection; the point spread function is our optimal matching filter.
2.2 Example
First, I define a Moffat PSF:
moffat <- function(r,fwhm,beta=4.765) {
hwhm <- (fwhm/2.0)
alpha <- hwhm*(2^(1.0/beta)-1)^(-0.5)
p <- (1.0 + (r/alpha)^2.0)^(-1.0*beta)
p/sum(p)
}
Plot the raw signal, with the sky sigma equal to the square root of the sky and the flux in a Moffat profile with a total sky flux 5 times the sky sigma:
set.seed(53727)
npoints <- 64
sigma <- 4.0
flux <- 5*sigma
x <- (-1*npoints/2):(npoints/2)
h <- moffat(x,4.25)
n <- rnorm(x,sd=sigma)
b <- sigma^2+flux*h
g <- b+n
png("sample.png",height=256,width=640)
par(mar=c(1.1,1.1,0.1,8),bty="l",xaxt="n",yaxt="n",mgp=c(0,0,0))
plot(x,g,type="s",ylab="flux")
mtext("Data with noise",4,0,las=1)
lines(x,b,col="red",type="s")
mtext("Data w/o noise",4,0,las=1,padj=2,col="red")
dev.off()
Our representation of the pixel transfer function is centered on the center of the array rather than the edge, so the default R convolution shifts everything. Define a new convolution to center it on the 0 index.
recenter <- function(x) {
npts <- length(x)
hnpts <- ceiling(npts/2)
z <- c(x[(hnpts+1):npts],x[1:(hnpts)])
}
Now, apply the convolution
fhat <- recenter(convolve(g,h,type="circular"))
png("matched.png",height=256,width=640)
par(mar=c(1.1,1.1,0.1,8),bty="l",xaxt="n",yaxt="n",mgp=c(0,0,0))
plot(x,g,col="gray",type="s")
mtext("Data with noise",4,0,padj=-2,las=1,col="gray")
lines(x,fhat,type="s")
mtext("Matched filtered",4,0,padj=0,las=1)
lines(x,b,col="green",type="s")
mtext("Data w/o noise",4,0,las=1,padj=2,col="green")
dev.off()
2.3 The Matched filter in Fourier space
Frequency filtering alters both the signal and noise of an image, but because the power spectrum of the signal and noise differ, the filtering will alter the signal and noise differently. To detect an image, we need the signal from a given object to be as high as possible relative to the noise; in other words, we need to apply a filter that maximizes the signal to noise of the objects near the detection limit.
Near the detection limit, the noise will be approximately uniform, and the power spectrum of the noise will be flat; the Fourier transform of the noise will be of equal magnitude at all frequencies, and have a mean of zero. The Fourier transform of the signal, however, will be the Fourier transform of the object, multiplied by the Fourier transform of the PSF. In the case of a point source, the power spectrum will be the shape of the PSF, with a magnitude that varies with the flux of the object.
So, our goal is to derive a matched filter that maximizes the signal to noise. In other words, our filtered image:
should be a least squares fit, minimizing:
This is the Wiener filter, taking considering only low signal to noise objects as part of the signal. In Fourier space
If we are looking at low signal to noise objects (|B|2 < < |N|2 ) and we are looking at a random distribution of point sources (F(ν) is flat), then
and the matched filter is a mirror reflection of the point spread function, normalized by some factor (which doesn't affect the effectiveness of the matched filter).
2.4 A Fuzzy-minded explanation
What, then, does this mean? A pixel in the "true" image f(x) has contributions from each F(ν) Fourier components. Convolution by the point spread function by the instrument, atmosphere, etc. reduces this amount to H(ν)F(ν), and noise N(ν) is added to the result, giving us our observed image g(x), whose Fourier transform is:
Our best estimate of F(ν) will be
Our image g(x), however, does not have the same signal to noise in each frequency component. If we weight each frequency component by its (S/N)2 rather than doing a direct inverse transform, then we get
which is just the inverse transform of our match filtered image.
2.5 An example
Going back to our previous example,
GG <- abs(fft(g))^2
BB <- abs(fft(b))^2
NN <- abs(fft(n))^2
FFhat <- abs(fft(fhat))^2
k <- 1:(npoints/2)
png("matchedps.png",height=256,width=640)
par(mar=c(4,4,0.1,0.1),bty="l")
plot(k,GG[k],type="s",log="y",lwd=3,xlab="wavenumber",ylab="power")
lines(k,FFhat[k],type="s",col="blue")
lines(k,BB[k],type="s",col="green")
lines(k,NN[k],type="s",col="red")
leglabels<-c("Original data","Matched filtered data","Noise-free data", "Noise")
legcols<-c("black","blue","green","red")
legend(20,1e+06,leglabels,col=legcols,lty=1,bty="n")
dev.off()
Of course, when we are doing our fit, we do not know specifically what all the noise values are, so we use our expectation value for our magnitude of the noise:
NNe <- rep(0,times=npoints+1)+npoints*(sigma^2)
GGe <- rep(0,times=npoints+1)+BB+NNe
png("matchedpsue.png",height=256,width=640)
par(mar=c(4,4,0.1,0.1),bty="l")
plot(k,GG[k],type="s",log="y",lwd=2,col="gray",xlab="wavenumber",ylab="power")
lines(k,FFhat[k],type="s",col="blue")
lines(k,GGe[k],type="s",lwd=2)
lines(k,BB[k],type="s",col="green")
lines(k,NNe[k],type="s",col="red")
leglabels<-c("Original data","Matched filtered data","Expected signal+noise", "Noise-free signal", "Expected noise")
legcols<-c("gray","blue","black","green","red")
legend(20,1e+06,leglabels,col=legcols,lty=1,bty="n")
dev.off()
