  We propose and implement a fast, universally applicable method for extracting
the angular power spectrum C_l from CMB temperature maps by first estimating
the correlation function \xi(\theta). Our procedure recovers the C_l's using
N^2 (but potentially N logN), operations, where N is the number of pixels. This
is in contrast with standard maximum likelihood techniques which require N^3
operations. Our method makes no special assumptions about the map, unlike
present fast techniques which rely on symmetries of the underlying noise
matrix, sky coverage, scanning strategy, and geometry. This enables for the
first time the analysis of megapixel maps without symmetries. The key element
of our technique is the accurate multipole decomposition of \xi(\theta). The
C_l error bars and cross-correlations are found by a Monte-Carlo approach. We
applied our technique to a large number of simulated maps with Boomerang sky
coverage in 81000 pixels. We used a diagonal noise matrix, with approximately
the same amplitude as Boomerang. These studies demonstrate that our technique
provides an unbiased estimator of the C_l's. Even though our method is
approximate, the error bars obtained are nearly optimal, and converged only
after few tens of Monte-Carlo realizations. Our method is directly applicable
for the non-diagonal noise matrix. This, and other generalizations, such as
minimum variance weighting schemes, polarization, and higher order statistics
are also discussed.

Astrophysics