© ESO, 2012

1. Introduction

Measurements of the cosmic microwave background (CMB) anisotropies are powerful cosmological probes. In the currently favored cosmological model, with the nearly Gaussian-distributed curvature perturbations, almost the entire statistical information is contained in the CMB angular power spectrum. The observed quantity on the sky is generally the CMB temperature anisotropy $\mathrm{\Theta }\mathrm{\left(}\begin{array}{}{}_{\mathit{􏿻}}\\ \mathit{p}\end{array}\mathrm{\right)}$ in direction $\begin{array}{}{\mathit{􏿻}}_{}\\ \mathit{p}\end{array}$, which is described as $\mathit{T}\mathrm{\left(}\begin{array}{}{}_{\mathit{􏿻}}\\ \mathit{p}\end{array}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{CMB}}\mathrm{[}\mathrm{1}\mathrm{+}\mathrm{\Theta }\mathrm{(}\begin{array}{}{}_{\mathit{􏿻}}\\ \mathit{p}\end{array}\mathrm{\left)}\mathrm{\right]}$. This field is expanded on the spherical harmonic functions as $\begin{array}{ccc}& & \mathrm{\Theta }\mathrm{\left(}\begin{array}{}\mathit{􏿻}\\ \mathit{p}\end{array}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell }\mathrm{=}\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }}^{\mathit{\ell }}\mathit{a}\mathrm{\left[}\mathit{\ell ,m}\mathrm{\right]}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\begin{array}{}\mathit{􏿻}\\ \mathit{p}\end{array}\mathrm{\right)}\mathit{,}\\ & & \mathrm{where}\mathit{a}\mathrm{\left[}\mathit{\ell ,m}\mathrm{\right]}\mathrm{=}{\mathrm{\int }}_{S^{\mathrm{2}}}\mathrm{\Theta }\mathrm{\left(}\begin{array}{}\mathit{􏿻}\\ \mathit{p}\end{array}\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\left(}\begin{array}{}\mathit{􏿻}\\ \mathit{p}\end{array}\mathrm{\right)}\mathrm{d}\begin{array}{}\mathit{􏿻}\\ \mathit{p}\end{array}\mathit{,}\end{array}$where S² ⊂ R³ is the unit sphere, ℓ is the multipole moment, which is related to the angular size on the sky as ℓ ~ 180°/θ and m is the phase ranging from  − ℓ to ℓ. The a[ℓ,m]  are the spherical harmonic coefficients of the (noise-free) observed sky. For a Gaussian random field, the mean and covariance are sufficient statistics, meaning that they carry the whole statistical information of the field. If the random field has zero mean, E(a₀₀) = 0 and the expansion can be started at ℓ = 2, neglecting the dipole terms, i.e., ℓ = 1¹. For , the triangular array (a[ℓ,m] )_ℓ,m represents zero-mean, complex-valued random coefficients, with variance $E\mathrm{\left(}\mathrm{\right|}\mathit{a}\mathrm{\left[}\mathit{\ell ,m}\mathrm{\right]}{\mathrm{|}}^{\mathrm{2}}\mathrm{)}\mathrm{=}\mathit{C}\mathrm{[}\mathit{\ell }\mathrm{]}\mathit{>}\mathrm{0}\mathit{,}$(3)where C[ℓ]  is the CMB angular power spectrum, which only depends on ℓ because of the isotropy assumption. Therefore, from Eq. (3), an unbiased estimator of C[ℓ]  is given by the empirical power spectrum $\begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}}{\left|\mathit{a}\mathrm{\left[}\mathit{\ell ,m}\mathrm{\right]}\right|}^{\mathrm{2}}\mathit{.}$(4)Furthermore, becasue the random field is stationary, the spherical harmonic coefficients are uncorrelated, $E\mathrm{\left(}\mathit{a}\mathrm{\right[}\mathit{\ell ,m}\mathrm{\left]}{\mathit{a}}^{\mathrm{\ast }}\mathrm{\right[}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}{\mathit{m}}^{\mathrm{\prime }}\mathrm{\left]}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathit{.}$(5)Since they are Gaussian, they are also independent. The angular power spectrum depends on the cosmological parameters through an angular transfer function T_ℓ(k) as $\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{\int }\frac{\mathrm{d}\mathit{k}}{\mathit{k}}{\mathit{T}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(6)where k defines the scale and P(k) is the primordial matter power spectrum.

Making accurate measurements of this power spectrum has been one of the main goals of cosmology in the past two decades. We have seen a range of ground- and balloon-based experiments, such as Acbar (Reichardt et al. 2009) and CBI (Readhead et al. 2004), as well as satellite experiments, such as WMAP (Bennett et al. 2003) and the recently launched satellite Planck². All these experiments produce a temperature map of the sky from which the CMB power spectrum is obtained. The estimation of the power spectrum from CMB experiments is important because this spectrum is a way to estimate the cosmological parameters that describe the Universe. General methods for extracting this spectrum from a N_pix-map, with nonuniform coverage and correlated noise are quite expensive and time-consuming especially for Planck, where we will be dealing with N_pix = 5 × 10⁷. All these experiments estimate the CMB angular power spectrum from a sky map, which is a realization of the underlying true power spectrum; no matter how much the experiments improve, we are still limited to an accuracy within the cosmic variance. This means that even if we had a perfect experiment (i.e. with zero instrumental noise), we would not be able to recover a perfect power spectrum, because of the cosmic variance limit.

In this paper we investigate the possibility of estimating the true underlying power spectrum from a realized spectrum; an estimation of the true power spectrum without a need to know the cosmological parameters. For this we exploit the sparsity properties of the CMB power spectrum, and capitalize on it to propose an estimator of the theoretical power spectrum.

The idea of sparsity in different dictionaries has been used previously. For example, Mukherjee & Wang (2004) used wavelets to estimate the level of non-Gaussianity in the first year of the WMAP data. In Mukherjee & Wang (2003) wavelets were used to estimate the primordial power spectrum. Faÿ et al. (2008) have used wavelets to estimate the power spectrum from a CMB map, as an alternative to the MASTER method of Hivon et al. (2002). There have also been previous attempts to smooth the CMB power spectrum, for instance spline-fitting, where the smoothed power spectrum was used for visual aids (Oh et al. 1999). Also, Aghamousa et al. (2012) have smoothed WMAP 1-3-5-7 power spectra to investigate their evolution.

We used the sparsity of the CMB power spectrum as a key ingredient to estimate the theoretical power spectrum without having to know the cosmological parameters; this estimate will not belong to a set of possible theoretical power spectra (i.e., all C[ℓ]  that can be obtained by CAMB³ by varying the cosmological parameters). Instead, this estimation should be useful for other applications, such as:

• Monte Carlo: we may need to conduct Monte Carlo simulations insome applications without assuming the cosmologicalparameters.

• Wiener filtering: Wiener filtering is often used to filter the CMB map and requires the theoretical power spectrum as an input. We may not wish to assume any cosmology at this stage of the processing.

• Some estimators (weak lensing, integrated SachsÐWolfe (ISW), etc.) require the theoretical power spectrum to be known. Using a data-based estimate of the theoretical C[ℓ]  could be an interesting alternative, or at least a good first guess in an iterative scheme where the theoretical C[ℓ]  is required for determining the cosmological parameters.

Paper content

Section 2 introduces the concepts of sparse representation and its applications to the CMB power spectrum. In Sect. 3 we explain how sparsity is used to propose an estimator of the theoretical power spectrum. Experimental results are described and discussed in Sect. 4, and Sect. 5 is devoted to a discussion and comparison of the moving average estimator. In Sect. 6, conclusions are drawn and potential perspectives are stated.

2. Sparsity of the CMB power spectrum

2.1. A brief tour of sparsity

A signal X = (X[1] ,...,X[N] ) considered as a vector in R^N, is said to be sparse if most of its entries are equal to zero. If k number of the N samples are not equal to zero, where k ≪ N, then the signal is said to be k-sparse. If only a few of the entries have high values and the rest are zero or close to zero, the signal is said to be weakly sparse (or compressible). With a slight abuse of terminology we call compressible signals sparse. Generally, signals are not sparse in direct space, but can be sparsified by transforming them to another domain. For example, sin(x) is 1-sparse in the Fourier domain, while it is clearly not sparse in the original one. In the so-called sparsity synthesis model, a signal can be represented as the linear expansion $\mathit{X}\mathrm{=}\mathrm{\Phi }\mathit{\alpha }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{T}}{\mathit{\phi }}_{\mathit{i}}\mathit{\alpha }\mathrm{\left[}\mathit{i}\mathrm{\right]}\mathit{,}$(7)where α[i]  are the synthesis coefficients of X, Φ = (φ₁,   ...,   φ_T) is the dictionary, and φ_i are called the atoms (elementary waveforms) of the dictionary Φ. In the language of linear algebra, the dictionary Φ is a N × T matrix whose columns are the normalized atoms, supposed here to be normalized to a unit ℓ₂-norm, i.e. $\mathrm{\forall }\mathit{i}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,T}\mathrm{\right]}\mathit{,}\parallel {\mathit{\phi }}_{\mathit{i}}\parallel \begin{array}{}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}{\sum }_{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\left|{\mathit{\phi }}_{\mathit{i}}\mathrm{\left[}\mathit{n}\mathrm{\right]}\right|}^{\mathrm{2}}\mathrm{=}\mathrm{1}$⁴. A function can be decomposed in many dictionaries, but the best dictionary is that with the sparsest (most economical) representation of the signal. In practice, it is convenient to use dictionaries with fast implicit transform (such as Fourier transform, wavelet transform, etc.) which allow us to directly obtain the coefficients and reconstruct the signal from these coefficients using fast algorithms that run in linear or almost linear time (unlike matrix-vector multiplications). The Fourier, wavelet, and discrete cosine transforms are certainly the most well known dictionaries. A comprehensive account of sparsity and its applications can be found in the monograph Starck et al. (2010).

2.2. Which dictionary for the theoretical CMB power spectrum?

We investigated the sparsity of the CMB power spectrum in two different dictionaries, both having a fast implicit transform: the wavelet transform (WT) and the discrete cosine transform (DCT).

Fig. 1 Theoretical CMB power spectrum along with the reconstructed power spectra, using the DCT and WT dictionaries. The panels show the reconstructions for different fractions of the coefficients. The inner plots show the differences between the actual and the reconstructed power spectra. Both dictionaries suffer from boundary effects, but this is more severe for DCT because the corresponding atoms are not compactly supported. The power spectrum that is decomposed onto the two dictionaries is in the form ℓ(ℓ + 1)C[ℓ] /2π.

Fig. 2 Non-linear approximation (NLA) error curves for the two dictionaries. Below 1% the DCT curve is dropping faster, which means it is performing better. However, past  ~ 2% the DCT curve flattens off, while WT decreases to  ~ 0 very quickly.

Figure 1 shows an angular power spectrum (calculated with CAMB, Lewis et al. 2000; with WMAP7, Larson et al. 2011, parameters) along with the DCT- and WT-reconstructed power spectra with a varying fraction of the largest transform coefficients retained in the reconstruction. The inner plots show the difference between the actual power spectrum and the reconstructed ones. Clearly, with only few coefficients the shape of the power spectrum is correctly reconstructed in both dictionaries. The height and the position of the peaks and troughs are very important here because the estimation of the cosmological parameters heavily relies on these characteristics of the power spectrum. The best domain would be the one with the sparsest representation and yet the most accurate representation of the power spectrum. Let C[ℓ] ^(M) be its best M-term approximation, i.e., obtained by reconstructing from the M-largest (in magnitude) coefficients of C[ℓ]  in a given domain. To compare the WT and DCT dictionaries, we plot the resulting non-linear approximation (NLA) error curve in Fig. 2, which shows the reconstruction error E_M as a function of M, the number of retained coefficients ${\mathit{E}}_{\mathit{M}}\mathrm{=}\frac{{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{-}{\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\mathrm{2}}}{{\parallel \mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\parallel }_{\mathrm{2}}}\mathrm{\times }\mathrm{100.}$(8)As M increases we approach the complete reconstruction, and the error reaches zero when all coefficients have been used. Usually the domain with the steepest E_M curve is the sparsest domain. In this case, though, both dictionaries show a very similar behavior. There is only a small window in the coefficients for which DCT performs better than WT. However, DCT flattens after using  ~ 1% of the coefficients and does not improve the reconstruction until many coefficients have been used.

Both dictionaries seem to suffer from boundary problems at low and high ℓs. This can be solved for high ℓs because one can always perform the reconstruction beyond the desired ℓ. For low ℓs it can be solved by different means, for instance by extrapolating the spectrum. Note that the boundary problems are more severe in the DCT domain than in the WT; this is because DCT atoms are not compactly supported.

Fig. 3 A simulated CMB power spectrum along with the reconstructed spectra, using the DCT and WT dictionaries. The black solid line is the true underlying power spectrum from which the simulations were made. The blue and red dots show the simulated and the reconstructed power spectra respectively. With only 1% of the coefficients, DCT can recover the input power spectrum (i.e. the black solid line) very well, recovering the peaks and troughs accurately. Unlike DCT, WT seems to have difficulties in recovering the peaks and troughs. The inner plots shows the NMSE curves.

Next we investigated the sparsity of a set of realized spectra in the two dictionaries. We simulated 100 maps from the theoretical power spectrum used above and estimated their power spectra using Eq. (4). As before, we decomposed each realization in the DCT and WT dictionaries and reconstructed keeping increasing fractions of the highest coefficients. At this stage, it is important to note that because we are dealing with the empirical power spectrum, we are no longer in an approximation setting but instead in an estimation one. Indeed, the empirical power spectrum can be considered a noisy version of the true one. Intuitively, reconstructing from a very small fraction of high coefficients will reject most of the noise (low estimator variance) but at the price of retaining only a small fraction of the true spectrum coefficients (strong bias). The opposite is true when a large proportion of coefficients is kept in the reconstruction. Therefore, there is a threshold value that will entail a bias-variance trade-off, hence minimizing the estimation risk. This is exactly the idea underlying thresholding estimators in sparsifying domains.

Fig. 4 Average of the DCT-reconstructed power spectra when different fractions of the most significant coefficients used in the reconstruction. The inner plot shows the NMSE curve for this average. The minimum of the curve is at less than 10% of the coefficients, meaning that the true CMB spectrum can be recovered with less than 10% of the coefficients while ensuring a good bias-variance trade-off.

This discussion is clearly illustrated by the inner plots of Fig. 3, which shows the normalized mean-square error (NMSE) defined as ${\mathrm{NMSE}}_{\mathit{M}}\mathrm{=}\frac{{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{-}{\begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\mathrm{2}}}{{\parallel \mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\parallel }_{\mathrm{2}}}\mathrm{\times }\mathrm{100}\mathit{,}$(9)as a function of the fraction of coefficients used in the reconstruction. The error is large when only a few coefficients are used. As more coefficients are included, one starts to recover the main (i.e., the general shape of the spectrum) features of the power spectrum. With more coefficients, more noise enters the estimation and the error increases again. The NMSE curve shows a clear minimum at which the underlying true power spectrum is best recovered.

Despite the differences in the performance of the two dictionaries, the minima of the NMSE have about the same proportions of the coefficients. This is because the NMSE reflects a global behavior. On the one hand, although the DCT can recover the features of the spectrum correctly, it is less smooth than WT. On the other hand, the WT cannot reconstruct the proper shape of the power spectrum, but provides a smoother estimate.

Figures 4 and 5 show the same results for an average over the 100 realizations. Evidently, on average DCT performs very well in recovering the features of the CMB spectrum. Indeed, the minimum of NMSE curve is at a lower proportion of coefficients for DCT than WT. This is because the small noisy features of the DCT-reconstructed spectra cancel out in the averaging, while the wrong recovery of the shape of the spectrum by WT does not. In a nutshell, DCT seems to perform better in reconstructing the true underlying CMB power spectrum from its realizations.

To summarize, from the above discussion, we conclude the following:

• the CMB power spectrum is very sparse in both DCT and WTdictionaries, although their sparsifying capabilities aredifferent;

• DCT recovers global features of spectrum (i.e., the peaks and troughs), while WT recovers localized features;

• WT recovers more localized (noisy) features than the global ones for realizations, while the DCT concentrates on the global features.

In the next section, these complementary capabilities of the DCT and WT transforms will be combined to propose a versatile way for adaptively estimating the theoretical power spectrum from a single realization of it.

3. Sparse reconstruction of the theoretical power spectrum

We begin with the simple model where the observed signal Y is contaminated by a zero-mean white Gaussian noise, Y = X + ε, where X is the signal of interest and . Sparse recovery with an analysis-type sparsity prior amounts to finding the solution of the following problem: $\underset{\mathit{X}}{\min }{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}{\mathrm{\Phi }}^{\mathit{T}}\mathit{X}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\mathrm{1}} \mathit{s}\mathit{.}\mathit{t}\mathit{.} {\parallel \mathit{Y}\mathrm{-}\mathit{X}\parallel }_{\mathrm{2}}\mathrm{\le }\mathit{\delta ,}$(10)where Φ^TX represents the transform coefficients of X in the dictionary Φ, and δ controls the fidelity to the data and obviously depends on the noise standard deviation σ.

We now turn to denoising the power spectrum from one empirical realization of it. In this case, however, the noise is highly non-Gaussian and needs to be treated differently. Indeed, as we will see in the next section, the empirical power spectrum will entail a multiplicative χ²-distributed noise with a number of degrees of freedom that depends on ℓ. That is, the noise has a variance profile that depends both on the true spectrum and ℓ. We therefore need to stabilize the noise on the empirical power spectrum prior to estimation, using a variance stabilization transform (VST). Hopefully, the latter will yield stabilized samples that have (asymptotically) constant variance, say 1, irrespective of the value of the input noise level.

Fig. 5 Same as Fig. 4, but for the WT domain. It seems that WT cannot perform well for simulated power spectra, compared to the DCT domain; the minimum of the NMSE curve is at more than 35% – while the curve is quite flat after  ~ 20%.

3.1. Variance stabilizing transform

In the statistical literature the problem of removing the noise from an empirical power spectrum is called periodogram denoising (Donoho 1993). Komm et al. (1999) approximated the noise with a correlated Gaussian noise model and derived a threshold at each wavelet scale using the MAD (median of absolute deviation) estimator. A more elegant approach was proposed in Donoho (1993) and Moulin (1994), where the so-called Wahba VST was used. This VST is defined as $\mathrm{𝒯}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}\left(\log \mathit{X}\mathrm{+}\mathit{\gamma }\right)\frac{\sqrt{\mathrm{6}}}{\mathit{\pi }}\mathit{,}$(11)where γ = 0.57721... is the Euler-Mascheroni constant. After the VST, the stabilized samples can be treated as if the noise contaminating them were white Gaussian noise with unit variance.

We will take a similar path here, generalizing the above approach to the case of the angular power spectrum. Indeed, from Eq. (4), one can show that under mild regularity assumptions on the true power spectrum, $\begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\stackrel{\mathit{d}}{\mathrm{\to }}\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathit{Z}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathit{,}\mathrm{where}\mathrm{\forall }\mathit{\ell }\mathrm{\ge }\mathrm{2}\mathit{,}\hspace{0.17em}\mathrm{2}\mathit{LZ}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{~}{\mathit{\chi }}_{\mathrm{2}\mathit{L}}^{\mathrm{2}}\mathit{,L}\mathrm{=}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathit{,}$(12)where $\stackrel{\mathit{d}}{\mathrm{\to }}$ means convergence in distribution. From Eq. (12), it is appealing then to take the logarithm to transform the multiplicative noise Z into an additive one. The resulting log -stabilized empirical power spectrum reads ${\mathit{C}}^{\mathit{s}}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{:}\mathrm{=}{\mathrm{𝒯}}_{\mathit{\ell }}\mathrm{\left(}\begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\right[}\mathit{\ell }\mathrm{\left]}\mathrm{\right)}\mathrm{=}\log \begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{-}{\mathit{\mu }}_{\mathit{L}}\mathrm{=}\log \mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{+}\mathit{\eta }\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathit{,}$(13)where η[ℓ] : = log Z[ℓ]  − μ_L, L = 2ℓ + 1. Using the asymptotic results from Bartlett & Kendall (1946) on the moments of log  − χ² variables, it can be shown that ${\mathit{\mu }}_{\mathit{L}}\mathrm{=}{\mathit{\psi }}_{\mathrm{0}}{}^{\mathrm{\right(}}\mathit{L}^{\mathrm{\left)}}\mathrm{-}\log \mathit{L}$, E(η[ℓ] ) = 0 and ${\mathit{\sigma }}_{\mathit{L}}^{\mathrm{2}}\mathrm{=}\mathrm{Var}\left[\mathit{\eta }\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\right]\mathrm{=}{\mathit{\psi }}_{\mathrm{1}}{}^{\mathrm{\right(}}\mathit{L}^{\mathrm{\left)}}$, where ψ_m(t) is the standard polygamma function, ${\mathit{\psi }}_{\mathit{m}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{d}}^{\mathit{m}\mathrm{+}\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathit{m}\mathrm{+}\mathrm{1}}}\log \mathrm{\Gamma }\mathrm{\left(}\mathit{t}\mathrm{\right)}$.

We can now consider the stabilized C^s[ℓ]  as noisy versions of the log C[ℓ] , where the noise is zero-mean additive and independent. Owing to the central limit theorem, the noise tends to Gaussian with variance ${\mathit{\sigma }}_{\mathit{L}}^{\mathrm{2}}$ as ℓ increases. At low ℓ, normality is only an approximation. Indeed, the noise η[ℓ]  has a probability density function of the form ${\mathit{p}}_{\mathit{\eta }}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{L}{\mathrm{)}}^{\mathit{L}}}{{\mathrm{2}}^{\mathit{L}}\mathrm{\Gamma }\mathrm{\left(}\mathit{L}\mathrm{\right)}}\exp \mathrm{[}\mathit{L}{\mathrm{(}\mathit{\ell }\mathrm{+}{\mathit{\mu }}_{\mathit{L}}\mathrm{-}{\mathrm{e}}^{\mathit{\ell }\mathrm{+}{\mathit{\mu }}_{\mathit{L}}}\mathrm{)}}^{\mathrm{]}}\mathit{,}$(14)which might be used to estimate the thresholds in the wavelet domain.

To standardize the noise, we slightly modify the VST (13)to the normalized form ${\mathit{C}}^{\mathit{s}}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{:}\mathrm{=}{\mathrm{𝒯}}_{\mathit{\ell }}\mathrm{\left(}\begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\right[}\mathit{\ell }\mathrm{\left]}\mathrm{\right)}\mathrm{=}\frac{\log \begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{-}{\mathit{\mu }}_{\mathit{L}}}{{\mathit{\sigma }}_{\mathit{L}}}\mathrm{=}{\mathit{X}}^{\mathit{s}}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{+}\mathit{\epsilon }\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathit{,}$(15)where now the noise ε[ℓ]  is a zero-mean (asymptotically) Gaussian with unit variance, and X^s[ℓ] : = log C[ℓ] /σ_L. It can be checked that the Wahba VST (11)is a specialization of Eq. (15) to L = 0.

In the following, we will use the operator notation for the VST that applies Eq. (15) entry-wise to each X[ℓ] , and ℛ(X) its inverse operator, i.e. ℛ(X): = (ℛ_ℓ(X[ℓ] ))_ℓ with ℛ_ℓ(X[ℓ] ) = exp(σ_LX[ℓ] ).

3.2. Signal detection in the wavelet domain

Without of loss of generality, we restrict our description here to the wavelet transform. The same approach applies to other sparsifying transforms, e.g., DCT.

To estimate the true CMB power spectrum from the wavelet transform, it is important to detect the wavelet coefficients that are “significant”, i.e. those wavelet coefficients whose absolute value is too high to be due to noise (cosmic variance + instrumental noise). Let w_j[ℓ]  be the wavelet coefficient of a signal Y at scale j and location ℓ. We define the multiresolution support M of Y as ${\mathit{M}}_{\mathit{j}}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}{\mathit{w}}_{\mathit{j}}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{issignificant,}\\ & \mathrm{if}{\mathit{w}}_{\mathit{j}}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{otherwise}\mathrm{.}\end{array}$(16)For Gaussian noise, it is easy to derive an estimate of the noise standard deviation σ_j at scale j from the noise standard deviation, which can be evaluated with good accuracy in an automated way (Starck & Murtagh 1998). To detect the significant wavelet coefficients, it suffices to compare the wavelet coefficients in magnitude |w_j[ℓ] | to a threshold level t_j. This threshold is generally taken to be equal to κσ_j, where κ ranges from 3 to 5. This means that a small magnitude compared to the threshold implies that the coefficients are very likely to be due to noise and are hence insignificant. This decision corresponds to the hard-thresholding operator $\begin{array}{ccc}\begin{array}{c}\\ \\ \end{array}& & \end{array}$(17)To summarize, the multiresolution support is obtained from the signal Y by computing the forward transform coefficients, applying hard thresholding, and recording the coordinates of the retained coefficients.

3.3. Power spectrum recovery algorithm

We now turn to the adaptive estimator of the true CMB power spectrum C[ℓ]  from its empirical estimate $\begin{array}{}{}_{\mathrm{􏽢}}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}$. Because we benefit from the (asymptotic) normality of the noise in the stabilized samples C^s[ℓ]  in Eq. (15), we are in position to easily construct the multiresolution support M of C^s as described in the previous section. Once the support M of significant coefficients has been determined, our goal is to reconstruct an estimate $\begin{array}{}{\mathrm{􏽥}}_{}\\ \mathit{X}\end{array}$ of the true power spectrum, known to be sparsely represented in some dictionary Φ(regularization), such that the significant transform coefficients of its stabilized version reproduce those of C^s (fidelity to data). Furthermore, because a power spectrum is a positive, a positivity constraint must be imposed. These requirements can be cast as seeking an estimate that solves the following constrained optimization problem: $\underset{\mathit{X}}{\min }\mathrm{\parallel }{\mathrm{\Phi }}^{\mathit{T}}\mathit{X}{\mathrm{\parallel }}_{\mathrm{1}} \mathrm{s}\mathit{.}\mathrm{t}\mathit{.} \{\begin{array}{c}\\ \\ \end{array}$(18)where  ⊙  stands for the Hadamard product (i.e., entry-wise multiplication) of two vectors. This problem has a global minimizer, which is bounded. However, in addition to non-smoothness of the l₁-norm and the constraints, the problem is also non-convex because of the VST operator . It is therefore far from easy to solve.

We propose the following scheme, which starts with an initial guess of the power spectrum X⁽⁰⁾ = 0, and then iterates for n = 0 to N_max − 1, $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\mathrm{􏽥}\\ \mathit{X}\end{array}\mathrm{=}\mathrm{ℛ}\mathrm{(}\mathrm{𝒯}\mathrm{(}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{)}\mathrm{+}\mathrm{\Phi }\mathit{M}\mathrm{\odot }{\mathrm{(}{\mathrm{\Phi }}^{\mathit{T}}{\mathrm{(}{\mathit{C}}^{\mathit{s}}\mathrm{-}\mathrm{𝒯}{\mathrm{(}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{)}}^{\mathrm{)}}}^{\mathrm{)}}}^{\mathrm{)}}\begin{array}{}\end{array}{\mathit{X}}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathrm{𝒫}}_{\mathrm{+}}\mathrm{(}\mathrm{\Phi } {\mathrm{ST}}_{{\mathit{\lambda }}_{\mathit{n}}}\mathrm{(}{\mathrm{\Phi }}^{\mathit{T}}\begin{array}{}\mathrm{􏽥}\\ \mathit{X}\end{array}{\mathrm{)}}^{\mathrm{)}}\mathit{,}\end{array}\begin{array}{}\end{array}\end{array}$(19)where denotes the projection onto the positive orthant and guarantees non-negativity of the spectrum estimator, ST_λn(w) = (ST_λn(w[i] ))_i is the soft-thresholding with threshold λ_n that applies term-by-term the shrinkage rule ${\mathrm{ST}}_{{\mathit{\lambda }}_{\mathit{n}}}\mathrm{\left(}\mathit{w}\mathrm{\right[}\mathit{i}\mathrm{\left]}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathrm{\left|}\mathit{w}\mathrm{\right[}\mathit{i}\mathrm{\left]}\mathrm{\right|}⩾{\mathit{\lambda }}_{\mathit{n}}\mathit{,}\\ & \mathrm{otherwise}\mathit{.}\end{array}$(20)Here, we have chosen a decreasing threshold with the iteration number n, λ_n = (N_max − n)/(N_max − 1). More details pertaining to this algorithm can be found in Starck et al. (2010).

Algorithm 1 summarizes the main steps of the sparse denoising algorithm. A similar approach was proposed in Starck et al. (2009) and Schmitt et al. (2010) for Poisson noise removal in 2D and 3D data sets.

3.4. Instrumental noise

In practice, the data are generally contaminated by an instrumental noise, and estimating the true CMB power spectrum C[ℓ]  from the empirical power spectrum $\begin{array}{}{}_{\mathrm{􏽢}}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}$ requires one to remove this instrumental noise. The instrumental noise is assumed to be stationary and independent of the CMB. We also assumed that we have access to the power spectrum of the noise, or we can compute the empirical power spectrum $\begin{array}{}{}^{\mathrm{􏽢}}\\ {\mathit{S}}_{\mathit{N}}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}$ of at least one realization, either from a JackKnife data map or from realistic instrumental noise simulations. The above algorithm can be adapted to handle this case after rewriting the optimizing problem as follows: $\underset{\mathit{X}}{\min }\mathrm{\parallel }{\mathrm{\Phi }}^{\mathit{T}}\mathit{X}{\mathrm{\parallel }}_{\mathrm{1}} \mathrm{s}\mathit{.}\mathrm{t}\mathit{.} \{\begin{array}{c}\\ \\ \end{array}$(21)

Thus, Eq. (19) becomes $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\mathrm{􏽥}\\ \mathit{X}\end{array}\mathrm{=}\mathrm{ℛ}(\mathrm{𝒯}\mathrm{(}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{+}\begin{array}{}\mathrm{􏽢}\\ {\mathit{S}}_{\mathit{N}}\end{array}\mathrm{)}\begin{array}{}\end{array}\mathrm{+}\mathrm{\Phi }\mathit{M}\mathrm{\odot }({\mathrm{\Phi }}^{\mathit{T}}({\mathit{C}}^{\mathit{s}}\mathrm{-}\mathrm{𝒯}\mathrm{(}{\mathit{X}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{+}\begin{array}{}\mathrm{􏽢}\\ {\mathit{S}}_{\mathit{N}}\end{array}\mathrm{)})))\mathrm{-}\begin{array}{}\mathrm{􏽢}\\ {\mathit{S}}_{\mathit{N}}\end{array}\begin{array}{}\end{array}{\mathit{X}}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathrm{𝒫}}_{\mathrm{+}}(\mathrm{\Phi }{\mathrm{ST}}_{{\mathit{\lambda }}_{\mathit{n}}}\mathrm{\left(}{\mathrm{\Phi }}^{\mathit{T}}\begin{array}{}\mathrm{􏽥}\\ \mathit{X}\end{array}\mathrm{\right)})\mathit{.}\end{array}\begin{array}{}\end{array}\end{array}$(22)Algorithm 1 can be modified accordingly.

Fig. 6 Theoretical CMB power spectrum (black line), the empirical power spectrum of one realization (blue dots), and the avegared of estimated power spectra (red line) using the TOUSI algorithm. The inner plots show a zoomed-in version.

Fig. 7 Cosmological parameters estimated from the true CMB power spectrum (black line) and the mean of the reconstructed power spectra (red line). The dashed line is the true input parameters, i.e., those used to calculate the theoretical power spectrum using CAMB.

3.5. Combining several dictionnaries

We have seen in Sect. 2.2 that the WT and DCT dictionaries had complementary benefits. Indeed, each dictionary is able to capture the features with shapes similar to its atoms well. More generally, assuming that we have D dictionaries Φ₁,...,Φ_D. Given a candidate signal Y, we can derive a support M_d associated to each dictionary Φ_d, for d ∈  { 1,...,D } . The optimization problem to solve now reads $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\underset{\mathit{X}}{\min }\mathrm{\parallel }{\mathrm{\Phi }}^{\mathit{T}}\mathit{X}{\mathrm{\parallel }}_{\mathrm{1}} \mathrm{s}\mathit{.}\mathrm{t}\mathit{.} \begin{array}{}\end{array}\{\begin{array}{c}\\ \\ \end{array}\end{array}\begin{array}{}\end{array}\end{array}$(23)Again, this is a challenging optimization problem. We propose to solve it by applying Eq. (22) successively and alternatively on each dictionary Φ_d. Algorithm 2 describes in detail the different steps.

4. Application: Monte Carlo simulations

We simulated 100 maps from a theoretical CMB power spectrum that was calculated by CAMB. The power spectra of these maps are equivalent to 100 realizations of the true CMB power spectrum – this realized spectrum is what we have access to in reality. Each of these 100 simulated spectra were run through the TOUSI algorithm, with the aim of recovering the theoretical spectrum from which these 100 spectra were simulated (i.e., the one that was calculated by CAMB).

Figure 6 shows an example of this; the empirical power spectrum of one realization (blue dots) that was fed into the TOUSI algorithm, the average of the 100 estimated power spectral using TOUSI (red line) and the input theoretical spectra (black line). The black line is the input theoretical spectrum that was calculated by CAMB, which is what we are trying to recover. The reconstruction of the peaks and troughs of the power spectrum by this algorithm is very impressive. This is very important because these features define the cosmological parameters. To additionally check the accuracy of these reconstructed spectra, we estimated a set of cosmological parameters from these spectra, using CosmoMC (Lewis & Bridle 2002). First, a set of cosmological parameters were estimated from the 100 simulated spectra. The results are shown in Fig. 7 as black solid lines. Then a “mean” reconstructed spectrum was calculated by averaging the 100 reconstructed spectra. This, in principle (i.e., if the algorithm has worked), should be an estimation of the true input spectrum with the same characteristics and the same cosmological parameters. To test this, we used this “mean” reconstructed spectrum to simulate another 100 maps and then 100 spectra. These simulated spectra were run through CosmoMC to estimate the same set of cosmological parameters. These are shown as red lines in Fig. 7. Clearly, the TOUSI algorithm can reconstruct the true underlying power spectrum with great accuracy in the cosmological parameters.

Fig. 8 Power spectrum estimate in the presence of instrumental noise. The blue dots show the empirical power spectrum of one realization with instrumental noise. Yellow dots show the estimated power spectrum of one of the simulated noise maps. Green dots show the the spectrum with the noise power spectrum removed. The black and red solid lines are the input and reconstructed power spectra, respectively. The inner plots show a zoomed-in version.

4.1. Data with instrumental noise

Here we present the performance of the TOUSI algorithm in the presence of instrumental noise. The noise maps were simulated using a theoretical (Planck level) noise power spectrum. They were added to the CMB maps simulated above and the power spectra of the combined maps were estimated using Eq. (4).

Figure 8 shows the reconstruction of the theoretical CMB spectrum in the presence of noise. The blue dots show the empirical power spectrum of one realization with instrumental noise. Yellow dots show the estimated power spectrum of one of the simulated noise maps. Green dots show the the spectrum with the noise power spectrum removed. The black and red solid lines are the input and reconstructed power spectra, respectively. The theoretical power spectrum can be reconstructed up to the point where the structure of the power spectrum has not been destroyed by the instrumental noise. In our case, having Planck level noise, this goes to ℓ up to 2500. Obviously, the TOUSI performs very well in reconstructing the input power spectrum even in the presence of instrumental noise.

4.2. Test on WMAP7 power spectrum

We tested our algorithm on real data. Figure 9 shows the application of our technique to the WMAP7 power spectrum. The method indeed works well up to ℓ of  ~800. After this point the instrumental noise becomes so dominant that the features of the spectrum are washed out and cannot be recovered.

5. Sparsity versus averaging

A very common approach to reduce the noise on the power spectrum is the moving average filter, i.e., average values in a given window, $\begin{array}{}\mathrm{􏽥}\\ {\mathit{C}}^{\mathit{A}}\end{array}\mathrm{\left[}\mathit{b}\mathrm{\right]}\mathrm{=}\frac{\mathrm{1}}{\mathit{b}\mathrm{(}\mathit{b}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\omega }}_{\mathit{b}}}\sum _{\mathit{\ell }\mathrm{=}\mathit{b}\mathrm{-}\frac{{\mathit{\omega }}_{\mathit{b}}}{\mathrm{2}}}^{\mathit{b}\mathrm{+}\frac{{\mathit{\omega }}_{\mathit{b}}}{\mathrm{2}}}\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\begin{array}{}\mathrm{􏽢}\\ \mathit{C}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathit{,}$(24)and the window size ω_b is increasing with ℓ. Here, we used window sizes of  { 1,2,5,10,20,50,100 } , respectively, for ℓ ranging from  { 2,11,31,151,421,1201,2501 }  to  { 10,30,150,420,1200,2500,3200 } , which have also been used in the framework of the Planck project in Leach et al. (2008).

Figure 10 shows a reconstruction of the power spectrum for TOUSI versus averaging. From the NMSE curve it is clear that our algorithm is much more efficient. Figure 11 shows the average error the 100 realizations as a function of ℓ$\mathit{E}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{=}\frac{\mathrm{1}}{\mathrm{100}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{100}}\mathrm{\parallel }\mathit{C}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}\mathrm{-}\begin{array}{}\mathrm{􏽥}\\ {\mathit{C}}_{\mathit{i}}\end{array}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}{\mathrm{\parallel }}_{\mathrm{2}}\mathit{,}$(25)where $\begin{array}{}{\mathrm{􏽥}}_{}\\ {\mathit{C}}_{\mathit{i}}\end{array}$ is the estimated power spectrum from the i-th realization. We display the errors for the spectra estimated by the empirical estimator (the realization, black dotted line), the averaging estimator (red dashed line) and TOUSI (solid blue line). The cosmic variance is overplotted as a solid black line. Evidently, the expected error is highly reduced when using the sparsity-based estimator.

Fig. 9 TOUSI algorithm applied to the WMAP7 power spectrum. The technique works well up to ℓ of  ~800, i.e., before the instrumental noise becomes dominant washes the features of the spectrum out.

Fig. 10 Mean denoised spectra with wavelets (blue solid line) and averaging (red dashed line) from 100 realizations. The inner plot shows the normalized error for both dictionaries.

Fig. 11 Mean error for the 100 realizations, for the realizations (black dotted line), the averaging denoising (red dashed line) and the sparse wavelet filtering (blue solid line). The inner plot shows a zoom between l = 2000 and l = 3000.

6. Conclusion

Measurements of the CMB anisotropies are powerful cosmological probes. In the currently favored cosmological model, with the nearly Gaussian-distributed curvature perturbations, almost the entire statistical information is contained in the CMB angular power spectrum. We have investigated the sparsity of the CMB power spectrum in two dictionaries, the DCT and WT. In both dictionaries the CMB power spectrum can be recovered with only a few percentages of the coefficients, meaning the spectrum is very sparse. The two dictionaries have different characteristics and can accommodate reconstructing different features of the spectra. The DCT can help recover the global features of the spectrum, while WT helps recover small localized features. The sparsity of the CMB spectrum in these two domains has helped us develop an algorithm, TOUSI, which estimates the true underlying power spectrum from a given realized spectrum. This algorithm uses the sparsity of the CMB power spectrum in both WT and DCT domains and takes the best from both worlds to obtain a highly accurate estimate from a single realization of the CMB power spectrum. This could be a replacement for CAMB if knowing the cosmological parameters is not necessary. The developed IDL code will be released with the next version of ISAP (Interactive Sparse astronomical data Analysis Packages) via the web site http://jstarck.free.fr/isap.html.

Acknowledgments

The authors thank Marian Douspis, Olivier Doré, and Amir Hajian for useful discussions. This work is supported by the European Research Council grant SparseAstro (ERC-228261).

References

All Figures

Fig. 1 Theoretical CMB power spectrum along with the reconstructed power spectra, using the DCT and WT dictionaries. The panels show the reconstructions for different fractions of the coefficients. The inner plots show the differences between the actual and the reconstructed power spectra. Both dictionaries suffer from boundary effects, but this is more severe for DCT because the corresponding atoms are not compactly supported. The power spectrum that is decomposed onto the two dictionaries is in the form ℓ(ℓ + 1)C[ℓ] /2π.

In the text

	Fig. 2 Non-linear approximation (NLA) error curves for the two dictionaries. Below 1% the DCT curve is dropping faster, which means it is performing better. However, past  ~ 2% the DCT curve flattens off, while WT decreases to  ~ 0 very quickly.
In the text

Fig. 3 A simulated CMB power spectrum along with the reconstructed spectra, using the DCT and WT dictionaries. The black solid line is the true underlying power spectrum from which the simulations were made. The blue and red dots show the simulated and the reconstructed power spectra respectively. With only 1% of the coefficients, DCT can recover the input power spectrum (i.e. the black solid line) very well, recovering the peaks and troughs accurately. Unlike DCT, WT seems to have difficulties in recovering the peaks and troughs. The inner plots shows the NMSE curves.

In the text

Fig. 4 Average of the DCT-reconstructed power spectra when different fractions of the most significant coefficients used in the reconstruction. The inner plot shows the NMSE curve for this average. The minimum of the curve is at less than 10% of the coefficients, meaning that the true CMB spectrum can be recovered with less than 10% of the coefficients while ensuring a good bias-variance trade-off.

In the text

	Fig. 5 Same as Fig. 4, but for the WT domain. It seems that WT cannot perform well for simulated power spectra, compared to the DCT domain; the minimum of the NMSE curve is at more than 35% – while the curve is quite flat after  ~ 20%.
In the text

	Fig. 6 Theoretical CMB power spectrum (black line), the empirical power spectrum of one realization (blue dots), and the avegared of estimated power spectra (red line) using the TOUSI algorithm. The inner plots show a zoomed-in version.
In the text

	Fig. 7 Cosmological parameters estimated from the true CMB power spectrum (black line) and the mean of the reconstructed power spectra (red line). The dashed line is the true input parameters, i.e., those used to calculate the theoretical power spectrum using CAMB.
In the text

Fig. 8 Power spectrum estimate in the presence of instrumental noise. The blue dots show the empirical power spectrum of one realization with instrumental noise. Yellow dots show the estimated power spectrum of one of the simulated noise maps. Green dots show the the spectrum with the noise power spectrum removed. The black and red solid lines are the input and reconstructed power spectra, respectively. The inner plots show a zoomed-in version.

In the text

	Fig. 9 TOUSI algorithm applied to the WMAP7 power spectrum. The technique works well up to ℓ of  ~800, i.e., before the instrumental noise becomes dominant washes the features of the spectrum out.
In the text

	Fig. 10 Mean denoised spectra with wavelets (blue solid line) and averaging (red dashed line) from 100 realizations. The inner plot shows the normalized error for both dictionaries.
In the text

	Fig. 11 Mean error for the 100 realizations, for the realizations (black dotted line), the averaging denoising (red dashed line) and the sparse wavelet filtering (blue solid line). The inner plot shows a zoom between l = 2000 and l = 3000.
In the text