text:
# Comparing filters for the detection of point sources

## 1 Introduction

The detection of localized signals or features on one-dimensional (1D) or two-dimensional images (2D) is one of the most challenging aspects of image analysis.

We are interested in the detection of compact sources (signal) embedded in a background (1D case). We assume that the profile of the source and the statistical properties of the background are known. Linear filtering of the data in order to eliminate partially the background is the primary goal. Several filters have been introduced in the literature to deal with the problem. Four examples are: the continuous ‘Mexican hat’ wavelet, difference of two Gaussians, matched filters and scale-adaptive filters. The two first cases are filters given ‘a priori’, adapted to the detection of point sources, whereas the matched filter is constructed taking into account the profile and background in order to get the maximum SNR at the source position. The scale-adaptive filter is constructed taking into account the previous properties and also the constrain to have a maximum in filtered space at the scale and source position. Hereinafter, we will identify the position of the possible sources with the local maxima.

The ‘Mexican hat’ wavelet used as a filter has been extensively applied during the last years to analyse optical, X-ray and microwave data. Optical images of galaxy fields have been analysed to detect voids and high-density structures in the first CfA redshift survey slice (Slezak et al. 1993). Microwave images have been analysed (Cayón et al. 2000; Vielva et al. 2001a) and combined with the maximum entropy method (Vielva et al. 2001b) to obtain catalogues of point sources from simulated maps at different frequencies that will be observed by the future Planck mission. On the other hand, the Mexican hat has also been used to detect X-ray sources (Damiani et al. 1997) and presently for the on-going XMM-Newton mission (Valtchanov et al. 2001).

Other useful filters include the so-called ‘matched’ filters. They optimize the signal-to-noise ratio and have been used mainly in signal processing. The generalization to two dimensions of the previous matched filters have been used recently to detect clusters of galaxies from optical imaging data (Postman et al. 1996; Kawasaki et al. 1998). In this approach the method uses galaxy positions, magnitudes (and photometric/spectroscopic redshifts if available) to find clusters and determine their redshift. Also they have been applied to microwave maps to detect clusters through the Sunyaev-Zeldovich effect, either on single or multifrequency maps (Herranz et al. 2002a, b)

Other remarkable filters that have been recently (Sanz et al. 2001) introduced in the astrophysics literature are the so-called scale-adaptive filters. Applications have been done to the detection of point sources in time-ordered data in the future Planck mission (Herranz et al. 2002c), detection of SZ-emission in single microwave maps and X-ray emission from clusters of galaxies in single X-ray maps (Herranz et al. 2002a). More recently, application for the detection of SZ-clusters in multifrequency maps representing the future observations by the Planck mission has been considered (Herranz et al. 2002b).

One interesting question relative to all of these filters is the $`\mathrm{𝑜𝑝𝑡𝑖𝑚𝑎𝑙𝑖𝑡𝑦}`$, that we define in terms of the following properties: confidence level of the detections, number of spurious sources which emerge in the process and number of real sources detected (detection limit and completitude magnitude). As we will see in this paper, the previous properties are not only related to the SNR gained in the filtering process but depend on the filtered momenta to 4th-order, the curvature of the source and the amplification in the 1D case. The combination of these quantities in a complicated way makes the decision on one filter totally dependent on the source profile and the background. This shows that, in addition to the amplification, one must take into account other quantities which also play an important role. The amplification was suggested by Vio et al. (2002) in order to compare filters. The final identification of the sources with a simple $`n\sigma `$ thresholding is not the whole story regarding optimality. Moreover, if the scale of the source needs to be estimated from the data, it must be identified ‘a posteriori’ on the SNR map when using matched filters, what introduces noise in the process. On the contrary, the $`\mathrm{𝑎𝑑𝑎𝑝𝑡𝑖𝑣𝑒}`$ filter allows one to get the sources straightforwardly on the filtered map.

In section 2, we introduce two useful quantities: number of maxima in a Gaussian background and define the detection problem. In section 3, we remark on the optimal statistic that allows one to define the region of acceptance, i. e. the confidence level of the detections. In section 4, we comment on the different filters to be compared (Mexican hat, matched and scale-adaptive filters). In section 5, we give some analytical and numerical results. In section 6, we describe the numerical simulations performed to test some theoretical aspects and give the main results. Finally, in section 7, we summarize the main results and applications of this paper.

## 2 Local peak detection

Let us assume a 1D background (e. g. one-dimensional scan on the celestial sphere or time ordered data set) represented by a Gaussian random field $`\xi (x)`$ with average value $`\xi (x)=0`$ and power spectrum $`P(q),q|Q|`$: $`\xi (Q)\xi ^{}(Q^{})=P(q)\delta _D(qq^{})`$, $`\xi (Q)`$ is the Fourier transform of $`\xi (x)`$ and $`\delta _D`$ is the 1D Dirac distribution. The distribution of maxima was studied by Rice (1954) in a pioneer article, the expected number density of maxima per intervals $`(x,x+dx)`$, $`(\nu ,\nu +d\nu )`$ and $`(\kappa ,\kappa +d\kappa )`$ is given by

$$n_b(\nu ,\kappa )=\frac{n_b}{\sqrt{2\pi }}\frac{\kappa }{\sqrt{1\rho ^2}}e^{\frac{1}{2(1\rho ^2)}(\nu ^2+\kappa ^22\rho \nu \kappa )},$$
(1)
$$n_b\frac{1}{2\pi \beta \sqrt{\rho }},\nu \frac{\xi }{\sigma _0},\kappa \frac{\xi ^{\prime \prime }}{\sigma _2},\beta \sqrt{\frac{\sigma _0}{\sigma _2}},\rho \frac{\sigma _1^2}{\sigma _0\sigma _2},$$

where $`\nu (\mathrm{},\mathrm{})`$ and $`\kappa (0,\mathrm{})`$ represent the normalized field and curvature, respectively, and $`n_b`$ is the expected total number density of maxima (i. e. number of maxima per unit interval $`dx`$). $`\sigma _n^2`$ is the moment of order $`2n`$ associated to the field.

If the original field is linear-filtered with a circularly-symmetric filter $`\mathrm{\Psi }(x;R,b)`$, dependent on $`2`$ parameters ($`R`$ defines a scaling whereas $`b`$ defines a translation)

$$\mathrm{\Psi }(x;R,b)=\frac{1}{R}\psi \left(\frac{|xb|}{R}\right),$$
(2)
we define the filtered field as

$$w(R,b)=𝑑x\xi (x)\mathrm{\Psi }(x;R,b).$$
(3)
Then, the moment of order $`n`$ of the linearly-filtered field is

$$\sigma _n^22_0^{\mathrm{}}𝑑qq^{2n}P(q)\psi ^2(Rq),$$
(4)
being $`P(q)`$ the power spectrum of the unfiltered field and $`\psi (Rq)`$ the Fourier transform of the circularly-symmetric linear filter.

Now, let us consider a Gaussian source (i. e. profile given by $`\tau (x)=e^{x^2/2R^2}`$) embedded in the previous background. Then, the expected number density of maxima per intervals $`(x,x+dx)`$, $`(\nu ,\nu +d\nu )`$ and $`(\kappa ,\kappa +d\kappa )`$, given a source of amplitude $`A`$ in such spatial interval, is given by

$`n(\nu ,\kappa |\nu _s)`$ $`=`$ $`{\displaystyle \frac{n_b}{\sqrt{2\pi }}}{\displaystyle \frac{\kappa }{\sqrt{1\rho ^2}}}`$ (5)
$`e^{\frac{(\nu \nu _s)^2+(\kappa \kappa _s)^22\rho (\nu \nu _s)(\kappa \kappa _s)}{2(1\rho ^2)}},`$
where $`\nu (\mathrm{},\mathrm{})`$ and $`\kappa (0,\mathrm{})`$, $`\nu _s=A/\sigma _0`$ is the normalized amplitude of the source and $`\kappa _s=A\tau _\psi ^{\prime \prime }/\sigma _2`$ is the normalized curvature of the filtered source. The last expression can be obtained as

$$\kappa _s=\nu _sy_s,y_s\beta ^2\tau _\psi ^{\prime \prime },\tau _\psi ^{\prime \prime }=2_0^{\mathrm{}}𝑑qq^2\tau (q)\psi (Rq).$$
(6)
We consider that the filter is normalized such that the amplitude of the source is the same after linear filtering: $`𝑑x\tau (x)\mathrm{\Psi }(x;R,b)=1`$.

The total number density $`n`$ in the presence of a local source is obtained integrating eqn. (5)

$$n=n_be^{\frac{\kappa _s^2}{2}}\left[1+B\left(\frac{\kappa _s}{\sqrt{2}}\right)\right],B(x)\sqrt{\pi }xe^{x^2}erfc(x).$$
(7)
## 3 Optimal statistics

We want to make a decision between filters based on $`\mathrm{𝑜𝑝𝑡𝑖𝑚𝑎𝑙𝑖𝑡𝑦}`$. We will assume that it includes the following properties: a) confidence level of the detections, b) number of spurious sources which emerge in the process and c) number of real sources detected. As we will see, the previous properties are not only related to the SNR gained in the filtering process but depend on the filtered momenta to 4th-order, the curvature of the source and the amplification (1D case).

We will distinguish two cases: I) Simple detection, consisting of detecting the presence of a signal $`s(x)`$ in a background and II) Simple measurement, consisting of detecting the signal $`s(x)=A\tau (x)`$ and measuring its unknown amplitude $`A`$ (we assume the profile is given) in the presence of a background.

### 3.1 Simple detection

First, let us establish the confidence level of the detections assuming $`\mathrm{𝑠𝑖𝑚𝑝𝑙𝑒}`$ detection (the amplitude of a source $`A`$ is given and so we can calculate $`\nu _s`$): let us consider a local peak in the 1D data set characterised by the normalized amplitude and curvature $`(\nu ,\kappa )`$, if $`H_0:p.d.f.p(\nu ,\kappa |\nu _s)`$ represents the $`\mathrm{𝑛𝑢𝑙𝑙}`$ hypothesis that the peak is a source with normalized amplitude $`\nu _s`$ given the data $`(\nu ,\kappa )`$, and $`H_1:p.d.f.p(\nu ,\kappa |0)`$ represents the $`\mathrm{𝑎𝑙𝑡𝑒𝑟𝑛𝑎𝑡𝑖𝑣𝑒}`$ hypothesis that the peak is a maximum of the background, we can associate to any region $`R_{}(\nu ,\kappa )`$ two errors

$$\alpha =_R_{}𝑑\nu 𝑑\kappa p(\nu ,\kappa |0),1\beta =_R_{}𝑑\nu 𝑑\kappa p(\nu ,\kappa |\nu _s),$$
(8)
$`\alpha `$ is the false alarm probability or confidence level of the detection (i. e. it represents the probability of interpreting noise as signal) whereas $`\beta `$ is the false dismissal probability or equivalently, $`1\beta `$ is the power of the detection (i. e. $`\beta `$ represents the probability of interpreting signal as noise). $`R_{}`$ is called the $`\mathrm{𝑎𝑐𝑐𝑒𝑝𝑡𝑎𝑛𝑐𝑒}`$ region.

Clearly, the previous probabilities can be obtained from the number density of maxima given by eqns. (157) as

$$p_b(\nu ,\kappa )p(\nu ,\kappa |0)=\frac{n_b(\nu ,\kappa )}{n_b},p(\nu ,\kappa |\nu _s)=\frac{n(\nu ,\kappa |\nu _s)}{n},$$
(9)
where $`n`$ and $`n_b`$ are the total number density in the presence or absence of a local source, respectively.

We will assume as decision rule the Neyman-Pearson one (C1): the acceptance region $`R_{}`$ giving the highest power $`1\beta `$, for a given confidence level $`\alpha `$, is the region

$$L(\nu ,\kappa |\nu _s)\frac{p(\nu ,\kappa |\nu _s)}{p(\nu ,\kappa |0)}L_{},$$
(10)
where $`L_{}`$ is a constant. So, the decision rule is expressed by the likelihood ratio: if $`LL_{}`$ the signal is present, whereas if $`L<L_{}`$ the signal is absent. Once we have assumed the previous decision rule for simple detection, one can calculate the ROC-curves (‘receiver operating curves’ in the signal processing jargon): $`\beta =\beta (\alpha ;\nu _s)`$.

Now, we will introduce the significance $`s(\nu _s)`$ of the detection

$$s^2\frac{[N_{signal}N_{nosignal}]^2}{\sigma _{signal}^2+\sigma _{nosignal}^2},$$
(11)
where in the numerator appears the difference between the mean number of peaks in the presence and absence of signal: $`(1\beta )N`$ and $`\alpha N`$, respectively. In the denominator appear the variances in the presence and absence of signal: $`\beta (1\beta )N`$ and $`\alpha (1\alpha )N`$, respectively. The last quantities have been calculated taken into account that, in the absence of signal, the probability of detecting locally peaks in $`q`$ of $`N`$ realisations of the background is given by the binomial distribution: $`P_q=\left(\begin{array}{cc}N& \\ q& \end{array}\right)\alpha ^q(1\alpha )^{Nq}`$. Therefore, the significance is given by the function

$$s^2(\alpha )\frac{(1\beta \alpha )^2}{\beta (1\beta )+\alpha (1\alpha )},$$
(12)
The next step is to assume a $`\mathrm{𝑐𝑟𝑖𝑡𝑒𝑟𝑖𝑜𝑛}`$ to define the optimal confidence level (C2): maximizing $`s`$ respect to $`\alpha `$ (or equivalently $`L_{}`$). In this sense, for each concrete filter we are able to get the best conditions (i. e. the maximum power corresponding to the optimal confidence level in the sense of maximizing the significance of the detection, see Allen et al. 2002) for the filtered data to be analysed. Clearly, for each filter we have found a unique region of acceptance $`R_{}`$ given the amplitude of the source $`A`$.

### 3.2 Simple measurement

In this case the signal has an unknown parameter, the amplitude, that is measured. First, we are interested in the detection of the signal and then in the estimation of its amplitude. Therefore, we calculate the a posteriori p.d.f. of all the possible values of the parameter $`A`$ given the data $`(\nu ,\kappa )`$. In absence of any a priori information about the p.d.f. $`p(A)`$ we will assume uniformity in the interval $`[0,\nu _c]`$, i. e. we can integrate over $`\nu _s`$: $`p(\nu ,\kappa )=\frac{1}{\nu _c}_0^{\nu _c}𝑑\nu _sp(\nu ,\kappa |\nu _s)`$. With this new p.d.f., the likelihood ratio is defined as

$$L(\nu ,\kappa )\frac{p(\nu ,\kappa )}{p(\nu ,\kappa |0)}=\frac{1}{\nu _c}_0^{\nu _c}𝑑\nu _sL(\nu ,\kappa |\nu _s),$$
(13)
We express the decision rule giving the region of acceptance $`R_{}`$ as (criterion C1)

$$R_{}:L(\nu ,\kappa )L_{},$$
(14)
where $`L_{}`$ is a constant. Once we have assumed the previous decision rule for simple measurement, one can calculate the ROC-curves $`\beta =\beta (\alpha )`$, where one must use $`p_b(\nu ,\kappa )`$ and $`p(\nu ,\kappa )`$, respectively. By introducing the significance $`s`$ of the detection by eqn. (11), we adopt the same C2 criterion: maximizing $`s`$ respect to $`\alpha `$, to obtain the optimal confidence level.

Moreover, the number of spurious (corresponding to background fluctuations) sources $`n_b^{}`$ and real detections $`n^{}`$ are obtained by integration of $`n_b`$ and $`n(\nu ,\kappa )\frac{1}{\nu _c}_0^{\nu _c}𝑑\nu _sn(\nu ,\kappa |\nu _s)`$ in the region of acceptance $`R_{}`$.

On the other hand, regarding the estimation of the amplitude $`A`$, it is natural to define the most probable value of $`A`$ (i.e. the amplitude $`A`$ corresponding to the value $`\nu _s`$ where the likelihood $`L(\nu ,\kappa |\nu _s)`$ takes its maximum) as being the result of the measurement and consider this value as the $`\mathrm{𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑑}`$ value of $`A`$.

## 4 The filters

### 4.1 The scale-adaptive filter (SAF)

The idea of a scale-adaptive filter (or optimal pseudo-filter) has been recently introduced by the authors (Sanz et al. 2001). By introducing a circularly-symmetric filter, $`\mathrm{\Psi }(x;R,b)`$, we are going to express the conditions in order to obtain a scale-adaptive filter for the detection of the source $`s(x)`$ at the origin taking into account the fact that the source is characterised by a single scale $`R_o`$. The following conditions are assumed: $`(1)w(R_o,0)=s(0)A`$, i. e. $`w(R_o,0)`$ is an $`\mathrm{𝑢𝑛𝑏𝑖𝑎𝑠𝑒𝑑}`$ estimator of the amplitude of the source. $`(2)`$ the variance of $`w(R,b)`$ has a minimum at the scale $`R_o`$, i. e. it is an $`\mathrm{𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡}`$ estimator. $`(3)`$ $`w(R,b)`$ has a maximum at $`(R_o,0)`$. Then, the filter satisfying these conditions is given by the equation

$$\stackrel{~}{\psi }(q)\psi (R_oq)=\frac{1}{acb^2}\frac{\tau (q)}{P(q)}\left[nb+c(na+b)\frac{dln\tau }{dlnq}\right],$$
(15)
$$a𝑑q\frac{\tau ^2}{P},b𝑑qq\frac{\tau }{P}\frac{d\tau }{dq},c𝑑qq^2\frac{1}{P}\left(\frac{d\tau }{dq}\right)^2.$$
(16)
where $`\tau `$ is the profile of the source in Fourier space ($`s(x)=A\tau (x)`$). Generically, $`\mathrm{\Psi }`$ is not positive (i. e. the name $`\mathrm{𝑝𝑠𝑒𝑢𝑑𝑜}\mathrm{𝑓𝑖𝑙𝑡𝑒𝑟}`$) and $`\mathrm{\Psi }`$ does not define a continuous wavelet transform. Moreover, the filter $`\mathrm{𝑎𝑑𝑎𝑝𝑡𝑠}`$ to the source profile, the background and the scale of the source, i. e. the name $`\mathrm{𝑎𝑑𝑎𝑝𝑡𝑖𝑣𝑒}`$ filter.

The previous equations have been used by the authors to obtain the adaptive filter for a Gaussian and an exponential profile (Sanz et al. 2001) and a multiquadric profile (Herranz et al. 2002a,b). The previous SAF has been recently used by Herranz et al. (2002c) for point source detection and extraction from simulated Planck time-ordered data.

Assuming a scale-free power spectrum, $`P(q)=Dq^\gamma `$, and a Gaussian profile for the source, the previous set of equations lead to the filter

$$\psi _a(q)=\frac{1}{\mathrm{\Gamma }(m)}q^\gamma e^{\frac{1}{2}q^2}\left[1t+\frac{t}{m}q^2\right],$$
(17)
$$m\frac{1+\gamma }{2},t\frac{1\gamma }{2}.$$

In Fig.1 appears the SAF for different values of the spectral index $`\gamma =0,1,2`$. In this case the filter parameters $`\beta `$ and $`\rho `$ and the curvature of the source $`y_s`$ are given by

$`\beta ^2=R^2\sqrt{{\displaystyle \frac{1+\frac{t^2}{m}}{m(1+m)\left[1+\frac{t^2}{m}+\frac{2t(2+t)}{m^2}\right]}}},y_s={\displaystyle \frac{\beta ^2}{R^2}},`$
$`\rho =\sqrt{{\displaystyle \frac{m}{1+m}}}{\displaystyle \frac{1+\frac{t^2}{m}+\frac{2t}{m^2}}{\sqrt{\left(1+\frac{t^2}{m}\right)\left(1+\frac{t^2}{m}+\frac{2t(2+t)}{m^2}\right)}}}`$ (18)
In Fig. 2 appear the parameters $`\beta \sqrt{\rho }`$, $`\rho `$ and $`y_s`$ as a function of the spectral index $`\gamma `$ for the SAF.

The identification of sources as peaks above a certain threshold (e. g. $`3\sigma _0`$) in filter space gives a low probability of false detections (reliability) because if the background has a characteristic scale of variation different from the sources, then practically everything detected with our method is real.

For comparison with the filter developed in this subsection, we shall briefly introduce other couple of filters that have been extensively used in the literature: the Mexican hat wavelet and the matched filter.

### 4.2 The Mexican Hat wavelet (MHW)

The MHW on $`\mathrm{}`$ is defined to be proportional to the Laplacian of the Gaussian function: $`\psi _h(x)(1x^2)e^{x^2/2}`$. Thus, in Fourier space

$$\psi _h(q)=\frac{2}{\pi ^{1/2}}q^2e^{\frac{1}{2}q^2}.$$
(19)
In Fig. 1 appears the MHW compared to other filters.

In this case the filter parameters $`\beta `$ and $`\rho `$ and the curvature of the source $`y_s`$ are given by

$$\beta ^2=\frac{R^2}{\sqrt{(2+t)(3+t)}},\rho =\sqrt{\frac{2+t}{3+t}},y_s=\frac{3\beta ^2}{2R^2}.$$
(20)
The parameters $`\beta \sqrt{\rho },\rho `$ and $`y_s`$ as a function of the spectral index $`\gamma `$ are given in Fig. 2 for the MHW.

We comment that the generalization of this type of wavelet for two dimensions has been extensively used for point source detection in 2D images.

It should be noted that the MHW studied in this work is constructed at a fixed scale $`R`$ given by the profile of the source. This differs from the filter used in Vielva et al (2001a,b) where an optimal scale is determined from the data and used to construct the MHW. This optimal scale is chosen to give maximum amplification of the source and thus the MHW at the optimal scale (MHWo) will give a higher gain than the MHW at the scale of the source. Moreover, Vielva et al. perform a multiscale fit in order to estimate the amplitude of the source.

### 4.3 The matched Filter (MF)

If one removes condition (3) defining the SAF in the previous section, it is not difficult to find another type of filter after minimization of the variance (condition (2)) with the constrain (1)

$$\stackrel{~}{\psi }_m(q)\psi (R_oq)=\frac{1}{2a}\frac{\tau (q)}{P(q)}.$$
(21)
This will be called $`\mathrm{𝑚𝑎𝑡𝑐ℎ𝑒𝑑}`$ filter as is usual in the literature. Note that in general the matched and adaptive filters are different.

For the case of a Gaussian profile for the source and a scale-free power spectrum given by $`P(q)q^\gamma `$, the previous formula leads to the following matched filter

$$\psi _m(q)=\frac{1}{\mathrm{\Gamma }(m)}q^\gamma e^{\frac{1}{2}q^2},m\frac{1+\gamma }{2}.$$
(22)
In Fig. 1 appears the MF for different values of the spectral index $`\gamma =0,1,2`$. We remark that for $`\gamma =1`$ the adaptive filter and the matched filter coincide and for $`\gamma =2`$ the matched filter and the Mexican Hat wavelet are equal.

For the MF the parameters $`\beta `$ and $`\rho `$ and the curvature of the source $`y_s`$ are given by

$$\beta ^2=\frac{R^2}{\sqrt{m(1+m)}},\rho =\sqrt{\frac{m}{1+m}},y_s=\rho .$$
(23)
In Fig. 2 appear the parameters $`\beta \sqrt{\rho },\rho `$ and $`y_s`$ as a function of the spectral index $`\gamma `$ for the MF.

## 5 Analytical and numerical results

We will distinguish the two cases of simple detection and simple measurement in the applications that follow.

### 5.1 Simple detection

In this case, we try to detect sources of known amplitude $`A`$. We note that the same amplitude $`A`$ translates into different thresholds $`\nu _{s,i}=A/\sigma _i`$ (where $`i`$ refers to the SAF, MF or MHW) in the filtered map for the three considered filters, since they have different amplifications. The relation between the values of $`\nu _s`$ for each filter can be easily obtained taking into account their relative amplification, which is given in Fig. 3. In the following comparison, we will consider the same amplitude of the source for all the filters and its value will be given as a function of the dispersion of the map filtered by the $`SAF`$.

Using eqns. (9,10) one obtains for the likelihood $`L(\nu ,\kappa |\nu _s)`$

$$L(\nu ,\kappa |\nu _s)=e^{\nu _s\phi \frac{\nu _s^2}{2}\mu }\left[1+B\left(\frac{\nu _sy_s}{\sqrt{2}}\right)\right]^1,$$
(24)
$$\phi (\nu ,\kappa )\frac{1\rho y_s}{1\rho ^2}\nu +\frac{y_s\rho }{1\rho ^2}\kappa ,\mu \frac{(1\rho y_s)^2}{1\rho ^2}.$$
(25)
The region of acceptance $`R_{}`$ is defined by $`LL_{}`$ or equivalently $`\phi \phi _{}(\nu _s)`$ (we assume $`\nu _s0`$) with

$$\phi _{}=\frac{1}{\nu _s}\mathrm{ln}\left[L_{}\left(1+B\left(\frac{\nu _sy_s}{\sqrt{2}}\right)\right)\right]+\frac{\nu _s}{2}\mu $$
(26)
Clearly, the last constraint can be rewritten

$$R_{}:\nu \nu _{}(\kappa ;\phi _{}),\nu _{}\frac{\rho y_s}{1\rho y_s}\kappa +\frac{1\rho ^2}{1\rho y_s}\phi _{}.$$
(27)
In order to find the acceptance region, we need to maximize the significance $`s^2`$ with respect to $`\alpha `$ (or equivalently $`L_{}`$) for each filter. Fig. 4 shows the curve of $`s^2`$ versus $`\alpha `$ and versus $`L_{}`$ for the three filters with $`\nu _{s,\mathrm{SAF}}=3`$ and $`\gamma =0,1.4,2.2`$. It is remarkable that for all the cases the maximum of $`s^2`$ is found at $`L1`$ within a few per cent. This is also true independently of the source amplitude (for the studied range of $`\nu _{s,\mathrm{SAF}}`$ from 1 to 5). Therefore the chosen criterion indicates that we are in the best conditions to discriminate between the two hypothesis by using a simple acceptance region: the candidates are accepted as detections if the probability of having background plus source is greater or equal than having only background. The value of $`\alpha `$ that maximizes the significance depends mainly on the amplitude of the source and ranges from $`0.3`$ for $`\nu _{s,\mathrm{SAF}}=1`$ to $`0.01`$ for the highest amplitude considered, $`\nu _{s,\mathrm{SAF}}=5`$. This simply reflects the fact that the higher the amplitude of the source, the lower the number of spurious detections and vice versa. Taking into account the previous results we have calculated the acceptance region using $`L_{}=1`$ for all the cases studied in this section.

In the bottom right panel of Fig. 5 the acceptance region for the case $`\gamma =1.4`$ and $`\nu _{s,\mathrm{SAF}}=3`$ is given for the different filters. Those maxima with curvature and amplitude above the line are accepted as point sources, while those below are rejected. The slopes of the lines are easily understood by looking at the other three panels of the figure, that show the curvature distribution given a threshold $`\nu =\nu _s`$ of maxima of the filtered background (thin line) and of the background plus source (thick line). For the MF both distributions are identical, i.e., we can not differentiate between true and spurious detections based on their curvature. Therefore, the acceptance region is fixed by the amplification of the filter: if the maximum is found above a certain threshold it is accepted, otherwise is rejected. For the SAF, the distribution of curvature is shifted to higher values when a source is present. Thus, maxima with larger curvature are accepted at lower thresholds than those with smaller curvature, since the latter are more likely to be produced by the background. On the contrary, when the field is filtered with the MHW, the maxima produced by the source tend to have a smaller curvature than those of the background and the slope of the curve changes with respect to the SAF case.

Once $`R_{}`$ has been obtained, we can calculate the density number of spurious sources (corresponding to background fluctuations) $`n_b^{}`$ and real detections $`n^{}`$ by integrating eqns. (1,5) in the region of acceptance:

$`n_b^{}`$ $`=`$ $`{\displaystyle \frac{n_b}{2}}[\mathrm{erfc}\left({\displaystyle \frac{\sqrt{1\rho ^2\phi _{}}}{\sqrt{2}(1\rho y_s)}}\right)`$
$`+\sqrt{2}ay_s\mathrm{e}^{a^2\phi _{}^2}\mathrm{erfc}({\displaystyle \frac{\sqrt{1\rho ^2}}{1\rho y_s}}y_sa\phi _{})],`$
$`a`$ $`=`$ $`\sqrt{{\displaystyle \frac{1\rho ^2}{2(12\rho y_s+y_s^2)}}},`$ (28)
$`n^{}`$ $`=`$ $`{\displaystyle \frac{n_b}{2}}{\displaystyle _0^{\mathrm{}}}𝑑\kappa \kappa e^{\frac{1}{2}(\kappa y_s\nu _s)^2}\mathrm{erfc}(z),`$
$`z`$ $`=`$ $`{\displaystyle \frac{\sqrt{1\rho ^2}}{\sqrt{2}(1\rho y_s)}}\left[\phi _{}y_s\kappa \nu _s{\displaystyle \frac{(1\rho y_s)^2}{1\rho ^2}}\right]`$ (29)
The number densities $`n_b^{},n^{}`$, the ratio $`r=n^{}/n_b^{}`$, $`\alpha `$ and $`1\beta `$ have been plotted in Fig. 6 versus the spectral index $`\gamma `$ for a source amplitude corresponding to $`\nu _{s,\mathrm{SAF}}=3`$ for the SAF, MF and MHW. In order to compare the performance of the different filters, we have also plotted the relative difference $`D`$ of the detection ratio with respect to the MF, which is defined as

$$D=\frac{r_ir_{\mathrm{MF}}}{r_{\mathrm{MF}}}\times 100$$
(30)
where the subindex $`i`$ refers to the different filters. Therefore a positive value of $`D`$ indicates that the corresponding filter has a better detection ratio than the MF. Three different regions can be seen in relation to the performance of the filters. For $`\gamma 01`$ the MF clearly outperforms the SAF and the MHW. The SAF works better than the other two filters in the range $`\gamma 11.6`$. Finally, at the highest values of $`\gamma `$ the MHW has the best performance of the three considered filters, although it is only slightly better than the MF. This behaviour is qualitatively similar for other source amplitudes, although the performance of the MF in comparison with the other filters is improved for bigger amplitudes and gets worse for smaller ones (see Fig. 7).

Although a detailed comparison with the MHWo is beyond the scope of this work, we have also plotted the detection ratio relative to the MF obtained by the MHWo for the case $`\nu _{s,\mathrm{SAF}}=3`$ using the region of acceptance defined by $`L1`$ (see middle right panel of Fig. 6). The MHWo clearly outperforms the MHW at the scale of the source in the considered range of spectral indexes. In fact, this filter gives the best detection ratio for high values of $`\gamma `$.

### 5.2 Simple measurement

In this case, the likelihood ratio is given by

$$L(\nu ,\kappa )=\frac{1}{\nu _c}_0^{\nu _c}𝑑\nu _se^{\nu _s\phi \frac{\nu _s^2}{2}\mu }\left[1+B\left(\frac{\nu _sy_s}{\sqrt{2}}\right)\right]^1.$$
(31)
As in the previous case, we need to maximize the significance with respect to $`L_{}`$ in order to get the acceptance region $`LL_{}`$. The values of $`L_{}`$ that maximize the significance are plotted versus the spectral index in Fig. 8 for $`\nu _{c,\mathrm{SAF}}=5`$. In this case, the values of $`L_{}`$ are in the range of $`1.5`$ to $`1.8`$ depending on the filter and spectral index. It also depends on the assumed value of $`\nu _c`$, higher values of this parameter lead to larger values of $`L_{}`$. This has an effect on the values of $`L_{}`$ for the different filters, since the same amplitude of the source translates into different thresholds in the filtered maps and therefore we have a different $`\nu _c`$ for each filter. Thus the relative difference in the value of $`L_{}`$ between filters seen on Fig. 8 are closely related to the relative amplification (see Fig. 3). The corresponding acceptance regions for $`\gamma =0,1.4,2.2`$ are given in Fig. 9.

The density number of spurious sources (corresponding to background fluctuations) $`n_b^{}`$ and real detections $`n^{}`$ are obtained by integration of $`n_b(\nu ,\kappa )`$ and $`n(\nu ,\kappa )\frac{1}{\nu _c}_0^{\nu _c}𝑑\nu _sn(\nu ,\kappa |\nu _s)`$ in the region of acceptance $`R_{}`$ respectively.

The quantities $`n_b^{}`$, $`n^{}`$, $`n^{}/n_b^{}`$, D, $`\alpha `$ and $`1\beta `$ are given in Fig. 10 as a function of the spectral index. The acceptance region has been obtained using the value of $`L_{}`$ (given in Fig. 8) obtained from the maximization of the significance in each case. We see that the situation is qualitatively similar to the simple detection case regarding the relative difference in the detection ratio. Again the MF outperforms the other filters at low $`\gamma `$, the SAF has the best behaviour at intermediate values whereas the MHW gives a higher ratio at large $`\gamma `$.

On the other hand, a minimization of the likelihood ratio $`L(\nu ,\kappa |\nu _s)`$, given by eqn. (24), allows an estimation of the amplitude of the source $`\widehat{\nu }_s`$ to be obtained as a function of the data $`(\nu ,\kappa )`$. The estimation of the amplitude is then given as the solution of the equation

$$(y_s^2+\mu )\nu _s\phi +\frac{1}{\nu _s}\frac{B\left(\frac{\nu _sy_s}{\sqrt{2}}\right)}{1+B\left(\frac{\nu _sy_s}{\sqrt{2}}\right)}=0,$$
(32)
where the function $`B`$ is given by eqn. (7) and $`\mu `$ and $`\phi `$ by eqn. (25). A confidence level for the estimation of the amplitude of the source can also be obtained from eqn. (24).

## 6 NUMERICAL SIMULATIONS: RESULTS

Let us consider how the ideas presented in the previous sections apply to real cases. In order to do this, we simulated a set of one-dimensional ‘images’ containing a Gaussian background characterised by a power spectrum $`P(q)q^\gamma `$. Two different experiments have been carried out. In the first set of simulations Gaussian sources of known width and amplitude were introduced. This case corresponds with the *simple detection* scheme presented in section 3.1, in which the knowledge of the amplitude of the source (or, equivalently, its normalized amplitude $`\nu _s`$) can be used to determine the region of acceptance. In the second set of simulations the Gaussian sources had amplitudes that were drawn from a uniform distribution between a minimum and a maximum value. These simulations were used to test the *simple measurement* scenario, in which the question is not only whether the source is detected or not, but also if its amplitude can be estimated as well.

### 6.1 Simulations in the simple detection scenario

In order to test the performance of different filters in the simple detection scenario, we simulated a set of ‘images’ with $`N=1024`$ pixels each and a Gaussian background characterised by a power spectrum $`P(q)q^\gamma `$. For each simulated image, a source with a Gaussian profile of width $`FWHM=5`$ pixels and known amplitude $`A_o`$ was placed at the central pixel. Then, the image was filtered with a matched filter, a scale-adaptive filter and a Mexican Hat wavelet. The matched and scale-adaptive filters were chosen taking into account the source width and the background index as in eqns. (15) and (21), whereas the Mexican Hat had the same width than the Gaussian profile of the source. The normalized amplitude of the filtered source is $`\nu _{s,X}=A_o/\sigma _X`$, where $`X`$ refers to the filter in consideration (SAF, MF or MHW). After filtering, we look for a maximum at the position of the source we introduced. Here we implicitly assume that we know the exact position of the source, which is true for our simulations but will not be the case in a realistic observation. Fixing the position of the source allows us to clearly illustrate the behaviour of the optimal statistics described above. A short comment on position uncertainties will be presented in section . If there exists such a maximum, we measure its normalized amplitude $`\nu `$ and curvature $`\kappa `$. With these quantities and the previously known $`\nu _s`$ it is easy to determine whether the maximum lies inside a certain region of acceptance $`R_{}`$ or not. We adopt the criterion C1 introduced in section 3.1: $`R_{}:L(\nu ,\kappa )L_{}`$, and choose $`L_{}=1`$ (for a justification for this particular value see section 5). If the maximum is in the region of acceptance, we count it as a *valid detection* and, if not, as a *rejection* (a false dismissal). The number of rejections is directly related with the probability $`\beta `$ of eqn. (8).

On the other hand, let us consider a simulation with the same background, but without any source, and filter it in the same way as before. Now, let us repeat the same procedure than before: we look for a local maximum at the central pixel of the ‘image’ and, in case it exists, we apply the selection mechanism, that is, make the (wrong) hypothesis that a source of normalized amplitude $`\nu _s`$ is present in the simulation, and determine the region of acceptance as before. If the maximum is not inside the region of acceptance, it has been safely rejected. On the contrary, if it lies inside the region of acceptance, the noise will be interpreted as a signal and a false detection will occur. The number of false detections is directly related with the probability $`\alpha `$ of eqn. (8).

In order to test the behaviour of the probabilities of detecting a signal, rejecting a signal and obtaining a false detection as functions of the source amplitude $`A_o`$ and the filter, we tested 10 different values of $`A_o`$. Three different background regimes were tested: background index $`\gamma =0`$, corresponding to a case in which the MF is expected to detect with better reliability than the other two filters, $`\gamma =1.4`$, corresponding to a case in which the SAF is the most reliable, and $`\gamma =2.2`$, where the most reliable filter is the MHW (see section 5). The lower and higher amplitude limits $`A_o^{min}`$ and $`A_o^{max}`$ were chosen so that after filtering with the SAF the normalized amplitudes are $`\nu _{s,SAF}^{min}1`$ and $`\nu _{s,SAF}^{max}5`$. For each value of $`A_o`$ we performed 10000 simulations containing a source in the central pixel and other 10000 simulations containing only background. Each simulation was filtered using the SAF, the MF and the MHW. For each pair of filtered simulations (one with source and one without it) the following quantities were recorded: the number of times that a source was properly detected $`N_d`$, the number of times a source produced a maximum but it was rejected by the selection criterion $`N_r`$ and the number of spurious detections (due to the background) that happened in the simulations without source, $`N_b`$. The total number of maxima found when a source is present is $`N_t=N_d+N_r`$. Another interesting quantity is the ratio $`r=N_d/N_b`$, that gives us an idea of the practical reliability of the detection.

#### 6.1.1 The case $`\gamma =0`$

The results of the simulations with $`\gamma =0`$ are shown in figure 11. The results of the simulations, indicated with points, are compared with the theoretical curves obtained in section 5. Solid lines and filled circles make reference to the scale-adaptive filter. Dotted lines and open circles make reference to the matched filter. Dashed lines and asterisks make reference to the Mexican Hat wavelet.

The first panel on the top and the left of figure 11 shows the number of detections $`N_d`$ (local maxima due to sources that are inside the region of acceptance) in the 10000 simulations for the three filters. $`N_d`$ is directly related with the number density of maxima $`n^{}`$, being the relationship between the two quantities dependent on the sizes of the pixel and the filter scale $`R`$. Since the matched filter has the greatest gain, it gives the highest number of detections. The number of detections found in the simulations is lower than the expected value. This result is not unexpected due to the following reason: the theoretical expected number of detections is calculated by multiplying the number density of maxima in the presence of a local source $`n^{}`$ by the number of simulations (10000 in this case) and the pixel size (in units of the filter width $`R`$). This product can, eventually, give more than one theoretical detection per simulation and pixel, if the pixel size and $`\nu _s`$ are big enough. Obviously, in a real simulation one can only find zero or one maxima per pixel. This effect is more conspicuous for high $`\nu _s`$ and, as we will see later, for higher $`\gamma `$ values, producing apparent paradoxes such as the prediction, in some cases, of a number of detections greater than the number of simulations we have. Therefore, the theoretical curve must be considered in the best case only as a upper limit for the real number of detections.

The panel on the top and the right of figure 11 shows the number of spurious detections $`N_b`$ found in the 10000 simulations versus the expected value. Again, $`N_b`$ is directly proportional to the number density of maxima due to the background $`n_b^{}`$. The dots were obtained by counting the number of accepted background maxima considering all the pixels in absence of local source and then computing the average number of spurious detections per pixel, and the statistical $`1\sigma `$ error bars of such average. The agreement between the expected values and the results from the simulations is remarkable. The matched filter gives the lowest gross number of spurious detections.

The left and the right lower panels in figure 11 show the ratio between the number of background maxima interpreted as source ($`N_b`$) and the total number of maxima due to the background, and the ratio between source detections ($`N_d`$) and the total number of maxima in presence of a source ($`N_t=N_d+N_r`$), respectively. The first ratio gives us an estimation of the false alarm probability $`\alpha `$ whereas the second one gives the power of the detection $`1\beta `$. The results obtained with the SAF and the MF are strikingly similar to the expected values, whereas the concordance is slightly worse in the case of the MHW. Note that the problem of the theoretical curves mentioned above does not affect these ratios because the terms that come from the number of simulations and the pixel area cancel in the division.

More interesting is the ratio between the number of effective detections and the number of spurious ones. This quantity is shown in the middle panel on the left of figure 11. Since the proportionality between $`N_d`$ and $`n^{}`$ is the same than the one between $`N_b`$ and $`n_b^{}`$, this ratio is a good estimator of the ratio $`r=n^{}/n_b^{}`$. The inverse of this ratio gives the percentage of spurious detections that is found in a number $`N`$ of detected maxima. Again, the agreement between theory and simulations is excellent. Note that the greatest values of $`r=N_d/N_b`$ correspond to the matched filter. That means that, if what we are looking for is *reliability* in the detection, the matched filter is the best choice in this case $`\gamma =0`$.

In order to quantify the differences between filters regarding reliability (in the sense mentioned above), it is useful the quantity $`D`$ defined in equation (30). This $`D`$ can be understood as the per cent relative difference between filters. In the middle panel on the right of figure 11 it is shown the value of $`D`$. The differences between the MF and the other filters are $`50\%`$ for the SAF and $`80\%`$ for the MHW at intermediate $`\nu _s`$ values.

Gaussian backgrounds with $`\gamma =0`$ correspond to uncorrelated (white) noise and are very common in many fields of physics. For instance, the background of many astronomical images is dominated by uncorrelated Poissonian noise, which is well approximated by this kind of noise when averaging over many observations.

#### 6.1.2 The case $`\gamma =1.4`$

The results of the simulations with $`\gamma =1.4`$ are shown in figure 12. As in the case before, the results from the simulations are compared with the theoretical curves obtained in section 5. The different lines and panels are the same as those in figure 11. The agreement between the theoretical expectations and the results is very good.

The number of detections $`N_d`$ at a given normalized amplitude $`\nu _s`$ is higher for the three filters than in the case $`\gamma =0`$. The filtering is more efficient at enhancing sources when $`\gamma `$ increases due to the fact that at high $`\gamma `$ the power of the background concentrates in large scale structures that are more easily removed by the filters.

In this case the performance of the three filters is similar in number of detections as well as in the probabilities $`\alpha `$ and $`1\beta `$. The SAF, however, produces a lower number of spurious detections due to the background. Therefore, the ratio $`r`$ between authentic and spurious detections is better for the case of the SAF (although the MF performs almost as well in this aspect). The relative difference $`D`$ between SAF and MF in this case takes values around $`10\%`$. The MHW performs significantly worse than the other two filters in this case.

We remark that noise index in the interval $`1\gamma 1.6`$ are not rare in some areas of Astronomy. For example, in the scanning of the sky of many CMB experiments (MAP, Planck, etc) backgrounds with noises in this range may appear due to the combination of CMB and Galactic foregrounds with the scanning $`1/f`$ noise.

#### 6.1.3 The case $`\gamma =2.2`$

The results of the simulations with $`\gamma =2.2`$ are shown in figure 13. As in the two cases before, the results from the simulations are compared with the expectations obtained in section 5. The different lines and panels are the same as those in figure 11.

This case corresponds with the region in which the MHW outperforms the other two in reliability. The matched filter, however, performs almost as well as the MHW. The relative differences $`D`$ in the detection ratio are of a few percent between the MHW and the MF. The SAF performs significantly worse in this case.

The agreement between the results of the simulations and the theoretical predictions is good again, though not so good as in the previous cases. The reason is that in practice it is not possible to simulate perfectly a $`P(q)q^\gamma `$ background with $`\gamma >0`$ because such backgrounds diverge when $`q`$ tends to $`0`$. The discrepancy between the ideal power spectrum and the simulated one gets worse when $`\gamma `$ increases.

Spectral indices around 2 correspond to smooth backgrounds and are less common than smaller indices. However, they appear in certain cases. For instance, microwave observations dominated by dust emission have $`\gamma 2`$.

#### 6.1.4 Determination of the position

In the simulations presented above, our aim was reproducing the theoretical scheme using simulations. Therefore, we have focused only on what happens in one pixel of the image, the pixel in which the source is located. However, in most real situations, the position of the source is not known and we will need to consider all the pixels of the image.

Note that we are assuming that the position of the maximum is a good estimator of the position of the source. In practice, this is not necessaraly true. In particular, the maximum can be shifted from the position of the source due to background fluctuations, although the probability of finding a peak due to source plus background decreases quickly as the distance to the source position increases. We have calculated the probability of a peak appearing at a certain distance $`d_p`$ (in pixel units) from the position of the simulated source for the case $`\gamma =1.4`$. We expect the results to be qualitatively similar for the other cases. At $`\nu _s=5.0`$ we found that the probability of the peak appearing at distance $`d_p=1`$ is approximately $`0.2`$ for the three considered filters. At $`\nu _s=3.0`$ this probability increases to $`0.4`$, whereas at $`\nu _s=1.0`$ it increases to $`0.5`$. When we consider $`d_p=2`$ pixels we found these probabilities to be practicaly zero for cases $`\nu _s=5`$ and $`\nu _s=3`$, and 0.2 for $`\nu _s=1`$. For $`d_p3`$ the probability is almost zero for all cases. Therefore, in a realistic situation where we need to consider all the pixels of the image, and in particular the neighbouring pixels of the source, the number of detections will slightly increase (in a very similar manner for all the considered filters). Regarding the number of spureous sources, they do not depend on the pixel position. Thus the spureous detections will increase linearly with the considered number of pixels. This behaviour is again the same for all the filters. Therefore the main conclusion is that while the global performance of the filters will be worse in a realistic case in which the position of the sources is not known, the relative performance of the filters will still be the same as that described in the previous sections.

Note that in our simulations the FWHM of the beam is 5 times the pixel size and that basically all maxima produced by a point source fall within the beam size. Therefore, by simply using the position of the maximum as an estimator of the position of the source, this last quantity can be determined with a precision given by the beam size of the experiment. More ellaborated techniques can be used to improve the determination of the position. For instance one could use the average of the points at half the peak intensity along the slopes of the image, or fitting these points to a certain profile. This is out of the scope of the present paper and remains for a future work.

### 6.2 Simulations in the simple measurement scenario

Now consider that the true amplitude of the source is unknown. This is the most frequent case in practise. In such a case, one must both detect and estimate the amplitude of the sources that are supposed to be embedded in the data. This problem can be confronted by adopting the methodology considered at the end of section 3.2: in absence of any a priori information about the p.d.f. $`p(A)`$ we assume uniformity in the interval $`[0,\nu _c]`$, where $`\nu _c`$ is some cut threshold over which it is not expected to find sources. An approximate value of $`\nu _c`$ can be obtained directly from the data under consideration simply by adopting $`\nu _c\left(|I_{max}|+|I_{min}|\right)/\sigma `$, where $`I_{max}`$ and $`I_{min}`$ are the maximum and minimum data values, respectively. In practise, it is very unlikely that a source with normalized amplitude greater than the previous value is present in the data. If the p.d.f. $`p(A)`$ is known, that information can be used instead of the previous uniform distribution. As explained in section 3.2, using the p.d.f. $`p(A)`$ it is possible to determine a region of acceptance for the detections and the normalized likelihood $`L(\nu ,\kappa |\nu _s)`$. By means of the normalized likelihood one can estimate the most probable value and confidence limits (‘error bars’) for the amplitude $`A`$ of any detected peak.

To test this methodology a set of simulations similar to the ones described in section 6.1 were performed. Half of the simulations had a Gaussian source of $`FWHM=5`$ pixels located at the central pixel, whereas the other half had only noisy background. Now we try to detect sources in both kind of simulations without knowing a priori the amplitude of the (hypothetical) source. In each case, a region of acceptance is determined as in eqns. (1314). If a peak is found and it lies inside the region of acceptance, its most probable amplitude is estimated by looking for the maximum value of the normalized likelihood $`L(\nu ,\kappa |\nu _s)`$ and the $`68\%`$ confidence limits (‘error bars’) are calculated using the normalized likelihood $`L(\nu ,\kappa |\nu _s)`$. The number of ‘detections’ found in the simulations without source indicates the probability of spurious detection in this scenario.

The amplitude of the simulated sources takes values between 0 and a maximum value $`A_{max}`$ so that $`A_{max}/\sigma _{SAF}=5`$, being $`\sigma _{SAF}`$ the average value of $`\sigma _0`$ over 1000 simulations filtered with the SAF. We take this value as the value of $`\nu _c`$ for the calculations. For the sake of simplicity and without any loss of generality, we consider 25 amplitude bins in the interval $`[0,\nu _c]`$ instead of a continuous sampling. A number of 5000 simulations with source and 5000 simulations without source were performed for every amplitude bin. As in the previous section, three cases were considered: background indexes $`\gamma =0`$, $`\gamma =1.4`$ and $`\gamma =2.2`$.

#### 6.2.1 The case $`\gamma =0`$

In the simple measurement scenario we need to take two aspects into account. On the one hand, there is the pure detection aspect of the process, that is, how many sources are detected, and how many spurious detections due to background are accepted by our decision criterion. On the other hand, there is the issue of how well the parameters of the sources (i.e. their amplitude) are estimated.

The three first rows of table 1 show the relevant results for the simulations with $`\gamma =0`$ and the three filters. The theoretical expectation for the number of simulations we performed is shown in parenthesis. The number of detections is slightly lower in the simulations than in the theoretical prediction, as we have seen in section 6.1, and for the same reasons that were explained there. The number of spurious detections $`N_b`$ was slightly higher than the theoretical expectation. The lower number of detections and higher number of spurious sources lead to ratios $`r=N_d/N_r`$ that are worse than in the theoretical case. The difference, however, is never greater than a $`20\%`$, and the theoretical behaviour is preserved: the MF gives the best ratio $`r`$, followed by the SAF and then by the MHW. The agreement in the power of $`1\beta `$ of the filters is very good for the three filters.

Regarding the estimation of the amplitudes, the results are shown in figure 14. The filled circles show the mean value of the estimated normalized amplitude as a function of the true normalized amplitude of the source. The estimation has been performed using eqn. (32). As a reference, the mean values of the directly measured amplitudes are shown by the open circles. Clearly, the estimation based in the maximization of the likelihood $`L(\nu ,\kappa |\nu _s)`$ is more accurate than the naive estimation using directly the measured values of the maxima. Error bars show the dispersion of the estimates around the mean value in the simulations (solid line for the estimation using eqn. (32) and dotted line for the naive estimation).

For high amplitudes the estimation of the amplitude is unbiased for the three filters. At low amplitudes there is a positive bias due to a selection effect: the sources are very faint and only in the cases when a background fluctuation enhances by chance the source the associated maximum is able to enter in the acceptance region $`R_{}`$. Using eqn. (32) we are able to estimate reasonably well amplitudes as low as $`2.5\sigma _{\mathrm{SAF}}`$.

The lower panel on the right of figure 14 show the mean values of the estimated amplitudes using eqn. (32) for comparison between the three filters. SAF is denoted by filled circles, whereas open circles and asterisks stand for MF and MHW, respectively. The three filters show a very similar efficiency at high amplitudes, whereas at low amplitudes the MF and the SAF work better than the MHW.

As noted in section 3.2, using the normalized likelihood $`L(\nu ,\kappa |\nu _s)`$ it is possible to give a confidence interval for the estimation of the amplitude. To illustrate this point, figure 15 shows the estimated amplitude and the $`68\%`$ confidence intervals calculated using eqn. (24) for the three filters in two different cases: when the source is in the limit of detection (left panels of the figure) and when the source has a good SNR after filtering (right). Only 30 examples were randomly chosen for each filter and amplitude bin in order to make the plot readable. For comparison, the mean value and statistical error bar that appears in figure 14 is shown at the right of each plot. The confidence intervals are significantly larger than the statistical error bars. Hence, they must be considered as upper limits to the true error. However, in a realistic case they are the only thing that one can safely say about the error distribution without having to resort to simulations.

#### 6.2.2 The case $`\gamma =1.4`$

The rows 4 to 6 in table 1 shows the results for the simple measurement scenario and $`\gamma =1.4`$. As expected, the total number of detections $`N_d`$ has increased with respect to the case $`\gamma =0`$, while the number of spurious detections $`N_b`$ has decreased. The relative difference between both numbers and their theoretical expectations are below the $`13\%`$ in all the cases. The lowest number of spurious detections corresponds to the SAF. Moreover, the best ratio $`r=N_d/N_b`$ corresponds to the SAF as well. The differences between the predicted $`r`$ and the ratio given by the simulations is below the $`3\%`$ for the SAF and the MF, and below the $`6\%`$ in the case of the MHW. The agreement in the probability $`1\beta `$ is very good for the three cases.

Figure 16 shows the results of the estimation of the amplitude in the case $`\gamma =1.4`$. The meaning of the panels, points and lines is the same than in figure 14. Again, the estimation of the amplitude using eqn. (32) is better than the simple estimation using directly the measured intensity of the maxima. The three filters perform very similarly in this aspect. As in the previous case, the effective limit for a good estimation of the amplitude is around $`\nu _{s,\mathrm{SAF}}2.5`$. We would like to remark, however, that since the filters produce greater amplification in this case, the absolute (no normalised) amplitude one can detect and correctly estimate is in this case lower than in the case $`\gamma =0`$.

The estimation of confidence intervals works exactly in the same way as before.

#### 6.2.3 The case $`\gamma =2.2`$

The rows 7 to 9 in table 1 shows the results for the simple measurement scenario and $`\gamma =2.2`$. The total number of detections $`N_d`$ is greater than in previous cases, that is, the amplification (gain) of the filters is greater than for lower $`\gamma `$ values. In the same way, the number of spurious detections $`N_b`$ is smaller than before. The relative difference between both numbers and their theoretical expectations are below the $`15\%`$ in all the cases. As expected, the lowest number of spurious detections and the best ratio $`r`$ correspond to the MHW (but the MF performs almost equally well). Again, there is a remarkable agreement in the probability $`1\beta `$ between simulations and theoretical expectations.

Figure 17 shows the results of the estimation of the amplitude in the case $`\gamma =2.2`$. There are not qualitative differences with respect to the previous two cases. Note once more that having the same normalized amplitude $`\nu `$ than in the cases $`\gamma =0,1.4`$ does not mean having the same amplitude at the beginning. As the gain of the filter grows, we can reach fainter sources.

## 7 Conclusions

Filtering is a useful technique for the denoising and enhancing of compact signals in Astronomy and many other applications. However, the diversity of different kind of filters makes difficult to choose which one is the most convenient in any particular case. In this paper we have addressed the problem of the choice of filter for the detection and measurement of point sources in a noisy background. We focus on the one-dimensional case and assume that the sources have a Gaussian spatial profile and the noise can be modelled by an isotropic and homogeneous Gaussian random field. Then we compare three different filters that have been thoroughly used in the study of Cosmic Microwave Background radiation data: the Mexican Hat wavelet, the matched filter and the scale-adaptive filter. Although we have focused in these three filters, the methodology we describe can be easily applied to any linear filter.

Two main scenarios have been explored: the case in which one tries to detect a source of known amplitude but uncertain existence and location in the data (*simple detection*) and the case in which one tries to detect sources and determine their a priori unknown amplitude (*simple estimation*). In both cases, local peak detection is the first step to be undertaken. In the one-dimensional case, any local peak is characterised by its intensity (amplitude) $`\nu `$ and its curvature $`\kappa `$. Once a peak is detected, a *decision* about its validity must be done.

We use the Neyman-Pearson criterion to define a region of acceptance in the $`(\nu ,\kappa )`$ space that maximizes the *significance* of the detections. The limits of such region will depend on the properties of the background (namely, its power spectrum) and the filter that has been applied to the data.

In the simple measurement scenario, the Neyman-Pearson criterion is applied by calculating the a posteriori p.d.f. of all the possible values of the amplitude given the data $`(\nu ,\kappa )`$. In absence of any a priori information about the p.d.f. of the amplitudes of the sources $`p(A)`$, we assume uniformity between $`0`$ and a certain cut $`\nu _c`$ that can be easily guessed from the data. If we had a priori information about $`p(A)`$, it could be straightforwardly included in the formalism. Once the source is detected, its amplitude $`A`$ can be estimated by maximizing the likelihood $`L(\nu ,\kappa |\nu _s)`$ with respect to $`\nu _s`$. Convenient confidence intervals can be given for any detection using the previous normalized distribution.

We have compared the MHW, the SAF and the MF in the two scenarios presented above. Both analytical and empirical comparisons have been performed considering a generic background characterised by a power spectrum $`P(q)q^\gamma `$. In the analytical comparison we have shown how to derive useful formulae to predict the number density of local peaks that lie inside the Neyman-Pearson region of acceptance both in presence and in absence of source. With these quantities it is possible to compute the ratio between the number of true detected sources and the number of spurious detections, that is, a measure of the *reliability* of each filter in the task of detecting sources. We find that, regarding this last quantity, there are three different regions demarcated by the value of the index $`\gamma `$: for $`0\gamma <1`$ the MF outperforms the other two in the reliability sense as well as in total number of true detections. For $`\gamma =1`$ the MF and the SAF coincide. For $`1<\gamma 1.6`$ the SAF gives the best reliability. The relative difference with the MF is greater than a $`10\%`$. For $`\gamma =2`$ the MF and the MHW coincide. Finally, for $`\gamma 1.6`$ the MHW is the most reliable filter for the detection, although the matched filter performs almost equally well. However, the performance of the MHW can be improved by using it with an optimal scale (that can be estimated from the data) which gives maximum amplification of the source. These three regions are present and their limits are roughly the same both in the simple detection and in the simple measurement scenarios.

We have performed exhaustive simulations in order to test the previous ideas. The cases $`\gamma =0`$, $`1.4`$ and $`2.2`$ have been chosen with the aim of exploring the three regions described above. The results of the simulations totally agree with the analytical expected behaviour.

Regarding the estimation of sources of unknown amplitude, the three filters perform equally well. The estimation of the amplitude using the likelihood of the data is fairly better than the estimation using directly the measured intensity of the peaks. At low source amplitudes there is a positive bias due to a selection effect. With these filters and the likelihood estimator, it is possible to safely reach thresholds as low as $`2.5\sigma _{\mathrm{SAF}}`$ in the filtered data. Due to the amplification effect of the filters, the equivalent thresholds in the unfiltered maps can be really low, specially for high $`\gamma `$ indexes. As an example, in the case $`\gamma =2.2`$ the MF produced in our simulations a mean gain of $`7.5`$, this meaning that the ‘safe threshold’ of $`2.5`$ after filtering translates into a $`0.33`$ threshold before filtering.

The ideas presented in this paper can be generalized to more general filtering schemes and for two-dimensional data sets. In this case, other quantities such as the ellipticity of the peaks can be useful to establish decision criteria similar to the ones presented here. A future work will address this particular issue.

## Acknowledgements

The authors thank Patricio Vielva for useful discussions. RBB thanks the Ministerio de Ciencia y Tecnología and the Universidad de Cantabria for a Ramón y Cajal contract. DH acknowledges support from a UC postdoctoral fellowship during the months of March to October 2002 and from the European Community’s Human Potential Programme under contract HPRN-CT-2000-00124, CMBNET, since November $`1^{st}`$ 2002. We acknowledge partial support from the Spanish MCYT projects ESP2001-4542-PE and ESP2002-04141-C03-01 and from the EU Research Training Network ‘Cosmic Microwave Background in Europe for Theory and Data Analysis’.