Calculating the angular power spectrum of a section of sky To calculate the angular power spectrum $C_l$ of the whole sky, one uses the variance of the coefficients of the spherical harmonics in the temperature fluctuation field. I.e.
$$C_l = \frac{1}{2l+1}\sum^{l}_{m=-l} (a_{lm}^{*} a_{lm})$$
How does one calculate the $C_l$ of a sub-set of the sky? Does one still decompose the temperature field in this sub-set into spherical harmonics and use the corresponding $a_{lm}$? This would seem inappropriate to me as the spherical harmonics are defined across the whole sky (surface of sphere) and would surely be a bad basis for a decomposition of the temperature field from a small section of the sky?
 A: Indeed spherical harmonics are inappropriate, since they are not orthogonal on the restricted domain. This is particularly noticeable in small-scale surveys like ACT and BOOMERanG, but even "full-sky" surveys mask bad data. COBE for instance masked out the entire galactic plane, so the problem has been known since then.

The solution presented by Górsky (1994 ApJ, 430, L85) is to Gram-Schmidt orthogonalize the spherical harmonics according to the inner product on the restricted domain. The procedure is as follows:


*

*Order the spherical harmonics $Y_{\ell m}$ into a column vector $y$, with $Y_{\ell m}$ being element $\ell^2 + \ell + 1 + m$ of $y$.

*Apply weights to each $y_i$ based on pixelization effects, resulting in $\tilde{y}_i$.

*Construct the matrix $W$ according to $W_{ij} = \langle \tilde{y}_i, \tilde{y}_j \rangle_\text{cut sky}$.

*Choleski-decompose $W = L L^\mathrm{T}$.

*Define a new vector of basis functions $\psi_i$, where $\psi = L^{-1} \tilde{y}$.


The resulting set is orthonormal: $\langle \psi_i, \psi_j \rangle_\text{cut sky} = \delta_{ij}$. The relation between the bases is such that if $f = \sum_i a_i y_i$ and $f\rvert_\text{cut sky} = \sum_i c_i \psi_i$, then $c_i = L_{ji} a_j$ and $a_i = L^{-1}_{ji} c_j$. As the author points out, the lower-triangular nature of $L$ prevents the bases from scrambling each other too much.
My understanding is that this method does not scale particularly well with the number of pixels $N$ in the image (methods like this tend to be $\mathcal{O}(N^2)$ or $\mathcal{O}(N^3)$).

An alternative is to use the restricted domain naively and then undo the statistical average effect of non-orthogonality. This was done in Hivon et al. (2002 ApJ, 567, 2).
Your theory predicts the values $C_\ell = 1/(2\ell+1) \cdot \sum_{m=-\ell}^\ell \lvert a_{\ell m} \rvert ^2$ for some function $f = \sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell a_{\ell m} Y_{\ell m}$. You measure $f$ in some places, with varying certainty, and with noise and other instrumental effects, and you want to extract the best values for $C_\ell$ from your data.
The suggested procedure can be summarized as follows:


*

*Define your weighting function $W$ on the sky (which can vanish in places). 

*Do the naive spherical harmonic decomposition of your data: $\tilde{a}_{\ell m} = \int_\text{sky} f W Y_{\ell m}^* \,\mathrm{d}\Omega$.

*Bin over $m$ as usual: $\tilde{C}_\ell = 1/(2\ell+1) \cdot \sum_{m=-\ell}^\ell \lvert \tilde{a}_{\ell m} \rvert^2$.

*Decompose the weighting function: $w_{\ell m} = \int_\text{sky} W Y_{\ell m}^* \,\mathrm{d}\Omega$.

*Calculate the weighting function power spectrum $\mathcal{W}_\ell = 1/(2\ell+1) \cdot \sum_{m=-\ell}^\ell \lvert w_{\ell m} \rvert^2$.

*Calculate the matrix elements
 $$ M_{\ell_1\ell_2} = \frac{2\ell_2+1}{4\pi} \sum_{\ell_3} (2\ell_3+1) \mathcal{W}_{\ell_3} \begin{pmatrix} \ell_1 & \ell_2 & \ell_3 \\ 0 & 0 & 0 \end{pmatrix}^2 $$
 using the Wigner 3-j symbol.

*The theoretical prediction relates to the observed values in a statistically averaged sense: $\langle \tilde{C}_\ell \rangle = \sum_{\ell'} M_{\ell\ell'} \langle C_{\ell'} \rangle$.


The paper goes on to discuss where beam width, filtering, and noise fit into the analysis, modifying how step (7) is used to extract the maximum likelihood estimators for $C_\ell$ given $\tilde{C}_\ell$.
