  Szapudi et al (2001) introduced the method of estimating angular power
spectrum of the CMB sky via heuristically weighted correlation functions. Part
of the new technique is that all (co)variances are evaluated by massive Monte
Carlo simulations, therefore a fast way to measure correlation functions in a
high resolution map is essential. This letter presents a new algorithm to
calculate pixel space correlation functions via fast spherical harmonics
transforms. Our present implementation of the idea extracts correlations from a
MAP-like CMB map (HEALPix resolution of 512, i.e. $ \simeq 3 \times 10^6$
pixels) in about 5 minutes on a 500MHz computer, including $C_\ell$ inversion;
the analysis of one Planck-like map takes less then one hour. We use heuristic
window and noise weighting in pixel space, and include the possibility of
additional signal weighting as well, either in $\ell$ or pixel space. We apply
the new code to an ensemble of MAP simulations, to test the response of our
method to the inhomogenous sky coverage/noise of MAP. We show that the
resulting $C_\ell$'s are very close to the theoretical expectations. The
HEALPix based implementation of the method, SpICE (Spatially Inhomogenous
Correlation Estimator) will be available to the public from the authors.
