SANEPIC: A MapMaking Method for Timestream Data From Large Arrays
Abstract
We describe a mapmaking method which we have developed for the Balloonborne Large Aperture Submillimeter Telescope (BLAST) experiment, but which should have general application to data from other submillimeter arrays. Our method uses a Maximum Likelihood based approach, with several approximations, which allows images to be constructed using large amounts of data with fairly modest computer memory and processing requirements. This new approach, Signal And Noise Estimation Procedure Including Correlations (SANEPIC), builds upon several previous methods, but focuses specifically on the regime where there is a large number of detectors sampling the same map of the sky, and explicitly allowing for the the possibility of strong correlations between the detector timestreams. We provide real and simulated examples of how well this method performs compared with more simplistic mapmakers based on filtering. We discuss two separate implementations of SANEPIC: a bruteforce approach, in which the inverse pixelpixel covariance matrix is computed; and an iterative approach, which is much more efficient for large maps. SANEPIC has been successfully used to produce maps using data from the 2005 BLAST flight.
1 Introduction
The problem of optimal “mapmaking” or “image reconstruction” is complex and multifaceted, with the basic procedures and even the terminology differing dramatically between different subfields of astronomy. The method adopted depends on the form in which the data are gathered, and on the dominant source of systematic effects. From the optical to the nearIR one talks about combining “frames”, along with measurements of “darks” and “flatfields”. For the reduction of Cosmic Microwave Background (CMB) data the now conventional method is to start from the principle that there is a linear algebra approach to solving the Maximum Likelihood problem. However, this has only been feasible up until now because of the limited number of detectors in the typical CMB experiment, and the fact that correlated signals among the detectors can be effectively ignored. Because of the rapid development of large bolometer arrays, the question which arises is: how does one adapt the CMB approach to dealing with substantial numbers of detectors, and where there are significant crosscorrelations of noise between the detector timestreams.
Data from the Balloonborne Large Aperture Submillimeter Telescope (BLAST, Devlin et al. 2004) represent a new challenge for bolometric timestreamtomap algorithms. Recent CMB experiments which use detectors similar to those used on BLAST, such as BOOMERANG (Crill et al., 2003) and Archeops (Benoit et al., 2002), only use a handful of separate bolometers. Furthermore, these experiments’ offaxis designs lead to small correlations between detectors. Consequently, the correlations could be ignored at the mapmaking stage, and each detector timestream could be treated as an independent subset of the data. This has changed for BLAST, which has up to 139 detectors per band, with significant correlations induced by the onaxis design, as well as the higher frequencies of the observations. Just by itself, the large number of channels increases the impact of even small time stream correlations, the contribution from which does not integrate down with increasing number of detectors, unlike the uncorrelated noise. The high level of correlation (largely induced by temperature drifts in the obscuring secondary supports in BLAST) make it important that the correlations be handled carefully in the mapmaking process.
This paper describes a Signal And Noise Estimation Procedure Including Correlations (SANEPIC) which has been developed for the analysis of BLAST. This algorithm will also have application to many next generation experiments which will involve both noise correlations between channels (including correlations from the atmosphere) and very large numbers of detectors. This includes the next generation of largerformat arrays for use in ground and balloonbased instruments at microwave and millimeter wavelengths, which will have typically thousands of detectors.
This paper is arranged as follows. We next describe the pertinent aspects of BLAST. In Section 3, the longest section, we set out our basic mapmaking method. Simulations we use for testing the method are described in Section 4, with results presented in Section 5, demonstrating the benefit of accounting for the correlated noise in BLASTlike observations. Finally some of the maps obtained from the June 2005 BLAST flight are presented in Section 6.
2 BLAST observations of the submillimeter sky
The mapmaking procedure presented in this paper has been used to analyze the data from the BLAST 2005 flight. We will use BLAST as a specific example for the application of SANEPIC throughout the paper. We describe BLAST in Section 2.1, we then summarize the preprocessing of the data prior to mapmaking in Section 2.2. In Section 2.3, we derive a model of the data that we will use for the mapmaking.
2.1 BLAST instrument and observations
The Balloonborne LargeAperture Submillimeter Telescope incorporates a 2m primary mirror, and large format bolometer arrays operating at 250, 350, and 500m, designed to have 144, 96 and 48 bolometers, respectively (of which 139, 88 and 43, respectively, were used). The instrument is described in detail in Pascale et al. (2007). The low atmospheric opacity at the operating altitude of km allows BLAST to map the sky very quickly compared to groundbased experiments and to conduct large area shallow surveys as well as very deep surveys of the sky (Devlin et al., 2004). The BLAST wavelengths are near the peak of the spectral energy distribution of cold galactic dust, which gives BLAST the ability to conduct unique extragalactic and Galactic submillimeter surveys with high spatial resolution and sensitivity. BLAST thus enables studies of the distribution of very highredshift galaxies and of star forming regions in our Galaxy.
The typical observing strategy consists of scanning the telescope back and forth in azimuth, covering the entire field by slowly varying the elevation. Crosslinking of the data is assured by scanning the same field at another time of the day. Typical scanning strategies are given in Pascale et al. (2007).
The first scientific flight of BLAST took place in June 2005 from the Esrange Arctic base in Sweden to the Canadian Arctic. A total of 100 hours of data were taken in a variety of Galactic fields. They include a star forming region (Vulpecula) over , described in Chapin et al. (2007), three other fields of similar size in the Galactic Plane (which will be the focus of future papers), an integration towards the ELAISN1 field (see Oliver et al. 2000), the Cas A supernovae remnant over about (Hargrave et al. in preparation), and several compact Galactic and extragalactic sources (Truch et al., 2007). Hereafter we refer to these as the BLAST05 data, to distinguish them from the data taken during the December 2006 Antarctic flight.
2.2 Timeordered data preprocessing
The processing of BLAST data from detector timestreams to the final map product involves several steps prior to mapmaking. Each of these steps is designed to remove a particular (or several) artifact(s) from the data, and sometimes requires iterating, since some effects need to be removed simultaneously. In the following, we summarize the main processing stages leading to the timeordered segments which are used as inputs for the mapmaking process.
We start by identifying events in the data which are sharply localized in time, such as spikes from cosmic rays hits and other spurious sources. We use a method which allows us to discriminate between the different events depending on their signature in the data. Spikes which involve only a single sample are flagged and the corrupted samples are replaced by the average value of the samples in the vicinity. The data are deconvolved from the lowpass filter applied by the readout electronics. This filter has a frequency cut off of approximately Hz and is designed to avoid high frequency noise aliasing. The deconvolution is performed in Fourier space. In addition, we have applied a lowpass cosine filter which limits the noise power from increasing at very high frequency (above Hz) due to deconvolution. We have checked that the noise power spectrum is relatively flat after these deconvolution operations. Finally, cosmic ray hits and other localized artifacts in the data timestreams are detected and cut out. In order to avoid biasing our data products by having systematically more false event detections located where the sky is bright (e.g., when scanning a point source), we iterate this process – we make maps starting from data which has been cleaned using the process described above, subtract the maps from each original dataset, and reprocess the data. The maps calculated at this intermediate stage are obtained by simply rebinning the data into pixels after strongly highpass filtering. The filter applied is a Butterworth filter with a frequency cut off of Hz, which is of the order of the knee frequency of the noise. Even if this operation suppresses most of the intermediate to large scales of the sky signal in the maps, it does not very much alter the signal from localized sources, at least to the level of accuracy needed at this stage, and it has the advantage of removing most of the stripes in the maps due to noise. We have verified that the resulting bias for bright calibrators is less than 1%.
About 2% of the data from the BLAST05 flight was removed due to cosmic ray events. Most of the events affected a single detector timestream, although some events affected a whole array at the same time. Detected spikes (from cosmic rays but also other spurious effects) are flagged over small time intervals of typically 1 second in the data. One second is too large an interval to simply ignore, so the corrupted data need to be replaced with random noise generated in a way that as much as possible preserves the statistical properties of the data. In the later mapmaking stage, which is partly performed in Fourier space, we assume continuity of the data. We could perform this gapfilling with a constrained noise realization (see e.g., Stompor et al. 2002). However, since the gap intervals are significantly smaller than the inverse of the knee frequency of the noise power spectrum (Hz), the noise can be well approximated by the sum of a white component plus a straight line of some slope across the gap. Specifically, we generate white noise in each gap with a standard deviation measured from the data in the vicinity of the gap, and add a baseline with the parameters fit using 20 samples on each side. The white noise generation is for restoring as best as possible the stationarity of the data (generated samples are not reprojected to the map at the end).
After having filled the gaps in the timestreams, we filter out very low frequency drifts which are poorly accounted for in the mapmaking procedure. A fifthorder polynomial is fit to the data and removed from each data segment in order to reduce fluctuations on timescales larger than the length of the considered segment which, depending on the specific case, varies from 30 minutes to a few hours. These fluctuations are poorly described because of the limited number of Fourier modes, and would cause leakage at all timescales (for instance a gradient in the timestream is described by a wide range of Fourier modes), degrading the efficiency of mapmaking. Note that we experimented with various polynomials and other effective highpass filters; we found that the results were not very sensitive to precisely how this is done, but some such filtering is certainly required. The degree of the polynomial was chosen empirically as a compromise between suppressing the artifacts and keeping the large scale signals in the final maps. Using simulations, we have checked that the effect on the transfer function of the signal in the final map is weak. The resulting data segments are then corrected for the timevarying calibration (see Truch et al. 2007) using measurements of an internal calibration lamp (Pascale et al., 2007). Finally, the data segments are apodized at the edges over samples and are highpass filtered at Hz with a Butterworth filter. This filtering has very little effect on the final maps, since modes in the data at lower frequencies that are cut in this way are not expected to contribute significantly to the signals. We discuss the choices of filters further in Section 3.2.
Accurate pointing reconstruction is a complicated procedure for balloonborne telescopes, and this affects the mapmaking task through the pointing matrix (defined in Section 2.3). The pointing reconstruction procedure is described in detail in Pascale et al. (2007). The next important step in reducing the BLAST data, which is calibration of the detectors, is detailed in Truch et al. (2007).
2.3 Model of the data
Having performed the cleaning procedure described in the previous subsection, the resulting data timestreams can be modeled very accurately as the sum of pure signal and pure noise contributions. The data for detector observing at a given wavelength and at time sample can be written as
(1) 
where labels the pixels in the final map, is the pointing matrix for bolometer (whose elements, indexed with time and pixel , give the weight of the contribution of pixel to the sample at time for bolometer ), is the signal amplitude at pixel , and the noise amplitude at time for bolometer is . Summation over repeated indices is assumed here.
Ideally, the element of the pointing matrix is equal to the value of the beam response , where points to the pixel location, is the location of the beam center at time , and is a rotation matrix which depends on the rotation angle at time between the telescope and sky coordinate systems. In principle one could then recover a map of the sky deconvolved with the instrumental Point Spread Function (PSF). However, in practice this would be unacceptably noisy, as well as computationally intractable, because of the prohibitive volume of data. Even although might have mostly zero elements, it is nevertheless a huge matrix. It may be feasible to deconvolve a nontrivial beam response (e.g. like the BLAST05 beams, as shown in Truch et al. 2007) at the same time performing the mapmaking step, through an approximate treatment of the sparse pointing matrix. But we did not pursue that approach here.
In the simple case where the beam is symmetric, the mapmaking problem becomes tractable, provided one restricts oneself to reconstructing a map of the instrumentconvolved sky. We can then consider in Equation 1 as the map of the sky convolved with the beam, and consequently indicates simply where the detector points in the sky at a given time. In this case, is an ultimately sparse matrix with, in the BLAST case, simply a 1 in a single entry of each row. This approach, which has been conventional for CMB mapmaking (although with some adaptation for chopped data), is what we use in the following analysis. It gives no loss of information provided that the map pixels are sufficiently small, and one can simply assign all the flux from a bolometer’s to the map pixel to which it points at each time interval. The requirement for accuracy is that the pixel size is smaller by a factor three or more than the Full Width at Half Maximum (FWHM) of the instrumental PSF; in that case the additional convolution with the pixel shape gives negligible loss of angular resolution.
The noise term in the model represents the sum of all contributions to the timestreams which do not reproject on the sky. This will in general include instrumental noise, fluctuations in atmospheric emission and other loading and unrecognized cosmic ray hits. In general, some of those noise contributions will induce strong correlations between detector timestreams. In this paper, we adopt a very general model of the noise where the noise covariance matrix:
(2) 
(for bolometer indices , and time indices , ) has possibly nonzero elements even for . A key assumption, as we will see later, is that the noise is Gaussian (so that is sufficient to describe all the statistical properties of the noise) and stationary (constraining , see Sections 3.2.1 and 3.4).
In the specific case of BLAST05 observations, a very significant correlation of the noise is found in the timestreams, and we have shown that a more constraining model provides a very good description of the data. An independent component analysis (Delabrouille, Cardoso, & Patanchon, 2003) of the data enabled us to find that the noise and its correlations can be described to a high degree of accuracy by a noise component which is not correlated between detectors, together with a single commonmode component seen by all the detectors at a given wavelength (some correlations is also seen between detectors from different wavelength but we have chosen to treat each wavelength independently). The common part of the noise is instantaneous, meaning that the same commonmode noise is seen at the same time by all the detectors. In our model, the noise term in Equation 1 is then decomposed as:
(3) 
where the first term is the noise which is uncorrelated between detectors and the second term represents the commonmode component of the noise, rescaled by an amplitude parameter which depends on the detector but not on time. This model can be generalized easily to deal with multiple noise components in timestreams.
In the following section, we present a method to reconstruct given the data, in the framework of the linear model (Equation 1) and in the presence of correlated noise between detector timestreams.
3 Mapmaking method
3.1 Maximum Likelihood mapmaking
The use of Maximum Likelihood mapmaking techniques has been developed by many authors for application to Cosmic Microwave Background datasets (Wright et al., 1996; Tegmark, 1997; Borrill, 1999; Prunet et al., 2000; Tegmark et al., 2000; Ferreira & Jaffe, 2000; Doré et al., 2001; Natoli et al., 2001; Prunet et al., 2001; Dupac & Giard, 2002; Stompor et al., 2002; Hinshaw et al., 2003; Yvon & Mayet, 2005). Some other approach are more specific to destriping for Plancklike scanning strategies (Delabrouille, 1998; Maino et al., 2002; Keihänen et al., 2005; de Gasperis et al., 2005; MacíasPérez et al., 2007; Poutanen et al., 2006; Ashdown et al., 2007). Since there is already a large number of publications on the topic, here we present only a very brief overview of the approach of Maximum Likelihood mapmaking techniques.
Assuming the simple linear model given by Equation 1, the loglikelihood of the data can be calculated under the assumption that the noise is Gaussian and stationary. The solution is
(4) 
where is the noise covariance matrix in the time domain, and denotes transpose. Maximizing the above equation with respect to the map parameters (suppressing the pixel indices here for convenience) leads to the following well known estimator:
(5) 
The inverse pixelpixel covariance matrix of the noise in the map is the term in brackets in this equation, i.e.
(6) 
Computation of the solution to Equation 5 is far from trivial for most astronomical applications, due to the large amount of data, and hence this poses a difficult numerical challenge. The noise covariance matrix is a very large matrix of size the number of samples squared, which could easily be millions, while may be more reasonable in size but has no obvious symmetries, and so is still difficult to invert. Nevertheless, we have implemented a method aimed at finding the Maximum Likelihood solution given by Equation 5 when there are a large number of detectors and in the presence of possible correlations in the noise between different detector timestreams.
3.2 Implementation
In the simple case of dealing only with independent noise between detectors, our matrixinversion method is very similar to the MADCAP method, described in Stompor et al. (2002). However, we have developed our new approach to deal efficiently with multidetector data in the presence of correlated noise between detectors (described in detail in Section 3.4). In this section, we summarize the basic ideas for the simpler 1detector 1scan case.
In order to find the Maximum Likelihood solution of the map (Equation 5), we have developed two different algorithms. They both allow us to solve the linear system , with . The first approach explicitly computes the inverse pixelpixel covariance matrix , and we refer to this as the ‘bruteforce algorithm’. The second approach uses iterations which converge to the Maximum Likelihood map without the need for computing , and we refer to this as the ‘iterative algorithm’. Both approaches require as a first step the computation of the inverse of the timetime noise covariance matrix .
Inverse noise covariance matrix
In practice, even when we have knowledge of the statistical properties of the data as described by the power spectrum , the bruteforce inversion of is not tractable because of its enormous size – for a single BLAST detector observing for 10 hours at a datarate of Hz, the matrix has approximately elements. However, if we make the approximation that each data segment is “circulant”, meaning that the beginning and the end of a segment are connected without discontinuity and that there are no gaps in the data, then the matrix is also circulant (see Section 3.6 for a description of how we treat gaps in the data). Circulant matrices are much easier to store and to invert. With this approximation the matrix can be written
(7) 
where the correlation function between samples and depends only on the separation between the two samples. A circulant matrix has the property of having a diagonal matrix counterpart in Fourier space.
Let be the discrete Fourier operator, we have
(8) 
where denotes transpose conjugation, and is a diagonal matrix whose diagonal is described by the power spectrum of the data segment,
(9) 
The inverse of the noise covariance matrix is
(10) 
and because is a diagonal matrix, is also circulant, so that knowledge of only one row is enough to describe the entire matrix. Then the inverse covariance matrix can be written
(11) 
with
(12) 
where represents the inverse Fourier transform. The inverse of the covariance matrix can then be computed directly using the power spectrum of the data. This is a very fast operation (, with the number of samples), and does not require large memory since only one row of the matrix is computed.
The approximation that the data segments are ringshaped or circulant might seem unreasonable, but in the end this only has an effect on a small fraction of the matrix. For data with largescale correlations (data described by a power spectrum, for instance), the approximation implies an assumption that the two edges of the segment are very correlated. In reality, there is obviously little correlation at long timescales compared to short timescales, so extra striping could be introduced in the maps if one steps across the two edges of the data segment. We have addressed this problem in two ways in order to avoid introducing artifacts in the final map. In the case where we explicitly compute the inverse pixelpixel covariance matrix (as in the bruteforce inversion algorithm), for the estimation of entering in the computation of , we constrain for , with being the number of samples in the segment; this is sometimes known as “the MADCAP approximation” in the literature. It is important to note that this approximation cannot be used for the computation of , because it is performed partly in Fourier space (see Section 3.2.3). Instead we have apodized the data at the edges, and we have removed a loworder polynomial (5th order in practice) to reduce fluctuations having typical timescales of order (or larger than) the data segment (see Section 2.2). This is reasonable, since those scales are not well described by a limited number of Fourier modes, and hence will always be hard to reconstruct.
Inverse pixelpixel covariance matrix
The computation of the inverse pixelpixel covariance matrix, which is described in this section, is required only in the bruteforce inversion algorithm, or for an accurate error estimation in the pixel domain.
Since the pointing matrix has only one nonzero element per row in our simple model (one data sample is associated with a single map pixel), the matrix multiplication requires a single loop going across all the nonzero elements of . For most cases, a dominant faction of the mapmaking computing time will be devoted to this operation. If the data are only correlated within a typical length , we have the property: for , and is a banddiagonal matrix (elements separated from the diagonal by more than are negligible). The number of elements to go through in the loop is of the order of, but smaller than , which is hopefully much smaller than the size of the matrix itself.
Unfortunately, if the noise is described by a power spectrum of the form , the correlation length of the noise is basically of the order of the length of the whole dataset. However, the amplitude of the correlation is decreasing for very long timescales and becomes negligible beyond a certain scale. The function can then be artificially set to zero for , with chosen empirically so that the correlation of the noise is low enough for scales longer than , and also that there is very little constraint on the signal at those scales. For the specific case of BLAST observations, we find that s is a good compromise, as illustrated in Figure 1.
The noise strongly dominates the signal in BLAST observations for scales longer than (corresponding to frequencies smaller than Hz), since: (1) this frequency is well below the knee frequency of the noise power spectrum; and (2) there is very little signal at frequencies smaller than Hz, which is more than ten times smaller than the scanning frequency. Therefore we have used s for the computation of the noise covariance matrix. The impact of fixing for in the initial power spectrum is shown in Figure 1. We can see a tight relation between the power spectrum and the inverse covariance matrix, since getting for can be obtained only by modifying the power spectrum at frequencies smaller that 1/.
Computation of
In the computation of , which is necessary for both of our algorithms, the multiplication is performed in Fourier space where the noise covariance matrix is diagonal. We obtain this vector by dividing the Fourier transform of the data by the power spectrum of the noise. Another way to represent this is by considering that since is circulant, is a convolution operation. Assuming that this model holds, the resulting data vector contains whitened noise. The remaining operation just performs the addition of the filtered data sample onto the pixels of (i.e. the map), and hence is very fast.
In the case of BLAST observations, since we are not attempting to recover useful information from timescales larger than s, we perform a highpass prefiltering of the input data at Hz.
Matrix inversion algorithm
In the matrix inversion algorithm, is directly computed, as described in Section 3.2.2. The next step is to solve the linear system . For small maps, in which can be stored in memory, we perform a Cholesky decomposition. For larger maps, we write the matrix to disk and perform an iterative inversion of the system using a conjugate gradient method with preconditioner.
This algorithm allows one to easily perform multiple MonteCarlo simulations, since the pixelpixel covariance matrix, being independent on the data realization, can be computed once (this is assuming that the noise power spectrum is known and does not have to be estimated from the data themselves), and used for all the realizations of simulated data. Another advantage of this approach is that it allows for the exact computation of the errors in the map, which are given by the covariance matrix.
However, the matrix inversion algorithm is generally slower than the iterative algorithm which we will discuss next, and can be used only for relatively small maps; we found that we were limited to around 200,000 pixels if the matrix is written to disk, or less than 20,000 pixels if the matrix is stored in memory.
Iterative algorithm
We now present an iterative algorithm based on conjugate gradient with preconditioner to obtain the Maximum Likelihood solution for the map (a similar algorithm has been used in Ashdown et al. (2007)). Let us rewrite Equation 5, which relates the best estimate of the map with the data, after multiplying both sides by the pixelpixel covariance matrix:
(13) 
If we define as an estimate of the map at iteration , the conjugate gradient method allows to solve the linear system by minimizing iteratively the following criterion:
(14) 
where
(15) 
This criterion is indeed minimum and equal to zero if is the Maximum Likelihood solution.
One can interpret as the weighted variance of the difference between two map vectors. The first of these vectors, , is the inverse pixelpixel covariance matrix times the current estimate of the map, while the second vector, , is a map constructed by simply coadding prewhitened data. We decide that convergence is reached when the quantity (which is much easier to compute than , and also converges to zero) gets smaller than a predefined value. In practice the number of iterations required for convergence is of the order of 100.
The conjugate gradient method is not described here, since it
is a fairly standard numerical tool, and the interested reader can find
many descriptions in the literature.
One advantage of this iterative algorithm is that the computation of the full pixelpixel covariance matrix is not required, and the operation can be done step by step. Indeed, we start by computing , which is an estimate of a “signal” timestream. This operation is equivalent to scanning over the current estimate of the map using the pointing solution. The subsequent operation (which should now be familiar), is carried out in Fourier space, as described in Section 3.2.3 (without applying any extra filtering), and in Section 3.4 for the case of correlated noise between detectors.
This iterative approach is in general much faster than the bruteforce inversion approach, because the most timeconsuming operations are performed in Fourier space. It also requires less memory, since is not explicitly computed. Of course, if there are found to be (or known to be) nontrivial correlations in , then it may have to be calculated explicitly, hence requiring the bruteforce approach. However, provided the pixelpixel correlations only involve relatively few pixels, it should be possible to calculate a restricted part of (or perhaps an approximation for) in a modified iterative approach. A related concept is discussed in Section 3.8.
3.3 Multiscan, multidetector case
In the previous subsection, we presented the general method for the simple case where only a single continuous observation is considered. We now describe how we combine observations from different detectors at the same wavelength, as well as different data segments obtained over different “visits” during the flight, where by “visit” we mean a period in the data which starts after a sufficiently long gap, or after the observation of a different region of the sky.
For convenience, the data vector in our model (Equation 1) now contains all the individual data segments from different detectors, and also within a single channel, concatenated end to end. The noise vector in Equation 1 is defined in a similar manner. The matrix in Equation 1 is then the result of stacking individual pointing matrices. The Maximum Likelihood solution is also written as in Equation 5, with becoming the full covariance matrix of the noise, including all the channels and data chunks. To start with, let us assume that there is no correlation of the noise between data segments. This is a very good assumption if we consider data segments obtained over different visits, but is certainly not a good approximation for segments obtained simultaneously with different channels, since we found a very strong commonmode noise between detectors. We will consider the simple nocorrelation case first, and then in the next subsection (Section 3.4) we will generalize the mapmaking method to account for noise correlations from different detectors.
In the absence of correlations between data segments, the timetime noise covariance matrix is blockdiagonal and each block can be inverted separately. Defining as the subcovariance matrix for the data in segment number , and as the subpointing matrix going from the map to the data segment , the inverse pixelpixel covariance matrix can be written as
(16) 
The computation time for obtaining this matrix is proportional to the number of data segments. In this simple case the computation of can be written
(17) 
where the computation of each term is fast and can be performed partly in Fourier space.
3.4 Detectordetector correlated noise
We now allow for the presence of correlations in the noise between different detectors. In the case of multiple visits, the noise covariance matrix, , still has null crossterms for samples from two different data visits. Therefore, if the data vector is sorted by visit, then is block diagonal and each block contains the correlation coefficients between all the detectors for the samples within the time interval defined as a single visit. Each visit can be treated independently, since the submatrices can be inverted separately, and Equation 16 is still valid, but in this case is the label for the blocks in . In the following, we therefore focus on a single visit, and consider observations by all the detectors; the generalization to multiple visits should be clear.
To simplify the notation, let denote the noise covariance matrix for the visit being considered, with and being the data and noise vectors, respectively, containing the timestream segments for all the detectors put end to end, and being a block of of size , corresponding to the noise correlations between detectors and . Let us define the multichannel Fourier transform operator such that
(18) 
with containing endtoend Fourier transforms of each data segment. is a blockdiagonal matrix, and each block is the Fourier transform operator for one data segment.
In Fourier space, the noise covariance matrix can be written
(19) 
If we consider a single block of the noise covariance matrix for detectors and , we obtain:
(20) 
Under the assumption that the data are stationary and continuous at the edges (see Section 3.2 for a discussion), is a circulant matrix, since each element depends only on the time interval . is then a diagonal matrix with the diagonal given by the crosspower spectrum of the noise between detectors and :
(21)  
Here is the noise covariance matrix of size for a given mode , where is the total number of detectors. The computation of the inverse of is straightforward, since each Fourier mode can be treated independently. If is the inverse noise covariance matrix for mode , the same relation as in Equation 21 applies between and all .
From Equation 19, we can calculate the inverse covariance matrix of the noise in real space:
(22) 
Then, a block of between detectors and can be written
(23) 
Because is a diagonal matrix (as discussed previously) is circulant, and is related to the inverse of the matrix containing the cross and autopower spectra of the noise:
(24) 
From this relation, we can see that in practice is relatively easy to construct, since each of its blocks (referring to each pair of detectors) is a circulant matrix, so only a row of each submatrix needs to be calculated using the Fast Fourier Transform. Finally, in the case the noise covariance matrix is used for multiplication in real space (as in the bruteforce algorithm), the same approximation described in Section 3.2 is performed on each block of , i.e., for (or for if ). The global structure of the final inverse noise covariance matrix is illustrated in Figure 2.
In our model of the data (Equation 3) in which all the correlations between detectors are described by a single commonmode, can be related to the power spectra of the commonmode and of the uncorrelated part of the noise between detectors:
(25) 
This relation can be generalized if multiple commonmode components are present: would then be a mixing matrix and becomes the covariance matrix of these noise components.
Having computed the inverse noise covariance matrix of the timestreams, we can express (using Equation 6) the noise covariance matrix in the pixel domain:
(26) 
where, as before, and label the detectors. The computation time for is now proportional to the number of detectors squared. The calculation of is straightforward:
(27) 
and this is fast, since (as already shown) the operation is a convolution, which can be performed in Fourier space. One can see from Equation 24 that can be expressed directly with respect to the cross and autopower spectra of the noise:
(28) 
The formalism presented above can also be generalized easily to deal with detectors operating at different wavelenths. The map vector could be merging different maps at different wavelengths and the noise matrix would account for all the correlations of the noise between detectors. The joint multiband mapmaking would be suitable in practice when some contaminations from thermal fluctuations in the instrument or atmospheric emission are present, because they correlate the noise at all wavelengths.
3.5 Noise power spectrum
The Maximum Likelihood solution for the final map depends on the noise power spectra for each dataset (through in Equation 5), which are assumed to be perfectly known. However, in practice the noise power spectra have to be inferred from the data themselves and some uncertainties are associated with this iterative process.
In practice a first (approximate) estimate of the noise power spectrum can be obtained by rebinning the power spectra of each data segment, neglecting the contribution of the astrophysical signal in the timelines. Indeed, for most of the fields observed with BLAST, and in particular for blank extragalactic fields (like the ELAISN1 field in BLAST’s flight from Sweden), the noise is highly dominant over the sky signal at all frequencies. However, this is not true for measurements of bright regions in the Galactic Plane. Therefore, for the first iteration’s noise estimate we focus on the data taken for about 6 hours while scanning the deepest extragalactic field (ELAISN1) and use this to estimate the noise power spectra, which then become the noise input for making the first set of maps of all of our fields. This approach is not entirely satisfactory, since the noise is not stationary during the flight – the noise power spectra are seen to vary over long timescales, although they are quite constant within a single visit of each field. This nonstationarity appears most obvious when there are variations of the scanning strategy between different visits (a change of the scanning frequency induces a variation of the location of peaks in the power spectrum), but can also occur due to variations of detector loading or the detector bias being changed during the flight.
Because of the observed nonstationarity of the noise, we would like an independent estimate of the noise power spectra for each visit for each of the observed fields. After starting from a first estimate of the noise power spectra based on ELAISN1 data as described above, we adopt an iterative approach between maps and noise power spectra. At each iteration, the estimated maps are subtracted from the data, prior to noise power spectrum estimation. We can summarize our procedure in the following steps:

Estimate from using a field known to have little signal;

Estimate from ;

Reestimate from Equation 5 using , and iterate on these last two steps until convergence is achieved.
We stop iterating when the noise power spectra do not vary by more than 1 per thousand from iteration to iteration (and find that in practice only 3 to 6 iterations are necessary to reach convergence).
We now focus on how we estimate the noise power spectra in steps 1 and 3. For the simplest case, where no correlations are assumed between detectors, we simply compute a binaveraged power spectrum for each data segment:
(29) 
where, for bin number , , labels the data segment, and is the data vector from which an estimate of the map has already been subtracted. We have chosen logarithmic spacing between bins and an estimate of for each mode is obtained by logarithmic interpolation, which leads to a smooth power spectrum estimate.
In the more complicated case where correlations between detectors are assumed to be important and therefore are not neglected, we must estimate, for every iteration, each cross and autopower spectrum of the data between detectors, , which enter into the computation of the inverse noise covariance matrix (Equation 24). Each crosspower spectrum could be directly estimated as in Equation 29 (using and in the formulae for detectors and ), but instead we choose to reduce the number of parameters to estimate at each step, by assuming that the data are described by a common mode between detectors plus independent noise (Equation 3). In the framework of this model, the expected cross and autopower spectra depend directly on the following parameters: , the amplitude of the common mode in each channel; , the power spectrum of the commonmode part; and , the power spectrum of the noise component, which is independent between detectors. The relation between the model of and the parameters has been shown in Equation 25. These parameters are typically not known a priori, and must be measured using the data themselves. We use a blind ‘component separation’ method developed for an entirely different problem in Delabrouille, Cardoso, & Patanchon (2003). This allows us to obtain a single estimate of all the parameters described previously, by simultaneously using all the observed timestreams of a given field for all the detectors in a specified channel (i.e. at a single frequency). The method is known to be the Maximum Likelihood solution for a Gaussian and stationary model of both the noise and the commonmode. The cross and autopower spectra are then computed following Equation 25, using these same estimated parameters. Figure 3 shows the estimated noise power spectra in a sample of BLAST data for one representative detector using three hours of timestreams during scans of the ELAISN1 field (which is known to be essentially devoid of signal).
The autopower spectrum is shown, as well as its decomposition in terms of the commonmode power spectrum and the uncorrelated noise power spectrum.
Because our estimate of the converged map of the sky is not perfect and contains contributions from residual noise, then in subtracting a simulated signal timeline from the data to estimate the noise power spectrum we reintroduce some noise to the data which could potentially bias our estimate of the noise power spectrum (Ferreira & Jaffe, 2000; Hamilton, 2003). However this effect is greatly reduced by the large redundancy in each pixel of the final maps, as a result of the many repeated scans and the large number of detectors at each wavelength. The bias can be neglected to first order for BLAST, since the noise level per pixel in the final map is much smaller than the noise in individual detector timestreams.
3.6 Dealing with gaps in the data
In order to derive the formalism presented so far, we have assumed that each data segment is stationary and hence consists of a continuous series of data points. However, we have seen in Section 2.2 that the BLAST05 data contain multiple gaps of typical size less than one second. The amount of data in these gaps is a few percent of the total. In order to reasonably restore the continuity of the data we have filled the gaps with random noise, as described in Section 2.2. The data samples generated in the gaps are not reprojected into the final map but are directed to “dummy” pixels. In principal, the optimal approach would be to create one dummy pixel per flagged data sample, avoiding the possibility of several simulated samples falling on the same pixel (through the rebinning of flagged pixels we do not want to introduce any spurious constraints for the mapmaking process which could arise from adding crossings over different time intervals). However, the approach of using one dummy pixel per gap sample is impractical, because the total number of dummy pixels would be excessive for typical BLAST observations, and the pixelpixel covariance matrix (which has size the total number of pixels squared) would be prohibitive to store and compute.
We have adopted a simpler approach, which does not lead to the mathematically exact solution, but comes very close (as has been shown in simulations). This consists of rejecting (for the computation of ) all the elements of associated with flagged samples. This is equivalent to removing from the rows and columns corresponding to the dummy pixels before the inversion of the matrix, as opposed to after inversion, which would be the correct treatment discussed above. This approach is also equivalent to assuming null offdiagonal terms for those rows and columns. However, such dummy pixels are obviously correlated to some degree with the real pixels in the map, and hence this cannot be entirely correct.
Nevertheless, we have verified that this approximation has a small impact (a few percent only) on the final map for most of the angular scales which are sampled, although we found some differences at very large angular scales, sizes of the order of the map size. For now, we have not put much effort into recovering those very large scales because they are subject to other effects, as discussed in Section 5.2.2. The minimal impact on the final map has been verified using pure signal simulations and by comparing the results obtained between our simple approach and the correct mapmaking solution. This approach works to a high degree of accuracy on most scales because the gaps are small and do not introduce important discontinuities in the timestreams.
We have used this simple approach because it gives sufficiently accurate results over the relevant angular scales, while being simple and fast to implement. However, another iterative procedure could be adopted, which would lead to the exact solution. In this approach, we define two maps. The first map “A”is made from only the uncorrupted (i.e. real sky) samples, while the second map “B” is obtained from projecting the simulated (i.e. for gapfilling) samples. The difficulty arises in deciding what to do when simulated data from different scans or detectors fall in the same pixel – one might want the “generated signal” (and not necessarily the noise) to be identical in both measurements, in order to satisfy the mapmaking hypothesis. If this condition is not satisfied, then some artifacts may be introduced into both maps. Here is a solution to this problem:

generate a first set of maps “A” and “B” after filling the gaps in the timestreams with white noise + a linear baseline.

fill the gaps in the data with the best estimate of the signal in map “A” together with white noise + a baseline which is fitted in the gap vicinity of the data “minus” signal timestream.

recompute the maps “A” and “B”. Step 2 ensures that the signal is the same for each generated sample falling in the same pixel of map “B”.
This approach can also be coupled with the procedure for estimating the noise power spectra described in Section 3.5. Preliminary results indicate that this approach works in practice. Detailed studies will be presented in a future publication.
3.7 Pixel constraints
For some specific observed fields we may have strong priors about the sky emission at a given location. For instance, we know that over some regions the astronomical signal should vary very smoothly or should be very weak with respect to the noise, at least outside some localized region. This is the case in particular when we map bright extragalactic sources in order to calibrate the detectors and estimate the beams; in these cases, regions beyond some predefined distance from the beam center can be assumed to have null flux (or a constant relative flux in the map, since we do not have access to the DC level in maps). If we really have strong prior knowledge that we are dealing with a bright localized region, then we can take a further drastic step – we can constrain the map to have the same value in some domains of the sky by defining a single pixel containing all the data samples falling in that region.
In practice, we define a small box centered at the source location and constrain the part of the map outside this box to have a constant value. This is a very efficient way to remove stripes from the map, since the extremities of all the paths across the map are readjusted. We have used this technique to make maps of the isolated calibrators observed by BLAST (Truch et al. 2007).
3.8 Error estimation
The variance of the noise in each pixel of the final map and its correlations are directly given by the pixelpixel covariance matrix . This is true given the following assumptions: that our model of the data holds, in particular that the noise is a purely Gaussian random process, which may not be the case in practice at low frequencies; and that our estimate of the samplesample noise covariance matrix is accurate enough that the errors do not propagate significantly into the final map. As already mentioned, we never explicitly compute the covariance matrix, but rather its inverse. The direct inversion would take a prohibitive computation time for most applications. However, to first order, we can obtain an estimate of the errors by inverting the diagonal part of only, neglecting the offdiagonal terms. This is equivalent to assuming that the noise in the final map is white. We have checked with the help of simulations that this very simple approximation is accurate to better than 10 percent for all our BLAST05 fields, even for those with very poor crosslinking.
Provided that the size of the map is reasonably small, so that we are able to explicitly calculate , we can obtain an accurate estimate of the errors for a limited (and small) number of pixels in the map. The variance for pixel , and the covariance with respect to the other pixels of the map, can be computed by solving the linear system:
(30) 
where is a unitary vector with a single 1 for pixel . If the Cholesky decomposition of the matrix has already been performed for the mapmaking procedure, the computation of Equation 30 is relatively fast and hence can be carried out for a grid of nonadjacent pixels, for example. This can be used to check the validity of the error prediction approximation described in the previous paragraph.
3.9 Computational requirements
For the bruteforce inversion algorithm (in which the full inverse pixelpixel covariance matrix is computed), five minutes of computation with a single GHz processor are needed to process 2 hours of data from a single detector at a rate of Hz. The computational time is proportional to the number of samples if this is longer than the assumed correlation length of the noise in the data (which has been evaluated to be s in BLAST05 timestreams). If noise correlations between detectors are also to be accounted for, the computational time is proportional to the square of the number of detectors.
Most of the computing time is spent on calculating . Inversion of the linear system to estimate the map is relatively fast (a few minutes to a few hours for maps of several square degrees in size).
For the iterative algorithm, about two minutes of computational time is required for two hours of data (under the same conditions described above). This assumes that 100 iterations are necessary to reach convergence (which is a realistic number for most applications), and the algorithm scales with . The situation is much better than for the brutforce inversion algorithm if correlations between detectors are included. In that case, the algorithm scales with the square of the number of detectors if this exceeds about 40. If there are fewer than about 40 detectors then the algorithm scales linearly with the number of detectors. As an example, if there are 100 detectors, then including noise correlations between the detectors increases the computational time by a factor of four with respect to the “nocorrelation” case. The full processing of 10 hours of BLAST05 data at all wavelengths, including detectordetector correlations, can be done with a single processor in a few days.
In terms of memory, the bruteforce inversion algorithm requires storage of the full matrix. However, for the iterative algorithm, only vectors of the size of the maps need to be kept in memory, which is much less demanding.
4 Simulations for testing SANEPIC
We now focus on the application of SANEPIC to data. Our aim is to develop tests to validate our method using simulated BLAST observations. We derive conclusions about how well low frequency noise in the maps can be reduced, depending on observational parameters such as scanning strategy, and we compare the results obtained with those from simpler methods based on filtering the data, e.g., commonmode subtraction. In this section, we describe the simulations performed to test the SANEPIC method.
We have generated several different sets of simulations of BLAST timestreams. Each set of simulated timestreams, representing one particular observed field, is generated for all the BLAST detectors used for the analysis of real data (132 at 250m, 78 at 350m and 39 at 500m) and is the sum of simulated astrophysical signal, independent noise, and commonmode noise between all detectors. The noise is generated randomly with Gaussian statistics, given fixed power spectra derived from real BLAST05 data. Figure 3 shows an example of the power spectrum of the noise in the data used as input to the simulations for one of the fields. The part of the noise which is independent between detectors is generated for every detector timestream and has a power spectrum well described on average by a relatively flat plateau for frequencies larger than about Hz, and by a part proportional to for smaller frequencies (these characteristics vary slightly from detector to detector). Our knowledge of the real bolometer noise power spectrum at low frequency is limited by the very dominant common mode. In simulations, the commonmode noise is generated once for all detectors and has a power spectrum very well fit by a power law with an index equal to approximately 2.5, together with some broad peaks, the largest being at the scanning frequency (the amplitude of the peak depends on the scanning strategy and the observed field). The commonmode power spectrum has an amplitude such that it reaches the level of the independent noise at about Hz (see Figure 3). The generated commonmode timestream is multiplied by an amplitude factor which varies from detector to detector by and is added to the simulated detector timestreams. The amplitude factors used for the simulations have been estimated from the data themselves.
In order to represent the astrophysical signal, we have simulated simple maps of diffuse emission with a power spectrum proportional to , as for typical Galactic cirrus emission (e.g. MivilleDeschênes et al. 2007). Maps are generated following Gaussian statistics with a resolution of 1, much higher than the typical pixel sizes in the final maps, in order to reduce artifacts due to repixelization. The amplitude of the fluctuations of the simulated map is chosen to match the expected level of signal in each observed field. The simulated maps are scanned using BLAST05 pointing, and pure signal timestreams are generated for each detector. Signal and noise timestreams are added at the end of the procedure (but see Section 5 for an explanation of why this operation is not always carried out).
We have generated two sets of simulations which correspond to two different fields observed by BLAST. We selected two fields that were observed with very different scanning strategies, since the performance of the mapmaking procedure is very dependent on scanning strategy; this allows us to test SANEPIC in two very different configurations. In the first case the scanning was performed mainly in a single direction over a short time interval, while in the second case the field was observed several times during the flight at different scanning angles, to achieve significant crosslinking in the map.
The first dataset uses observations of the Cas A supernova remnant emission which comprises about 20 minutes of data. BLAST observations of this field and derived conclusions will be described in detail in Hargrave et al. (in preparation). The rectangular region mapped has a size of the order of and was scanned two times back and forth over a short time interval. We have generated simulations corresponding to all the detectors at 250m (we used a total of 132 detectors).
The second dataset reproduces the observations of the intermediate velocity cloud IVC G86.5+59.6 (hereafter ‘G86’). Simulations include four different visits of the field performed during the flight at very different time intervals (ranging from a few hours to more than a day). Each continuous observation segment has a size which varies from one to two hours. Two scanning directions are dominant, which form an angle close to 50. The region covered has a size of about on the sky. Simulations for this field are performed specifically for all the 500m detectors (41 detectors used). Similar Monte Carlo simulations at 250m would have taken a factor of 10 longer, while we believe that the conclusion would remain unchanged.
A total of 20 sets of simulations of the observations for each field have been performed. For each set we vary the realization of the noise and of the signal input map. About four hours of computing time are needed to create one realization of a full set of simulations of G86 with a single processor, compared to a few minutes for the Cas A simulations. This is using the precomputed full pixelpixel covariance matrix, which was also used to analyze the real data.
5 Results from simulations
We now present the results obtained with SANEPIC applied to the two sets of simulated data. In each case, we compare the final map with other maps obtained using simpler mapmaking procedures. For these tests we have assumed that the noise power spectra are perfectly known, rather than estimated separately from each dataset; in practice we fix the noise power spectra to be the ones from the simulations. We have verified that relaxing this constraint has almost negligible change on the final maps.
For these simulated datasets, we have applied some preprocessing of the timestreams before applying SANEPIC, just as we do for the real data. We have systematically removed a 5th order polynomial from each timestream segment and weakly highpass filtered the data, as described in Section 2.2. Finally, for the gapfilling we flag the simulated data at the same locations as in the real data in order to check the influence of the flagging procedure in the final maps.
In the following we have applied SANEPIC independently to pure noise timestreams (containing independent noise and commonmode noise, but without simulated astrophysical signal) and to pure signal timestreams. This procedure allows us to easily derive conclusions about the noise properties in the final maps, as well as about the signal, without biasing the results, because SANEPIC is a linear method (as shown by Equation 5). This is only strictly true if the noise power spectrum is fixed as done here, and not estimated simultaneously along with the maps. Then, applying SANEPIC on pure noise timestreams and pure signal timestreams independently and adding the two final maps is rigorously equivalent to applying SANEPIC on signal plus noise timestreams. This has been checked numerically, and we find that the difference is consistent with double floating precision error. An important consequence of this is that the properties of the noise in the final map are independent of the signaltonoise ratio.
5.1 Case without crosslinking
Noiseonly timestreams
We first study the maps resulting from the noiseonly timestreams in the configuration of Cas A observations. The chosen pixel size of the map is 25 and matches the pixel size of the maps discussed in Hargrave et al. (in preparation). We compare the noise maps obtained from three different procedures:

Case 1: use SANEPIC with the correct treatment of the correlated noise.

Case 2: use SANEPIC fixing the correlation of noise between detectors to zero and fixing the noise power spectrum for each detector to the power spectrum of the sum of uncorrelated noise and common mode. This procedure is very similar to more standard mapmakers in the literature (e.g., Stompor et al. 2002).

Case 3: make a simple reprojection of the data onto a pixelized map by simply averaging the data falling in each pixel, after having filtered the timestream data with the same very weak lowpass filter used for SANEPIC. This procedure is sometimes called “coaddition”.
Figure 4 shows computed noise maps for one of the realizations of the noise in each of the three cases.
As expected, the map obtained with the simple pixel binning approach contains a very large amount of low frequency noise, with strong striping visible along the scan direction. Residual low frequency noise can also be seen in the map obtained using SANEPIC without accounting for the noise correlations between detectors. We do not expect this method to be very efficient, since it is very nonoptimal in cases (such as this example) where a very large fraction of the noise is correlated between detectors. In contrast, the noise map obtained with SANEPIC is quite satisfactory, showing reduced power at low frequency as compared to the previous case. Nevertheless, some very weak excess power is seen in the crossscan direction. This is expected, since the map is not crosslinked, and very poor constraints can be put on the crossscan directions at low spatial frequencies (two positions in the map separated by more than the size of the array in the crossscan direction are observed far apart in time).
In order to quantify the level of low frequency noise in the maps, we compute the 1D power spectra of the maps, averaged over the 20 realizations of the simulated data. For the computation of power spectra, we take into account only the central part of each map, where the level of redundancy in the observations is high (we use only the highest signaltonoise region in the maps). To do so, we apply an apodized mask to the maps going smoothly from 0 at the edges to 1. Figure 5 shows the noise power spectra in the three cases.
The noise level in the simple reprojection map is obviously very poor at all scales. Both of the other mapmakers reach the white noise level for scales smaller than 3 and have excess power at larger angular scales. Nevertheless, the gain between full SANEPIC and SANEPIC without correlations is very important at all scales larger than about 2 and reaches a maximum value of about 10 at around 20 angular scales. An interesting fact is that the knee frequency of the noise power spectrum in the optimal case here corresponds to the inverse of the physical scale of the detector array in the crossscan direction (which is of the order of 6). Indeed, there are no observational redundancies on scales larger than the array in the crossscan direction in the absence of crosslinking in the map. Thus the very long timescale noise present in the timestreams is not efficiently removed and reprojects in the final map at large angular scales. This effect is also present along the scan direction, but with a lower amplitude as the map is scanned back and forth. The trend of the large angular scale power spectrum of the noise in the map just follows the trend of the low frequency noise power spectrum in the timestreams. We will see in Section 5.2 that this effect is reduced when there are multiple scanning directions in the map.
In order to determine the direction in which the noise power is strongest in the map, we have also computed the 2dimensional noise power spectrum. The map of the 2D power spectrum of the noise obtained with SANEPIC (noise correlations included) is shown in Figure 6.
The large bright spot around the center corresponds to a relatively isotropic component of correlated noise (at least at large angular scales). It contains a large fraction of the noise power at large angular scales (seen in the 1D power spectrum in Figure 5). A smaller, but significant fraction of the correlated noise is concentrated in directions perpendicular to the scan direction, as can be seen in the figure. As already discussed, the reason for this excess power is that the noise in the crossscan direction is poorly constrained. This crossscan component of the noise is significant all the way up to the pixel scale.
Signalonly timestreams
We now focus on the signalonly timestream simulations. In order to demonstrate the superior performance of SANEPIC relative to simpler methods based on data filtering, we compare with a mapmaking method which consists of the following: we first remove the whole array average from each detector timestream and then make maps using SANEPIC, assuming no correlations between detectors. Removing the array average reduces the signal to almost zero for scales larger than the array and so we expect no large scale structures to survive in the map. This SANEPIC “commonmode subtraction” method is still a better procedure than just reprojecting the data (after commonmode subtraction) with a well chosen filtering to suppress noise drifts (at Hz, for instance, since that corresponds to the knee frequency of the independent part of the noise). This latter method is commonly used in the submillimeter community and is referred to as “sky removal” in reduction of SCUBA data (Jenness et al., 1998).
Figure 7 shows the input map for one of the signal realizations (left panel), as well as the maps obtained with SANEPIC (correlations included, central panel) and with the commonmode subtraction method (right panel). Results are expected to be worse in the second case, because of the extra filtering and also because SANEPIC gives less weighting to modes at lower frequency which are more contaminated by independent noise.
We can see from Figure 7 that part of the very large scale fluctuations with sizes of the order of the map are removed using SANEPIC, but apart from those very large scales, the input map and the SANEPIC map look very similar. More differences can be seen in the map obtained with the commonmode subtraction method. This is quantified in Figure 8, which compares the 1D power spectra of the two output maps, averaged over 20 simulations and multiplied by . Recall that the input spectrum varies as and so deviations from a flat line are the result of the mapmaking reconstruction. Note that the vertical scale is linear in Figure 8.
With SANEPIC, the power of the reconstructed map decreases for scales larger than about 30. There are three reasons for this: the power spectrum is computed over only a small fraction of the sky (and for an apodized map), so that structures of the order of the map size are never fully recovered; there is weak filtering of the timestreams at Hz and through the 5th order polynomial subtraction; modes in the maps that are very weakly constrained in the mapmaking procedure tend not to be reconstructed through the matrix inversion procedure, since the matrix is very illconditioned and numerical problems occur. The last two effects are the dominant ones. As a result, modes which are preferentially filtered are those which lie perpendicular to the scan direction.
At angular scales smaller than about 2, the power slightly decreases due to the smoothing effect of the pixelization. The transfer function at those scales is well described by a sinc function.
Turning now to the commonmode subtraction method, the power in the map is significantly reduced for scales larger than about 10, and drops rapidly to zero. This is because the commonmode subtraction removes power on all scales larger than the array. On smaller scales, the filtering effect is relatively weak, and is reduced when the number of detectors increases.
For these particular simulations the commonmode subtraction method (using SANEPIC, but with no correlations) does not in fact perform very poorly compared to the SANEPIC optimal approach. This is because the observed field is small, with a size just a few times bigger than the array, and structure at scales smaller than the array size are not strongly affected. This particular map is also not crosslinked. However, the situation is different for large crosslinked maps like the Vulpecula field, as discussed in Section 6.2.
5.2 Case with crosslinking
Noiseonly timestreams
We now focus on the set of simulations of the G86 field at 500m. As in the previous example, we first examine the maps resulting from noiseonly simulated timestreams using three methods: optimal SANEPIC (with noise correlations taken into account); SANEPIC without considering noise correlations between detectors; and the simple coadd method. The chosen pixel size for the reconstructed maps is 1, which allows for inversion of the covariance matrix with a single processor and hence rapid Monte Carlo simulations. The conclusions drawn would remain unchanged if the pixel size was reduced.
Figure 9 shows the final noise maps in the three cases for one realization of the simulations.
The simple reprojection map is obviously very stripy and would be of little use as an estimate of the signal; nevertheless it helps to visualize the directions of scanning in the map. We can see two main directions covering the central region of the map, oriented at about 50 to each other. A third scanning direction is also visible but has a much smaller weight. The central part of the map is the crosslinked region, where we expect the more optimal mapmaking procedures to excel.
In the map obtained using SANEPIC without considering noise correlations between detectors (middle panel of Figure 9), some large scale noise is still visible in the map, even if almost no residual striping is apparent in the crosslinked region. Indeed, too much weight is given to the large timescales in the timestreams (which are basically common between all the detectors) as compared to the smaller timescales for which there are more independent measurements, because the degree of correlation of the noise is weaker. The residual noise at large angular scales is much weaker in the SANEPIC map in which we account for the proper correlations of the noise between detectors (left panel of Figure 9). The noise in this map looks particularly white.
The noise level in each simulated map is quantified in Figure 10, showing the 1D power spectrum of residual noise averaged over 20 simulations, and computed over the crosslinked region (which has a diameter of about 100).
While several orders of magnitude are gained in the noise power at all scales using the SANEPIC “no correlation” method as compared to the simple reprojection method, accounting for the correlations with SANEPIC allows us to further reduce the noise power on scales ranging from 20 to the size of the map by an additional factor of . Toward smaller scales, the gain between the SANEPIC correlation versus no correlation test cases is still very significant down to about 10, where the both methods start to approach the white noise level. In the optimal map obtained with SANEPIC the ratio between the noise power spectrum at large scales and the white noise level is around 20, which is relatively small. Figure 11 shows the 2D power spectra of the noise in the SANEPIC map averaged over 20 realizations. As expected, the large scale noise is more important in directions perpendicular to the scans.
Signalonly timestreams
We now focus on signalonly simulations for G86. As in the case of the Cas A configuration, we compare the performance of SANEPIC with respect to the simpler commonmode subtraction method. Figure 12 shows the pure signal input map for one realization of the simulations (left panel) compared with two recovered maps. The first output map is obtained with SANEPIC including correlations between detectors (middle panel), while the second is obtained by subtracting the common mode between all detectors, followed by applying SANEPIC, but neglecting noise correlations between detectors (right panel). We see that the very largest scales are not recovered by SANEPIC. This is because of the weak filtering applied to the timestreams and, to a greater extent, because of inversion problems with SANEPIC on scales of the order of the size of the map (these scales are very poorly constrained by the mapmaking procedure). In the commonmode subtraction method, only the very largest scales are suppressed by the mapmaking procedure itself, but the filtering effect is more dramatic and extends to much smaller scales.
The effective filtering is quantified in Figure 13, which shows the power spectra of output maps averaged over 20 simulations of pure signal timestreams. Again the power spectra are multiplied by for comparison purposes.
In this case SANEPIC works well on scales up to about half a degree, above which it fails to recover structure in the map; this limit corresponds to scales of about a quarter of the map and larger. The filtering effect is much more pronounced in the map produced with the commonmode subtraction method, being strong for all scales above around 14. This is of course expected, since subtracting the average of the array strongly reduces the signal on scales larger than the array size. Therefore, in order to recover large and intermediate scale structure in the maps, it is beneficial to use SANEPIC instead of other methods that are based on simply filtering the data.
5.3 Advantages of crosslinking
The relative level of residual noise at low spatial frequency in the maps is significantly reduced in the G86 observational configuration as compared to the Cas A configuration. The fundamental difference is that the G86 observations contain multiple (essentially two for most of the data) scanning directions, while the Cas A observations are realized with only two passes across the field in the same direction. Multiple scanning directions give a huge number of additional constraints for the mapmaking procedure. In particular, large scale structures in the map are much better recovered in directions parallel to scans, because the noise there is smaller. Thus having multiple scanning angles allows for recovery of the sky fluctuations for all directions, and ends up giving almost no weight to the individual loosely constrained crossscan kmodes.
Differences in the results for maps from these two example scanning strategies can be quantified in two ways. First of all, for the G86 scanning strategy, the transition between white noise and “excess” large scale noise in the map occurs at a scale around 10, while the same transition occurs at a scale of around 3 for the Cas A scanning strategy. Secondly, the ratio between large scale noise power and white noise power is larger by more than two orders of magnitude for Cas A than for G86. On the other hand, some caution should be taken to not overinterpret this comparison, because the pixel size we used is larger for the G86 map (1) as compared to the Cas A map (25), and therefore the number of crossings per pixel is greater for the G86 map. Nevertheless, this simulation exercise has demonstrated that crosslinking in the map is extremely beneficial, particularly for recovering the large scale structures in the map.
5.4 Mapmaking transfer function
When carrying out a complex data processing procedure, it is important to check whether the results are biased in any way. We have found that the transfer function of the mapmaking procedure, defined as the ratio between the amplitude of fluctuations in the output pure signal map relative to the input map, is not always exactly unity, even for intermediate and small angular scales in the map. For example, in the Cas A configuration at 250m the fluctuation amplitudes in the final map are reduced by 3% on average as compared to the input map, almost uniformly across spatial frequencies and directions. This global discrepancy reaches the level of 9% at 500m. Moreover, it is also present for the G86 configuration. We believe that this reduction is due to the fact that the pixelpixel covariance is illconditioned and numerical imprecision occurs in the matrix inversion. We find that the bias tends to be smaller when the number of detectors is larger and also when the number of constraints increases, like when we have multiple scanning directions, or when we map isolated bright sources (presented in Truch et al. 2007) and constrain the data outside a defined region to have a constant flux (see Section 3.7 for details of this procedure). Since this bias can be estimated using simulations, it is straightforward to correct for. We have found that it is not always important, e.g. for the large Galactic map in the Vulpecula region analyzed in Chapin et al. (2007) the bias is negligible.
6 Application to real BLAST05 data
Now that we have looked at signalonly and noiseonly simulations, we now turn to real BLAST data. In this section, we show maps of two example fields from the BLAST05 data which have been obtained using SANEPIC.
6.1 Cassiopeia A
Figure 14 shows the map obtained from the observations of the Cas A field at 250m using SANEPIC including full consideration of the noise correlated between detectors. Detailed analysis of the maps at the three wavelengths is described in Hargrave et al. (in preparation). The properties of the noise and the transfer function of the signal in the map have been studied in detail from Monte Carlo simulations in Section 5.1. Results from such simulations are used to characterize the map. The power spectrum of this map has been compared to results from simulations in Figure 5.
The map structures are relatively smooth, due to the BLAST05 point spread function, which has a width of the order of 3, causing the drop in the 1D power spectrum below those spatial scales. The signal clearly dominates over noise on angular scales larger than about 3, and the diffuse structure should be reliable up to a large fraction of the overall map size.
6.2 Vulpecula region
Another field observed during the BLAST05 campaign is centered in the Galactic Plane close to the open cluster NGC 6823 in the constellation of Vulpecula. The region mapped has a size of about and was chosen for its highmass star formation activity. Complete analysis of this observed field is presented in Chapin et al. (2007).
A few hours of these data were taken at different time intervals during the flight. By design, this field has been observed with very different scanning directions, and is therefore it should be possible to recover diffuse large scale structures.
The map of the observed region at 250m obtained with SANEPIC is shown in Figure 15 (left panel). For comparison (right panel), we have computed another map using a much simpler method which consists of removing the array average signal at each timestep for all of the timestreams and reprojecting the data onto the map after filtering. This is like the “sky removal” procedure often carried out for groundbased submillimeter data (see also Section 5.1.2). One can see that it suppresses almost all the diffuse structure in the map.
No residual striping is visible in the map obtained with SANEPIC (left panel of Figure 15), mainly due to the presence of multiple scanning directions for this field. Large scale structures in the map are successfully recovered with SANEPIC, as can easily be seen by comparing with the right panel of Figure 15. This recovery applies to scales which are significantly larger than the array size. However, the resulting effective filtering after applying the “array average subtraction” method induces negative signals near bright sources in the map, while no such filtering effect is seen in the SANEPIC map (except perhaps near the edges of the map). This shows that optimal mapmaking methods (in the sense of least squares) like SANEPIC are better suited to recover point sources in the maps as well as diffuse structures.
7 Conclusions
Large format detector arrays operating at farIR and submillimeter wavelengths are becoming the norm, rather than the exception. Groundbased instruments are plagued with commonmode emission arising from the Earth’s atmosphere. And as we have found with BLAST, the same applies to high frequency balloonborne instruments, where we see correlated noise from thermal as well as atmospheric effects. There is an expectation that even upcoming satellites might be faced with similar issues, because of thermal variations in the spacecraft, for example. In general, we expect that correlated noise between detectors will be a major issue which all such experiments have to deal with, and we expect that the SANEPIC approach, which we have described here, will be widely applicable. Indeed, there is evidence from existing arrays (e.g., SHARCII and AzTEC) that once there are many detectors, there are multiple correlations between subsets of the detectors, as well as an overall commonmode term. Consequently, one sees correlations between contiguous blocks of detectors on the array, or sets of detectors which share amplifiers or are otherwise coupled through the electronics. Provided that these correlations can be investigated and their behavior modelled, it is straightforward to extend the SANEPIC approach to deal with several distinct sources of correlated noise. Hence we expect the SANEPIC approach to be applicable to future instruments such as SCUBA2, SPIRE, ACT, Planck HFI and others. There are also many experiments being planned which use large detector arrays to perform sensitive polarization measurements, and we see no reason why the SANEPIC approach could not also be extended to polarimetry.
Footnotes
 affiliation: Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
 affiliation: Laboratoire APC, 10, rue Alice Domon et Léonie Duquet 75205, Paris, France
 affiliation: patanchon@apc.univparis7.fr
 affiliation: Department of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff, CF24 3AA, UK
 affiliation: Jet Propulsion Laboratory, Pasadena, CA 911098099
 affiliation: Observational Cosmology, MS 5933, California Institute of Technology, Pasadena, CA 91125
 affiliation: Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
 affiliation: Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104
 affiliation: Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104
 affiliation: Department of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff, CF24 3AA, UK
 affiliation: Department of Physics, University of Miami, 1320 Campo Sano Drive, Carol Gables, FL 33146
 affiliation: Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
 affiliation: Department of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff, CF24 3AA, UK
 affiliation: Instituto Nacional de Astrofísica Óptica y Electrónica (INAOE), Aptdo. Postal 51 y 72000 Puebla, Mexico
 affiliation: Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104
 affiliation: Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
 affiliation: Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
 affiliation: Department of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON Â M5S 3H4, Canada
 affiliation: Department of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff, CF24 3AA, UK
 affiliation: Department of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON Â M5S 3H4, Canada
 affiliation: Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada
 affiliation: Istituto di Radioastronomia, Largo E. Fermi 5, I50125, Firenze, Italy
 affiliation: University of Puerto Rico, Rio Piedras Campus, Physics Dept., Box 23343, UPR station, San Juan, Puerto Rico
 affiliation: Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada
 affiliation: Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104
 affiliation: Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada
 affiliation: Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104
 affiliation: Department of Physics, Brown University, 182 Hope Street, Providence, RI 02912
 affiliation: Department of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff, CF24 3AA, UK
 affiliation: Department of Physics, Brown University, 182 Hope Street, Providence, RI 02912
 affiliation: Department of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON Â M5S 3H4, Canada
 affiliation: Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada
 slugcomment: To appear in the Astrophysical Journal
 e.g. the 1994 article by J.R. Shewchuk: http://www.cs.cmu.edu/ ~quakepapers/painlessconjugategradient.pdf
References
 Ashdown, M., et al. 2007a, A&A, 467, 761
 Ashdown, M., et al. 2007b, A&A, 471, 361
 Benoit A., et al, 2002, Astropart. Phys., 17, 101
 Borrill, J. D. 1999, preprint, astroph/9911389
 Chapin, E. L. C. et al. 2007, ApJ, submitted
 Crill B.P., et al., 2003, ApJS, 148, 527
 Delabrouille, J. 1998, A&A Supl. 127, 555
 Delabrouille, J., Cardoso, J. F, Patanchon, G. 2003, MNRAS, 346, 1089
 Devlin, M., et al. 2004, Proc. SPIE, 5498, 42
 Doré, O., et al. 2001, A&A, 374, 358
 Dupac, X., Giard, M., 2002, MNRAS, 330, 497
 Ferreira, P. G., Jaffe, A. H. 2000, MNRAS, 312, 89
 de Gasperis, G., et al. 2005, A&A, 436, 1159
 Hamilton, J.Ch. 2003, preprint, astroph/0310787
 Hinshaw, G., et al. 2003, ApJS, 148, 63
 Janssen, M. A., et al. 1996, preprint, astroph/9602009
 Jenness, T., Lightfoot, J. F., Holland, W. S. 1998, Proc. SPIE, 3357, 548
 Keihänen, E., KurkiSuonio, H., Poutanen, T. 2005, MNRAS, 360, 390
 MacíasPérez, J., et al., A&A, 467, 1313
 Maino, D., et al. 2002, A&A, 387, 356
 MivilleDeschênes, M.A., Lagache, G., Boulanger, F., Puget, J.L., 2007, A&A, 469, 595
 Natoli, P., et al. 2001, A&A, 372, 346
 Oliver, S., et al. 2000, MNRAS, 316, 749
 Pascale, E. et al. 2007, ApJ, submitted
 Poutanen, T., et al. 2006, MNRAS, 449, 1311
 Prunet, S., Netterfield, C. B., Hivon, E., Crill., B. P. 2000, in “Proceedings of the XXXV Rencontres de Moriond (http://moriond.in2p3.fr/J00/ProcMJ2000/), preprint, astroph/0006052
 Prunet, S., et al. 2001, in “Mining the sky”, ed. A.J. Banday, S. Zaroubi, M. Bartelmann, ESO Astrophysics Symposia, SpringerVerlag, p
 Stompor, R., et al. 2002, Phys. Rev. D, 022003
 Tegmark, M. 1997, ApJ, 480, 87
 Tegmark, M., et al., 2000, ApJ, 541, 535
 Truch, M. et al. 2007, ApJ, submitted
 Wright, E. L., Hinshaw, G., Bennett, C. L. 1996, ApJ, 458, 53
 Yvon, D., Mayet, F. 2005, A&A, 436, 729