A theoretical framework for the detection of
pointsources in Cosmic Microwave Background maps
Key Words.:
Methods: data analysis – Methods: statistical – Cosmology: cosmic microwave backgroundAbstract
Context:
Aims:The detection of pointsources in experimental microwave maps is a critical step in the analysis of the Cosmic Microwave Background (CMB) data. If not properly removed, these sources can have adverse effects on the estimation of the powerspectrum and/or the test of Gaussianity of the CMB component. In the literature, various techniques have been presented to extract point sources from an observed image but no general consensus about their real performance and properties has been reached. Their characteristics have been studied essentially through numerical simulations based on semiempirical models of the CMB and the Galactic foreground. Such models often have different levels of sophistication and/or are based on different physical assumptions (e.g. the number of Galactic components and level of the noise). Moreover, the application of a given technique to a set of data (either simulated or experimental) requires the tuning of one or more parameters that unavoidably is a subjective operation. Hence, a reliable comparison is difficult. What is missing is a statistical analysis of the properties of the proposed methodologies. This is the aim of the present paper.
Methods:The statistical properties of the detection techniques in the context of two different criteria, i.e. the NeymanPearson criterion and the maximization of the signaltonoise ratio (SNR), are analyzed through an analytical approach. Onedimensional as well twodimensional signals are considered. The case of multiple observing frequencies is also addressed.
Results:The conditions are fixed under which the techniques can work satisfactorily. Their limits are also illustrated and implementation details provided. We show that, exploiting some a priori information, it is possible to develop simple algorithms with performances not too far from those of more sophisticated but complex techniques. In this respect, a detection algorithm, tailored for high Galactic latitudes, is presented that could be useful in future groundbased experiments as, for example, the Atacama Large Millimeter/submillimeter Array (ALMA).
Conclusions:
1 Introduction
The detection of pointsources embedded in a noise background is a critical issue in the analysis of the experimental Cosmic Microwave Background (CMB) maps. The estimation of the powerspectrum of the CMB component and the test of its possible nonGaussian nature need the a priori detection and removal of these sources. In particular the former operation is rather delicate. Given its importance, this subject has been extensively considered in literature (see Herranz and Sanz 2008a, and references therein). However, no widespread accepted conclusion has been reached. The reason is that most of the detection techniques presented in literature lack a sufficiently rigorous theoretical background. The statistical characteristics are derived from numerical experiments only. Especially in the experiments where it is necessary to simulate maps at different observing frequencies, this approach is hazardous for two reasons: 1) the models used by the various authors in the simulations are not the same; 2) when a detection technique has to be applied to a set of data (either synthetic or real) the tuning of some parameters is unavoidably subjective. Because of this contrasting results appear in literature. A safer procedure consists in studying the techniques in a well defined theoretical framework. Although some of the conditions assumed to fix the theoretical context could be not realistic, anyhow it is possible to obtain a more objective comparison as well as some useful indications that help to understand what it is possible to expect when they are relaxed. The aim of the present paper is to provide a theoretical characterization of some of the detection techniques suited for both single and multiplefrequency CMB experiments. For ease of notation, initially the arguments will be developed for onedimensional signals . Later they will extended to the twodimensional situation. The paper is organized as follows: in Sec. 2 the problem of pointsource detection is described in detail for the case of a single spatial dimension. Here some standard material is presented with the aim to fix notation and formalism. In Sec. 3 the situation is considered when more signals are available. In Sec. 4 a technique tailored for observation at high Galactic latitude is presented that will be useful in some experiments planned in the near future with innovative instruments as the Atacama Large Millimeter/submillimeter Array (ALMA). Finally, the conclusions are given in Sec. 5. In appendix A, two procedures are presented that efficiently implement the techniques described in the text. These are extended to the twodimensional case in appendix B.
2 Formalization of the problem
The first step in the development of a detection technique is to fix the conditions that are pertinent to the problem of interest. In the case of CMB observations, the following conditions are commonly assumed:

The pointsources have a known spatial profile . The amplitude “” is a scalar quantity different from source to source, whereas is a function, due to the instrument beam, that is identical for all of them. Function is normalized in such a way that ;

The pointsources are embedded in a noisebackground, i.e. the observed signal is given by . In other words, noise is additive. Usually, this is a reasonable assumption;

Noise is the realization of a stationary stochastic process with known covariance matrix
(1) Actually, because of the Galactic contribution, especially at low Galactic latitudes, this hypothesis is not satisfied. However, it is assumed to hold locally. This allows the computation of statistics as the mean or the covariance matrix that, otherwise, should not be possible. Without loss of generality, it is assumed that .
Under these conditions, the detection problem consists in deciding whether is a pure noise (hypothesis ) or it contains also the contribution of a source (hypothesis ). In other words, the source detection problem is equivalent to a decision problem where two hypotheses hold:
(2) 
Under the probability density function of is given by whereas under by . At this point, it is necessary to fix the criterion to use for the detection. Clearly, one cannot hope to find all the sources present in a given signal. Hence, some choices are necessary. For example, one could decide that the nondetection or the misidentification of a bright source could be more important than that of a fainter one, or vice versa. A very common and effective criterion is the NeymanPearson criterion that consists in the maximization of the probability of detection under the constraint that the probability of false alarm (i.e., the probability of a false detection) does not exceed a fixed value . The NeymanPearson theorem (e.g., see Kay 1998) is a powerful tool that allows to design a decision process that pursues this aim: to maximize for a given , decide if the likelihood ratio (LR)
(3) 
where the threshold is found from
(4) 
The test of the ratio (3) is called the likelihood ratio test (LRT).
An important example of application of LRT is when noise is Gaussian with correlation function . Actually, in CMB experiments this condition is satisfied only for observations at high Galactic latitudes where the CMB emission and the instrumental noise are by far the dominant contributors. At lower latitudes, it is often assumed to hold locally. For example, the contribution to of components that in small patch of sky presents linear spatial trends are often approximated with stationary Gaussian processes with a steep spectrum (e.g. noises). In any case, even if the Gaussianity condition was unrealistic, it is often made anyway since it allows an analytical treatment of the problem of interest and the results can be used as a benchmark in the analysis of more complex scenarios. With Gaussian , it is
(5)  
(6) 
with
(7) 
The LRT is given by
(8) 
Hence, it results that has to be chosen when for the statistic (called NP detector) it is
(9) 
with such as
(10) 
i.e.,
(11) 
Here, with the standard Gaussian distribution function, is the corresponding inverse function. Equation (10) is due to the fact that is a Gaussian random variable with variance and expected values equal to zero under and under . For the same reason it is
(12) 
Equation (9) can be written in the form
(13) 
with
(14) 
From this equation appears that can be thought as a linear filter of signal , that is called matched filter (MF).
2.1 Some comments on the use of the matched filter in practical applications
There are some important points to stress with regard the MF when used in practical aplications. They are:

is a sufficient statistic (Kay 1998). Loosely speaking, this means that is able to summarize all the relevant information in the data concerning the decision (2). No other statistic can perform better. As a consequence, the claim that some filters (e.g. the Mexican hat wavelet, the scaleadaptive filters, the biparametric scaleadaptive filters…) are superior to MF according to the NeymanPearson criterion (Barreiro et al. 2003; LópezCaniego et al. 2005) is not correct. It is the result of the use of imprecise theoretical arguments (e.g. see Vio et al. 2004);

It is worth noticing that if the amplitude “” of the source is unknown, then Eq. (9) can be rewritten in the form
(15) with . In other words, a statistic is obtained that is independent of “”, i.e. also in the case that the amplitude of the source is unknown, maximizes for a fixed . The only consequence is that cannot be evaluated in advance. In principle this quantity could be evaluated a posteriori by using the maximum likelihood estimate of the amplitude, , but this is of little interest. More useful is that in real experiments one is typically interested in the detection of sources which have amplitudes characterized by a probability density function . In this case, once is fixed to a value and making to change “” across the domain of , the quantity , with as given by Eq. (12) and , provides an estimate of the fraction of undetected sources as function of their amplitude;

If and , then
(16) with any invertible linear operator (matrix). A useful consequence of this property is that if signal is convolved with a function (e.g., the beam of an instrument), this operation does not modify the optimality of MF. This fact could be useful in situations where more signals are available that are obtained with different point spread functions (see below);

The arguments above are developed under the implicit assumption that sources do not overlap and that their position is known in advance. In practical applications the first condition is satisfied – strictly speaking – only in very high resolution maps and it is assumed to be always valid for those sources above the confusion noise. The second statement, i.e. that the position of the sources is known in advance, is not true and the standard procedure consists in filtering by means of and in computing for the peaks in the resulting signal. Although there is no guarantee that a peak in the filtered signal marks the true position of a source even in the case this is effectively present, theoretical arguments as well as years of application in reallife problems have proved that this procedure is rather robust and able to provide excellent results;

If the Gaussianity of is relaxed, then MF is no longer optimal in the NeymanPearson sense. However, it remains optimal with respect to the signaltonoise ratio (SNR) ^{1}^{1}1Here, the quantity SNR is defined as the ratio between the squared amplitude of the filtered source with the variance of the filtered noise.. This means that, independently of the nature of the noise, MF provides the greatest amplification of the signal with respect to the noise. This can be easily verified through the minimization of the variance of the filtered noise with the constraint that (i.e. filter does not modifies the amplitude of the source), in formula ^{2}^{2}2We recall that the functions “” and “” provide the values of of for which the function has the smallest and greatest value, respectively.
(17) with a Lagrange multiplier. It is not difficult to see that
(18) (19) i.e. apart from a normalizing factor, is given by Eq. (14). Since is optimal as concerns the maximization of the SNR, no other methods can outperform it in this respect. For this reason, expedients as the introduction of additional constraints and/or of free parameters in (e.g. see Sanz et al. 2001; Herranz et al. 2002) has the only effect to reduce the detection performances (e.g. see Vio et al. 2002). One of the benefit in using is that the value of provides directly an unbiased estimate of the amplitude “” of the source. However, it is necessary to stress that in practical applications, where the true position of the source is not known and it is necessary to apply the procedure described above, this is no longer true;

When is a long signal, some computational problems come out. In fact, if the size of is , then is a matrix. Hence, the computation of the quantity can become quite expensive. The computational burden can be alleviated if , that is a Toeplitz matrix, is approximated with a circulant matrix . This is because can be diagonalized through
(20) Here, is the onedimensional Fourier matrix that is a complex, unitary, and symmetric matrix whose elements are given by
(21) Symbol “” denotes the onedimensional Fourier transform, , is the complex conjugate transpose of , and is a diagonal matrix containing the eigenvalues of (i.e., the powerspectrum of ). After that, since is a unitary matrix, one obtains that
(22) (23) Symbol “” denotes the elementwise multiplication, and a column vector containing the diagonal elements of the square matrix . In obtaining this result, the fact that has been used. In principle, all the terms in Eq. (23) could be computed quite efficiently by means of fast Fourier transform (FFT). In particular, array can be obtained by means of the reciprocal of the FFT of the autocorrelation function . Actually, in using , there is the problem that both and are implicitly assumed to be periodic signals with period . In general, this is not true. Hence, boundary effects are to the expected in the computation of the FFT. However, since typically has finite spatial support, these effects can be easily avoided by padding it with a sequence of leading zeros longer than the correlation length of the noise. This is visible in the last three panels in Fig. 1 where the array , that approximates the MF in Eq. (14), is shown when is constructed with the correlation function shown in the first panel and is a rectangular function (shown in the same panel) that is padded with an increasing number of leading zeros. When this number is sufficiently large, it is evident that . In principle, the same method could be applied to . However, because of the noise, usually this signal has no finite spatial support. Hence, in order to evaluate , often one wishes to compute the inverse FFT of , to remove a number of leading elements from the resulting array equal to that of the zeros used in the padding operation, and then to calculate the scalar product with .

In the case is a colored noise (i.e. is not a diagonal matrix), then the length of and should be longer than the correlation length of . This is clearly visible in Fig. 2 where the different performances of MF are compared when is the realization of a Gaussian process with a correlation length of about pixels and is a Gaussian with , dispersion set to three pixels that is computed, respectively, on , and pixels. It is evident that when the performance of MF becomes independent of this parameter. The comparison is based on the so called receiver operating characteristics (ROC) that is a plot of versus . This kind of plot is very useful since it allows a direct visualization of the detection performances of a given technique. More specifically, the ROC should always be well above a straight lineline since this corresponds to a detection performance identical to that of flipping a coin, ignoring all the data.
3 Extension to the mutiplefrequency case
In the context of CMB observations, there is a further complication in that there are signals , , , coming from the same sky area that are taken at different observing frequencies. Here, is the amplitude of the source at the th observing frequency, whereas is the corresponding spatial profile. For ease of notation, all the signals are assumed to have the same length . In general, the amplitudes as well as the profiles are different for different . However, if one sets
(24)  
(25)  
(26) 
it is possible to obtain a problem that is formally identical to that treated in the previous section. Hence, the MF is still given by Eqs. (13)(14) and is named multiple matched filter (MMF). The only difference with the classic MF is that now is a block matrix with Toeplitz blocks (BTB):
(27) 
i.e. each of the blocks is constituted by a Toeplitz matrix. In particular, provides the autocovariance matrix of the th noise, whereas , , the crosscovariance matrix between the th and the th ones.
In spite of these similarities, when some difficulties arise. In particular, cannot be written in a form equivalent to Eq. (15). This has important consequences in the fact that if the amplitudes are unknown, then cannot be computed. In other words, if the spectral characteristics of the radiation emitted by a source are not fixed, then the MMF is not applicable. This forces to resort to an approach based on the maximization of the total SNR of the filtered signals. Following this approach, model (17) has to be modified in the form
(28) 
where is an array, , and is a matrix
(29) 
with a array. Now, since and with and a diagonal matrix whose diagonal contains , it is trivial to show that Eq. (28) is equivalent to
(30) 
where and
(31) 
The solution is
(32) 
It is evident that, contrary to MMF, with it is possible to obtain a statistic ,
(33) 
that is independent of the unknown source amplitude . The price to pay is a reduced detection capability. In fact, better results should be obtainable if the SNR of each signal was maximized. However, this is an operation that requires the knowledge of . For , it is
(34) 
that again is a quantity independent of the source amplitude, whereas as expected the same is not true for
(35) 
In two recent works Herranz and Sanz (2008a) and Herranz et al. (2008b) have proposed a “new” class of filters, the so called matrix filters. Their idea is to filter separately each of the signals in such a way to obtain unbiased estimates of and at the same time to simultaneously minimize the total variance of the filtered signals. In the spatial domain, this problem can be written in the form
(36) 
with “” the Trace operator,
(37) 
a matrix of filters and
(38) 
a matrix of Lagrangian multipliers. The solution is
(39) 
Through it is possible to define a set of statistics
(40) 
that can be considered individually. However, the fact that
(41) 
indicates that and essentially represent the same filter. The only difference is that, after filtering signals , the latter does not compose them together as the former does. We call the modified multiple matched filter (MMMF).
In appendix A two procedures are presented that work in the Fourier domain and that allow a very efficient computation of both and . In appendix B the methods are extend to the twodimensional case.
3.1 Alternative techniques
A strategy to deal with multiplefrequency images that is alternative to the approaches presented above consists in composing signals together in a single array . After that, the classic MF could be applied. Without a priori information, the most obvious method is
(42) 
where is an operator (matrix) such as independent of “”. In this way, the effects due to the different instrumental beams can be avoided. As seen above, the use of does not represent a problem since it does not modify the statistic . This approach, that we indicate as summedimage matched filter (SMF), can be expected to provide results close to those of the MMF when the amplitudes as well as the level of the noises are similar. However, if this condition is not satisfied, its performances rapidly worsens especially if and , , have some degree of correlation. This last condition is typical of CMB signals where with a component independent of (i.e. the CMB contribution), and a whitenoise process due to the electronic of the instrument. In this case, it could be preferable a weighted sum as
(43) 
where
(44)  
(45) 
with . The first constraint implies that the contribution of in is completely removed, whereas the second one provides a normalizing factor. We indicate this method as weighted matched filter (WMF) and the particular case where with , as uniformly weighted matched filter (UWMF). This last corresponds to a situation where only one signal is used to eliminate the shared noise component , whereas the others are given an identical weight. The UWMF is not very effective, but it can be useful in absence of a priori information on the spectral properties of the sources.
The performance of these methods depends critically by many factors as, for example, the relative value of the amplitudes , the correlation lengths of , the relative importance of the noises as well as the degree of correlation between the different observing frequencies. This fact is evident in Fig. 3 that shows the vs. the amplitude when and and , respectively. Here , is a Gaussian with dispersion set to three pixels, with a zeromean, unitvariance Gaussian process whose autocorrelation function has Gaussian profile and dispersion set to ten pixels and finally , , two independent Gaussian whitenoise processes. Two different cases are considered for the noises . In the first the noises and have the same dispersion, i.e. , whereas in the second one and . When , the only possible weights for the WMF are either . This is the case considered by Chen and Wright (2008).
One indication that comes out from these examples is that, as expected, MMF outperforms all the other methods. Moreover, unless the level of noises are similar, MMMF outperforms SMF. Heuristically, these results can be explained by the fact that through MMF a sum of signals is computed that is weighted by means of both the intensities and the noise levels , whereas with MMMF the weighting is based on the noise levels only. In SMF there is no weighting at all. A different situation is for WMF. Here the two signals are subtracted. This operation provides a signal for which i) the correlated part in the noises is zeroed; ii) the total amplitude is ; iii) the variance of the instrumental noise is . Since , it is not difficult to realize that a benefit in the detection capability happens when . In the case case , similar arguments indicate that good results can be expected when a channel “” is available with . This is the case expected in CMB applications (see below).
4 An approach for CMB experiments at high Galactic latidude
In the near future some innovative groundbased experiments are planned for very high spatial resolution observations as, for example, with the Atacama Large Millimeter/submillimeter Array (ALMA). One important advantage of these experiments is that, contrary to the satellite observations, they will allow a certain control of the experimental conditions. Hence, a full exploitation of the capabilities of this facility (and other instruments) requires a careful planning of the observations. The recording of good quality data will allow an effective application of the chosen techniques for the detection of pointsources. For instance, with instruments as ALMA it is possible to plan observations at high latitude fields to map sources at very high spatial resolution in sky regions dedicated to CMB observations. In this case a detection technique is necessary that takes into account such a specific experimental situation.
In the context of pointsource detection, data can be thought as twodimensional discrete maps , each of them containing pixels, corresponding to different observing frequencies (channels), with the form
(46) 
Here, correspond to the contribution of the pointsources at the th frequency, whereas denotes the corresponding noise component. At high Galactic latitudes, the CMB component is expected to be the dominant one. Hence, may be modeled with
(47) 
where is the experimental noise corresponding to the th channel and the contribution of the CMB component that is frequencyindependent. The contribution of the pointsources is assumed in the form
(48) 
with the amplitude of the source to the th channel. According to Eq. (48), and without loss of generality, all the sources are assumed to have the same profile independently of the observing frequency. In fact, although in general this will not be true, as written in Sect. 2.1, it is possible to meet this condition by convolving the images with an appropriate kernel with no consequences. In the following it is assumed that the components are the realization of stationary, zeromean, stochastic processes.
The main feature of model (46)(47) is that the CMB contribution does not change with frequency. Hence, following the suggestion in Sec. 3.1, the WMF can be used. Of course, in order this approach be effective, it is necessary that the amplitude of a given pointsource is not the same in all the observing channels.
In the case , (i.e. only two maps are available), the only possible solution is or . However, for more degrees of freedom are available. This allows the selection of the weights in such a way that specific conditions are satisfied. In particular, one could wish that, after the linear composition of the maps, the quantity
(49) 
is maximized, i.e.
(50) 
Here, is the crosscovariance matrix of the noise processes whose th entry is given by , with the operator that transforms a matrix into a column array by stacking its columns one underneath the other. The rationale behind this choice is the the quantity is a measure of the amplitude of the pointsource in the weighted map with respect to the standard deviation of the measurement noise. The larger the more prominent is the pointsource with respect to the noise background; an attractive situation in problems of source detection. The maximization, via the Lagrange multipliers method, of with the constraint (44) provides the following system of nonlinear equations
(51) 
It can be solved through an iterative algorithm based on the Newton method
(52) 
where
(53) 
and
(54) 
with
(55)  
(56)  
(57)  
(58) 
Here, three points are of concern: a) the iteration can be initialized with a starting guess that contains random entries satisfying the constrains (44)(45); b) The normalization term in Eq. (52) implements the constraint (45). This cannot be done via a Lagrange multiplier since the array that maximizes the quantity can be determined unless a multiplicative constant. Hence, the Lagrange multiplier corresponding to the constraint (45) should be equal to zero; c) For the same reason, matrix is rank deficient, hence has to be understood as MoorePenrose pseudoinverse.
Once that the map has been produced, the MF can be used with no necessity to take into account the characteristics of the CMB. This is particularly useful in situations where only small patches of sky are available and hence the sizes of are much shorter than the correlation length of . As indicated earlier, we call this method the weighted matched filter (WMF).
As for MMF, the procedure described above needs that the array be specified. In fact, if the weights are computed for a source with a given , they will be optimal only for any other source with amplitude . However, it is the general trend that does matter. For example, satisfactory results can be expected for the sources that present steep spectra with similar behaviors (see below). This means that the optimization of can be carried out for subsets of sources for which the above condition approximately applies.
4.1 Numerical experiments
In this section we present some numerical experiments to test the performances of WMF. In particular, we consider a scenario where three different observing frequencies are available. We make the simplifying assumption that all the channels have the same pointspread function (PSF) which is a twodimensional circular symmetric Gaussian normalized to have a peak value equal to one and with a dispersion set to three pixels. In correspondence to the th channel the contribution due to a pointsource is , whereas the terms in the covariance matrices are given by
(59) 
Here is the variance of the instrumental noise, assumed of Gaussian whitenoise type, and is the covariance matrix of the CMB component sampled with a step of on a regular twodimensional grid. For , it is . This scenarios mimics that expected for the ”LowFrequency Instrument” mounted on the PLANCK satellite (Vio et al. 2003). The available data are assumed in form of square maps containing pixels each ^{3}^{3}3In practical CMB applications, the fact to working with square patches of sky is not a limit since, independently of the shape of the available maps, pointsource detection is typically carried out on small spatial windows sliding across the sky. This is for computational reasons as well as because the noise contaminating the maps has no uniform spatial characteristics.. This size is large enough, with respect to the correlation length of , to make results independent of it. For a power law spectrum as S= the amplitude of a pointsource at an observed frequency , given its amplitude at an observed frequency , expressed in Thermodynamic temperature can be written as:
(60) 
where
(61) 
(62) 
with the Planck constant, the Boltzmann constant, the CMB temperature, and is the spectral index.
When expressed in antenna temperature this latter takes
the value for the infrared sources and for the radio ones (dominated by synchrotron emission).
Note that this is an approximated expression, usually used in CMB experiment
to make a direct comparison between the sky temperature brightness and the sources brightness.
Fig. 4 compares the ROC for MMF and WMF for, respectively, the radio and the infrared sources case,
in a situation of relatively low SNR and with the level of noise that is the same for all the channels.
The amplitudes , computed through Eq. (60), correspond to the , , GHz observing frequencies.
We allow a 10% deviation from the amplitudes given in (60). Hence, in the same figures
the confidence envelopes are shown for both MMF and WMF that have
been obtained by computing the ROC for a set of one hundred arrays , ,
a Gaussian random array with mean zero and covariance matrix . For reference, the result obtainable with the UWMF is also plotted.
The improved performances of WMF with respect to UWMF is evident. Moreover, although as expected
the MMF is always superior to WMF, their performances are rather similar. The reason is that,
in the case of sources with steep spectra as given by Eq. (60), there are at least two observing frequencies, say “” and “”,
for which . As seen at the end of Sec. 3, this is the condition for WMF to work well.
5 Conclusions
In this work, the problem of pointsource detection in noise background has been addressed from a theoretical perspective. This allowed an objective comparison of the expected performance of the various techniques. From this comparison it is evident that “in se” no method is superior to the other ones. The differences are due to the amount of a priori information that they exploit. In other words, it appears that the effectiveness of a technique is not linked to its sophistication rather to the ability in using the available a priori information. In particular, the methods based on the NeymanPearson criterion can be expected to provide better performances than those based on the maximization of the signaltonoise ratio because they make use of the amplitude of the sources. For the same reason, an improvement in the detection capability can be expected when more signals are available that correspond to the same sky area taking at different observing frequencies. On the other hand, the a priori information that at high Galactic latitudes the dominant components are the electronic noise and the CMB, with this last independent of the observing frequency, suggests that a suboptimal approach based on an opportune linear combination of the signals that eliminates the CMB contribution could have a performance close to MF. Hence, it is useful in practical applications as it avoids the estimation of the covariance matrix (or the powerspectrum) of the CMB.
As last remark we would like to add the following. The superiority of a theoretical approach does not diminish the usefulness of the numerical experiments. It is, however, a bad habit to fix the characteristics of a statistical methodology only by means of numerical simulations. The application of a detection technique to a set of (either real or synthetic) data requires the tuning of some parameters that is a subjective operation. As shown in Sec. 2.1, MF requires that the position of the candidate source is known in advance, while in practical application this piece of information lacks. The common solution consists to filter with and then to compute the statistics for the peaks in the resulting signal, but, because of noise, there is no guarantee that a peak in the filtered signal marks the true position of a source even in the case this is effectively present. In a numerical experiment this implies the definition of a window, around the true position of the source, where a peak is assumed to identify a source candidate. The size of the window is a parameter that can have important consequences but there is not an objective criterion to fix it. The same holds also for the other techniques. Moreover, as remarked again in Sec. 2.1, (see also Fig. 2), the performance of MF depends on the relative size of the spatial region where the signal is sampled with respect to the correlation length of the noise. Hence the comparison of MF with other filters can give different results according to size of the patch of sky that are considered. For these reasons it is risky the use of classes of filters as the mexican hat wavelet family (LópezCaniego et al. 2006a, b) that lack any theoretical justification and whose effectiveness is supported only by numerical experiments.
References
 Barreiro et al. (2003) Barreiro, R.B., Sanz, J.L., Herranz, D., & MartinezGonzález, E. 2003, MNRAS, 342, 119
 Chen and Wright (2008) Chen, X., & Wright, E.L. 2008, ApJ, 681, 747
 Herranz et al. (2002) Herranz, D. & Sanz, J.L., Barreiro, R.B., & MartinezGonzález, E. 2002, ApJ, 580, 610
 Herranz and Sanz (2008a) Herranz, D. & Sanz, J.L. 2008a, arXiv:0808.0300v1
 Herranz et al. (2008b) Herranz, D., LópezCaniego, M., Sanz, J.L., & GonzálezNuevo, J. 2008b, arXiv:0808.2884v1
 Kay (1998) Kay, S.M. 1998, Fundamentals of Statistical Signal Processing: Detection Theory (London: Prentice Hall)
 LópezCaniego et al. (2005) LópezCaniego, M., Herranz, D., Barreiro, R.B., & Sanz, J.L. 2005, MNRAS, 359, 993
 LópezCaniego et al. (2006a) LópezCaniego, M. et al. 2006a, MNRAS, 369, 1063
 LópezCaniego et al. (2006b) LópezCaniego, M. et al. 2006b, MNRAS, 370, 2047
 Sanz et al. (2001) Sanz, J. L., Herranz, D., & MartinezGonzalez, E. 2001, ApJ, 552, 484
 Vio et al. (2002) Vio, R., Tenorio, L., & Wamsteker, W. 2002, A&A, 391, 789
 Vio et al. (2003) Vio, R. et al. 2003, A&A, 401, 389
 Vio et al. (2004) Vio, R., Andreani, P., & Wamsteker, W. 2004, A&A, 414, 17
Appendix A Efficient numerical implementations of the MMMF: one dimensional case
Filters and given, respectively, by Eqs. (32) and Eq.(39) can be efficiently computed in the Fourier domain. If in the BTB matrix given by Eq. (27) the Toeplitz blocks are approximated with circulant matrices , then each of them can be diagonalized by means of the matrix
(63) 
with “” the Kronecker product and the identity matrix. In this way a matrix
(64) 
is obtained where are diagonal blocks containing the eigenvalues of . The quantities and , can be computed firstly by solving the system of linear equations , and then through the product (or through the procedure described in Sec. 2.1). Since matrix is highly structured and sparse, the computational load is not excessive. This approach is suited to be implemented in highlevel programming languages as MATLAB that allow a friendly handling of arrays and matrices.
A more efficient algorithm, but suited for lowlevel programming languages as C or FORTRAN, is obtainable rearranging the elements of , , and according to the so called row rollout order, i.e.,
(65) 
with , and similarly for , and . After that, models. (30) and (36), as well as the corresponding solutions (32) and (39), still holds with and replaced, respectively, by
(66) 
where , and
(67) 
with
(68) 
Here, denotes the down circulant shifting operator that circularly down shifts the elements of a column array by positions. Now, it can be shown again (see Kay 1998, page 504) that, if one sets
(69) 
then
(70) 
where is a block diagonal matrix
(71) 
with
(72) 
, . Here, represents the crosspowerspectrum at frequency between and . Similarly, for a given signal , the array provides the corresponding FFT,
(73) 
From these considerations, it results that
(74) 
with
(75) 
Here, the advantage is represented by the fact that
(76) 
i.e., the inversion of , that is a large matrix with size , can be obtained through the inversion of a number of much smaller blocks. This fact, coupled with the structure of , allows to compute efficiently through blockmatrix operations. Equation (74) provides the discrete version of the result by Herranz and Sanz (2008a) and Herranz et al. (2008b). Finally, as for the spatial domain, it is
(77) 
Appendix B Extension of MF and MMMF to the twodimensional case
b.1 Singlefrequency observations
The extension of MF to the twodimensional signals and is conceptually trivial. If one sets
(78)  
(79)  
(80) 
formally the same problem is obtained as that given by Eq. (2). However, again some computational issues come out. The question is that, even for moderately sized signals, matrix becomes rapidly huge. In fact, if contains pixels, then is a matrix. Similarly to the onedimensional case, the computational burden can be alleviated by resorting to the Fourier domain. If is a rectangular map, then is a matrix of type block Toeplitz with Toeplitz blocks (BTTB), and it can approximated with a matrix of type block circulant with circulant blocks (BCCB). In fact, a BCCB matrix can be diagonalized through
(81) 
where , and is again a diagonal matrices containing the eigenvalues of . The good news is that these eigenvalues can be obtained through the application of the operator to the twodimensional FFT of the autocovariance function of . Since also is a unitary matrix, it is
(82)  
(83) 
Here, symbol “” now indicates the twodimensional FFT. Similar problems and solutions as in Sec. 2.1 hold concerning the fact that, with the use of , both and are implicitly assumed to be periodic functions along each dimension with period and , respectively. Once computed, filter , that is a array, can be converted in its original twodimensional form simply columwise reordering its elements in a matrix (i.e., the inverse operation of the operator).
b.2 Multiplefrequency observations
In the case of multifrequency observations, the situation becomes even worst since MF has to be applied to signals at the same time. Again, through
(84)  